library(fda)

##
## Slides 5-6:  fRegress.numeric:  Scalar Response
##
# Ramsay, Hooker, Graves, sec. 9.4.1

logAnnPrecip   = log10(apply(daily$precav,2,sum))

tempbasis65  = create.fourier.basis(c(0,365),65)
tempSmooth65 = smooth.basis(day.5, daily$tempav, tempbasis65)
tempfd65     = tempSmooth65$fd

templist      = vector("list",2)
templist[[1]] = rep(1,35)
templist[[2]] = tempfd65

conbasis   = create.constant.basis(c(0,365))
betabasis5 = create.fourier.basis(c(0,365),5)
betalist1  = vector("list",2)
betalist1[[1]] = conbasis
betalist1[[2]] = betabasis5

fRegressList1 = fRegress(logAnnPrecip, templist, betalist1)

betaestlist1  = fRegressList1$betaestlist
tempbetafd1   = betaestlist1[[2]]$fd

# Figure 9.1

plot(tempbetafd1, xlab="Day", ylab="Beta for temperature")

plot(tempbetafd1, xlab="", ylab="Beta for temperature", axes=FALSE,
     cex.lab=1.5)
axis(2, 0)
axisIntervals(cex.axis=1.5)
text(50, .0012, expression(hat(beta)(t)), cex=4)

##
## Slide 7.  fRegress.numeric: functional response, x = scalar
##

#  Section 10.1.1 Climate Region Effects on Temperature

(regions.         = unique(CanadianWeather$region)[c(3, 2, 4, 1)])
regions     = c("Canada", regions.)
p                = length(regions)
regionList       = vector("list", p)
names(regionList)= regions
regionList[[1]]  = c(rep(1,35),0)
# 35 locations plus a constraint
# that the region effects sum to zero

for (j in 2:p) {
  xj             = (CanadianWeather$region == regions.[j-1])
  regionList[[j]]= c(xj,1)
}

# tempfd from chapter 9

Lcoef       = c(0,(2*pi/365)^2,0)
harmaccelLfd= vec2Lfd(Lcoef, c(0,365))
tempbasis   = create.fourier.basis(c(0, 365), 65)
lambda      = 1e6
tempfdPar65 = fdPar(tempbasis, harmaccelLfd, lambda)
tempShifted = daily$tempav[dayOfYearShifted, ]
tempSmooth65= smooth.basis(day.5, tempShifted, tempfdPar65)
tempfd      = tempSmooth65$fd

coef    = tempfd$coef
coef36  = cbind(coef,matrix(0,65,1))
temp36fd= fd(coef36,tempbasis,tempfd$fdnames)

betabasis      = create.fourier.basis(c(0, 365), 11)
betafdPar      = fdPar(betabasis)
betaList       = vector("list",p)
names(betaList)= regions
for (j in 1:p) betaList[[j]] = betafdPar

fRegressList= fRegress(temp36fd, regionList, betaList)
betaestList = fRegressList$betaestlist
regionFit   = fRegressList$yhatfd

# Figure 10.1

op          = par(mfrow=c(2,3),cex=1)
for (j in 1:p) plot(betaestList[[j]]$fd, lwd=2,
                    xlab="Day (July 1 to June 30)",
                    ylab="", main=regions[j])
plot(regionFit, lwd=2, col=1, lty=1,
     xlab="Day (July 1 to June 30)", ylab="", main="Prediction")
par(op)

dayFine <- 1:365
beta.rgn <- lapply(betaestList[-1], predict, newdata=dayFine)
str(beta.rgn)

rgns <- ordered(regions., levels=regions.)
predTemp <- data.frame(day=dayFine,
                       region=rep(rgns, each=365),
                       beta.region=unlist(beta.rgn))
str(predTemp)
sapply(predTemp, length)

library(lattice)
xyplot(beta.region~day | region, predTemp, layout=c(4, 1), type='l')

xyplot(beta.region~day | region, predTemp, layout=c(4, 1), type='l',
       scales=list(x=list(at=seq(1, 7, 2)*365/8,
           labels=c("Q1", "Q2", "Q3", "Q4") )),
       xlab='')

##
## Slide 8.  fRegress.fdPar: Concurrent Functional Model
##

# Section 10.2.3 Knee Angle Predicted from Hip Angle
gaittime = seq(0.5,19.5,1)
gaitrange = c(0,20)
gaitfine = seq(0,20,len=101)

harmaccelLfd20 = vec2Lfd(c(0, (2*pi/20)^2, 0), rangeval=gaitrange)
gaitbasis = create.fourier.basis(gaitrange, nbasis=21)

#    GCV is minimized with lambda = 10^(-1.5).

gaitSmooth = smooth.basisPar(gaittime, gait,
       gaitbasis, Lfdobj=harmaccelLfd20, lambda=10^(-1.5))
gaitfd = gaitSmooth$fd

names(gaitfd$fdnames) = c("Normalized time", "Child", "Angle")
gaitfd$fdnames[[3]] = c("Hip", "Knee")

hipfd  = gaitfd[,1]
kneefd = gaitfd[,2]

kneefdMean = mean(kneefd)

xfdlist   = list(const=rep(1,39), hip=hipfd)
betafdPar = fdPar(gaitbasis, harmaccelLfd20)
betalist  = list(const=betafdPar, hip=betafdPar)

gaitRegress= fRegress(kneefd, xfdlist, betalist)

# Figure 10.7

op = par(mfrow=c(2,1))

# Intercept
betaestlist = gaitRegress$betaestlist
kneeIntercept = predict(betaestlist$const$fd, gaitfine)

# mean knee angle
kneeMean = predict(kneefdMean, gaitfine)

# Hip coefficient
hipCoef = predict(betaestlist$hip$fd, gaitfine)

# Squared multiple correlation
kneehatfd = gaitRegress$yhatfd$fd
#kneemat = predict(kneefd, gaittime)
kneehatmat = eval.fd(gaittime, kneehatfd)
#resmat. = kneemat - kneehatmat
resmat. = gait[,,'Knee Angle'] - kneehatmat
SigmaE = cov(t(resmat.))

kneefinemat   = eval.fd(gaitfine, kneefd)
kneemeanvec   = eval.fd(gaitfine, mean(kneefd))
kneehatfinemat= eval.fd(gaitfine, kneehatfd)
resmat        = kneefinemat - kneehatfinemat
ncurve        = dim(gait)[2]
resmat0 = kneefinemat - kneemeanvec %*% matrix(1,1,ncurve)
SSE0 = apply((resmat0)^2, 1, sum)
SSE1 = apply(resmat^2, 1, sum)
knee.R2 = (SSE0-SSE1)/SSE0

# Plot Hip Coefficient & Squared Multiple Correlation

ylim2=c(0, max(hipCoef, knee.R2))
plot(gaitfine, hipCoef, lwd=2, xlab='', ylab='', ylim=ylim2, type='l',
     main='Hip Coefficient and Squared Multiple Correlation')
abline(v=c(7.5, 14.7), lty='dashed')
lines(gaitfine, knee.R2, lty='dashed')

# done
par(op)

# Plot Hip Coefficient & Squared Multiple Correlation
# By itself for presentation

ylim2=c(0, max(hipCoef, knee.R2))
plot(gaitfine, hipCoef, lwd=2, xlab='gait cycle', ylab='',
     ylim=ylim2, type='l',
     main='Hip Coefficient and Squared Multiple Correlation')
abline(v=c(7.5, 14.7), lty='dashed')
lines(gaitfine, knee.R2, lty='dashed')
