library(fda)
library(ggplot2)


################################################################################
###                      FUNCTIONAL RESPONSE MODELS                         ####
################################################################################

# We could predict total precipitation from temperature. How about the 
# annual precipitation profile?

# We can also look at constant predictors -- see the weather demo for an ANOVA
# between different weather regions. 

#### Weather by Concurrent Model

data(CanadianWeather)

temp <- CanadianWeather$dailyAv[,,1]
precip <- CanadianWeather$dailyAv[,,3]  # logarithm

daytime <- (1:365)-0.5
day5 = seq(0,365,5)

dayrng <- c(0,365)

# Temperature smoothing

tempbasis <- create.fourier.basis(dayrng, 21)
tempSmooth <- smooth.basis(daytime, temp, tempbasis)
tempfd <- tempSmooth$fd

# log precipitation smoothing

precbasis <-  create.fourier.basis(dayrng, 11)
precSmooth <- smooth.basis(daytime, precip, precbasis)
precfd <- precSmooth$fd

# Smoothed Precipitation
x11()
plot(precfd)

# ggplot
# df <- data.frame(daytime = daytime, temp = c(temp), precip = c(precip),
#                  place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365))
# 
# df$precipfda <- c(eval.fd(precfd, evalarg = daytime))
# 
# x11()
# p0 <- ggplot (data = df, aes (x = daytime, y = precipfda, colour = place))
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = "Log10 of Precipitation") 
# p0 <- p0 + theme_bw ()
# print (p0)

# regression
prec.conc <- fRegress(precfd ~ tempfd)
yhatfd  <- prec.conc$yhatfdobj$fd  # Fitted smooths. 

# # Let's have a look at what we've got

# We can also look at a comparison between predicted and observed
# in Churchill and Regina

x11()
par(mfrow=c(1, 2))
plot(yhatfd[19], ylim = c(-0.45, 0.35))
lines(precfd[19], col = 2)

plot(yhatfd[20], ylim = c(-0.45, 0.35))
lines(precfd[20], col = 2)

# all
plot(yhatfd)
plot(precfd)

# ggplot
# Plot the estimate in Churchill and Regina
# df$preciphat <- c(eval.fd(yhatfd, evalarg = daytime))
# 
# x11()
# p0 <- ggplot (data = subset(df, is.element(place, c("Churchill", "Regina"))), aes (x = daytime, y = precipfda, group = place))
# p0 <- p0 + geom_line (color = "grey", size = 0.5)
# p0 <- p0 + geom_line (aes(y = preciphat, color = place))
# p0 <- p0 + labs (x = "Days", y = "Log10 of Precipitation")
# p0 <- p0 + facet_wrap( ~ place)
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)

# plot parameters 
x11()
par(mfrow=c(2,1))
plot(prec.conc$betaestlist[[1]])
plot(prec.conc$betaestlist[[2]])


# And compare observed with residuals
plot(precfd-yhatfd)
plot(precfd)
# 
# # Doesn't look like we really got very much. 


#### 2. Confidence Intervals

## In order to obtain confidence intervals, we can include the smoothing of 
## precipitation as a source of error. In order to do this, we need two things

y2cmap <- precSmooth$y2cMap

# This is the matrix that goes from the observations to the coefficients. 
# We now need a covariance matrix for the errors of the original observed 
# precipitation from the functional linear model

Errmat <- eval.fd(daytime, precfd) - eval.fd(daytime, yhatfd) 
SigmaE2 <- cov(t(Errmat))

# We can now run fRegress.stderr

conc.std <- fRegress.stderr(prec.conc, y2cmap, SigmaE2)

# And plot the results

plotbeta(prec.conc$betaestlist,conc.std$betastderrlist)

# There really doesn't appear to be much going on. 

# ggplot
# beta0 <- eval.fd(prec.conc$betaestlist[[1]]$fd, evalarg = daytime)
# beta1 <- eval.fd(prec.conc$betaestlist[[2]]$fd, evalarg = daytime)
# 
# betastderr0 <- eval.fd(conc.std$betastderrlist[[1]], evalarg = daytime)
# betastderr1 <- eval.fd(conc.std$betastderrlist[[2]], evalarg = daytime)
# 
# betadf <- data.frame(x = c(daytime, daytime), beta = c(beta0, beta1), 
#                      betastderr = c(betastderr0, betastderr1), 
#                      type = factor(rep(c("Intercept", "Slope"), each = length(daytime)), ordered = TRUE))
# 
# x11()
# p0 <- ggplot (data = betadf, aes (x = x, y = beta, group = type))
# p0 <- p0 + geom_ribbon (aes (ymin = beta - 1.96*betastderr, ymax = beta + 1.96*betastderr), fill = "grey")
# p0 <- p0 + geom_line (size = 1.5)
# p0 <- p0 + geom_hline (yintercept = 0, linetype = 2)
# p0 <- p0 + labs (x = "Days", y = "Beta(t)")
# p0 <- p0 + facet_wrap( ~ type, scales = "free_y")
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)


# prediction with CI

# bootstrap
library(MASS)
yhat <- eval.fd(daytime, yhatfd)
q.025 <- matrix(NA, nrow = nrow(yhat), ncol = ncol(yhat))
q.975 <- matrix(NA, nrow = nrow(yhat), ncol = ncol(yhat))

for (k in 1:35){
samp1000 <- mvrnorm(1000, mu = yhat[, k], Sigma = SigmaE2)
quants <- apply(samp1000, 2, quantile, probs = c(0.025, 0.975),  na.rm = TRUE)
q.025[, k] <- quants[1, ]
q.975[, k] <- quants[2, ]
}

# plot estimates in Churchill and Regina with CI
x11()
par(mfrow=c(1, 2))
plot(yhatfd[19], ylim = c(-0.45, 0.35))
x <- c(daytime, rev(daytime))
y <- c(q.025[, 19], rev(q.975[, 19]))
polygon(x, y, col="gray75", border="gray75")
lines(yhatfd[19])
lines(precfd[19], col = 2)

plot(yhatfd[20], ylim = c(-0.45, 0.35))
x <- c(daytime, rev(daytime))
y <- c(q.025[, 20], rev(q.975[, 20]))
polygon(x, y, col="gray75", border="gray75")
lines(yhatfd[20])
lines(precfd[20], col = 2)

# ggplot

# df$lower.pred <- c(q.025)
# df$upper.pred <- c(q.975)
# 
# x11()
# p0 <- ggplot (data = subset(df, is.element(place, c("Churchill", "Regina"))), aes (x = daytime, y = precipfda, group = place))
# p0 <- p0 + geom_ribbon (aes (ymin = lower.pred, ymax = upper.pred, fill = place), alpha = 0.3)
# p0 <- p0 + geom_line (color = "grey", size = 0.5)
# p0 <- p0 + geom_line (aes(y = preciphat, color = place))
# p0 <- p0 + labs (x = "Days", y = "Log10 of Precipitation")
# p0 <- p0 + facet_wrap( ~ place)
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)


# Assesing the fit
# plot residuals
x11()
plot(precfd-yhatfd)

# ggplot
# df.res <- data.frame(daytime = daytime, residuals = c(Errmat),
#                      place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365))
# 
# x11()
# p0 <- ggplot (data = df.res, aes (x = daytime, y = residuals, colour = place))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = "Residuals") 
# p0 <- p0 + theme_bw ()
# print (p0)

# Functional R^2
nomin <- mean.fd((precfd - yhatfd)^2)
denom <- mean.fd((center.fd(precfd))^2)

Rvalues <- 1- eval.fd(daytime, nomin)/eval.fd(daytime, denom)

plot(daytime, Rvalues, type = "l", ylim = c(0, 0.75))

# ggplot
# df.r2 <- data.frame(daytime = daytime, r2 = c(Rvalues))
# 
# x11()
# p0 <- ggplot (data = df.r2, aes (x = daytime, y = r2))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = expression(R^2)) 
# p0 <- p0 + theme_bw ()
# print (p0)


# F-Ratio

F.res = Fperm.fd(prec.conc$yfdPar$fd, prec.conc$xfdlist, prec.conc$betalist)

# ggplot

# popisky <- c("Observed statistic", "pointwise 0.05 critical value", "maximum 0.05 critical value")
# df.Fratio <- data.frame(daytime = F.res$argvals, Fratio = c(F.res$Fvals, F.res$qvals.pts, rep(F.res$qval, length(F.res$argvals))), 
#                         type = factor(rep(popisky, each = length(F.res$argvals)), ordered = TRUE, levels = popisky))
# 
# x11()
# p0 <- ggplot (data = df.Fratio, aes (x = daytime, y = Fratio, linetype = type))
# p0 <- p0 + geom_line ()
# p0 <- p0 + theme_bw ()
# p0 <- p0 + labs (x = "Days", y = "F-statistic", linetype = "") 
# print (p0)

#### Weather by Ramsey

# Functional ANOVA

# As the first, set covariates
regions <- unique(CanadianWeather$region)
p <- length(regions) + 1
regionList <- vector("list", p)
regionList[[1]] <- c(rep(1, 35), 0)
for (j in 2:p) {
  xj <- CanadianWeather$region == regions[j-1]
  regionList[[j]] <- c(xj, 1)
}

# The next step is to augment the temperature functional data object by a 36th observation
# that takes only zero values as required by additional condition Sum \beta_j(t) = 0

coef <- tempfd$coef
n.coef <- nrow(coef)
coef36 <- cbind(coef, matrix(0, n.coef, 1))
temp36fd <- fd(coef36, tempbasis, tempfd$fdnames)

# We now create functional parameter objects for each of the coefficient functions,
# using 11 Fourier basis functions for each.

betabasis <- create.fourier.basis(dayrng, 11)
betafdPar <- fdPar(betabasis)
betaList <- vector("list",p)
for (j in 1:p) betaList[[j]] <- betafdPar

# Now call fRegress, extract the coefficients and plot them, along with the predicted
# curves for the regions.

fRegressList <- fRegress(temp36fd, regionList, betaList)
betaestList <- fRegressList$betaestlist
regionFit <- fRegressList$yhatfd
regions <- c("Canada", regions)

TempRes36fd <- temp36fd - regionFit$fd   # residuals
TempRes36 <- eval.fd(daytime, TempRes36fd)

# regression

# Removing 36th observation
coefsres <- TempRes36fd$coefs[, -36]
TempResfd <- fd(coefsres, regionFit$fd$basis, regionFit$fd$fdnames)

coefsfit <- regionFit$fd$coefs[, -36]
regionFit1 <- fd(coefsfit, regionFit$fd$basis, regionFit$fd$fdnames)

prec.conc.mix <- fRegress(precfd ~ regionFit1 + TempResfd)

yhatmixfd  <- prec.conc.mix$yhatfdobj$fd  # Fitted smooths. 

# Comparison between predicted and observed
# in Churchill and Regina

x11()
par(mfrow=c(1, 2))
plot(yhatmixfd[19], ylim = c(-0.45, 0.35))
lines(precfd[19], col = 2)

plot(yhatmixfd[20], ylim = c(-0.45, 0.35))
lines(precfd[20], col = 2)

# ggplot

# dfmix <- df
# dfmix$preciphat <- c(eval.fd(yhatmixfd, evalarg = daytime))
# 
# x11()
# p0 <- ggplot (data = subset(dfmix, is.element(place, c("Churchill", "Regina"))), aes (x = daytime, y = precipfda, group = place))
# p0 <- p0 + geom_line (color = "grey", size = 0.5)
# p0 <- p0 + geom_line (aes(y = preciphat, color = place))
# p0 <- p0 + labs (x = "Days", y = "Log10 of Precipitation")
# p0 <- p0 + facet_wrap( ~ place)
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)

#### 2. Confidence Intervals

y2cmap = precSmooth$y2cMap

Errmat = eval.fd(daytime, precfd) - eval.fd(daytime, yhatmixfd) 
SigmaE2 = cov(t(Errmat))

conc.std.mix = fRegress.stderr(prec.conc.mix, y2cmap, SigmaE2)

# plot parameters - just Beta_2

plotbeta(prec.conc.mix$betaestlist[[2]], conc.std.mix$betastderrlist[[2]])

# ggplot
# beta0 <- eval.fd(prec.conc.mix$betaestlist[[1]]$fd, evalarg = daytime)
# beta1 <- eval.fd(prec.conc.mix$betaestlist[[2]]$fd, evalarg = daytime)
# beta2 <- eval.fd(prec.conc.mix$betaestlist[[3]]$fd, evalarg = daytime)
# 
# betastderr0 <- eval.fd(conc.std.mix$betastderrlist[[1]], evalarg = daytime)
# betastderr1 <- eval.fd(conc.std.mix$betastderrlist[[2]], evalarg = daytime)
# betastderr2 <- eval.fd(conc.std.mix$betastderrlist[[3]], evalarg = daytime)
# 
# betamixdf <- data.frame(x = c(daytime, daytime), beta = c(beta2), 
#                      betastderr = c(betastderr2), 
#                      type = factor(rep(c("Beta"), each = length(daytime)), ordered = TRUE))
# 
# x11()
# p0 <- ggplot (data = betamixdf, aes (x = x, y = beta, group = type))
# p0 <- p0 + geom_ribbon (aes (ymin = beta - 1.96*betastderr, ymax = beta + 1.96*betastderr), fill = "grey")
# p0 <- p0 + geom_line (size = 1.5)
# p0 <- p0 + geom_hline (yintercept = 0, linetype = 2)
# p0 <- p0 + labs (x = "Days", y = "Beta(t)")
# p0 <- p0 + facet_wrap( ~ type, scales = "free_y")
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)

# Assesing the fit
# plot residuals
x11()
plot(precfd-yhatmixfd)

# ggplot
# df.res.mix <- data.frame(daytime = daytime, residuals = c(Errmat),
#                      place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365))
# 
# x11()
# p0 <- ggplot (data = df.res.mix, aes (x = daytime, y = residuals, colour = place))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = "Residuals") 
# p0 <- p0 + theme_bw ()
# print (p0)


# Functional R^2
nomin <- mean.fd((precfd - yhatmixfd)^2)
denom <- mean.fd((center.fd(precfd))^2)

Rvalues <- 1 - eval.fd(daytime, nomin)/eval.fd(daytime, denom)
x11()
plot(daytime, Rvalues, type = "l", ylim = c(0, 0.8))

# ggplot
# df.r2 <- data.frame(daytime = daytime, r2 = c(Rvalues))
# 
# x11()
# p0 <- ggplot (data = df.r2, aes (x = daytime, y = r2))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = expression(R^2)) 
# p0 <- p0 + theme_bw ()
# print (p0)


# F-Ratio

F.res = Fperm.fd(prec.conc.mix$yfdPar$fd, prec.conc.mix$xfdlist, prec.conc.mix$betalist)

# ggplot
# popisky <- c("Observed statistic", "pointwise 0.05 critical value", "maximum 0.05 critical value")
# df.Fratio <- data.frame(daytime = F.res$argvals, Fratio = c(F.res$Fvals, F.res$qvals.pts, rep(F.res$qval, length(F.res$argvals))), 
#                         type = factor(rep(popisky, each = length(F.res$argvals)), ordered = TRUE, levels = popisky))
# 
# x11()
# p0 <- ggplot (data = df.Fratio, aes (x = daytime, y = Fratio, linetype = type))
# p0 <- p0 + geom_line ()
# p0 <- p0 + theme_bw ()
# p0 <- p0 + labs (x = "Days", y = "F-statistic", linetype = "") 
# print (p0)

### Fully functional Model

# Canadian Weather

# Temperature smoothing
tempbasis <- create.fourier.basis(dayrng, 21) #65
tempSmooth <- smooth.basis(daytime, temp, tempbasis)
tempfd <- tempSmooth$fd

# log precipitation smoothing
precbasis <-  create.fourier.basis(dayrng, 11) #65
precSmooth <- smooth.basis(daytime, precip, precbasis)
precfd <- precSmooth$fd

# Set up for the list of regression coefficient fdPar objects
nbasis <- 15 # 31

PrecBetaBasis <- create.fourier.basis(dayrng, nbasis)
PrecBeta0Par <- fdPar(PrecBetaBasis)
PrecBeta1fd  <- bifd(matrix(0, nbasis, nbasis), PrecBetaBasis, PrecBetaBasis)
PrecBeta1Par <- bifdPar(PrecBeta1fd)
PrecBetaList <- list(PrecBeta0Par, PrecBeta1Par)

#  Do the regression analysis
Prec.linmod <- linmod(tempfd, precfd, PrecBetaList)

Prec.beta1mat <- eval.bifd(day5, day5, Prec.linmod$beta1estbifd)
yhatfullfd <- Prec.linmod$yhatfdobj

# Plot Beta
persp(day5, day5, Prec.beta1mat,
      xlab="age", ylab="age",zlab="beta(s,t)",
      cex.lab=1.5,cex.axis=1.5)

# with RGL
# library(rgl)
# open3d()
# persp3d (day5, day5, Prec.beta1mat, ylab = "time s [days]", xlab = "time t [days]", zlab = "beta(t,s)",
#           col = "lightblue")

# Plot the estimate in Churchill and Regina

x11()
par(mfrow=c(1, 2))
plot(yhatfullfd[19], ylim = c(-0.5, 0.35))
lines(precfd[19], col = 2)

plot(yhatfullfd[20], ylim = c(-0.5, 0.35))
lines(precfd[20], col = 2)

# ggplot

# 
# df$preciphatfull <- c(eval.fd(yhatfullfd, evalarg = daytime))
# 
# x11()
# p0 <- ggplot (data = subset(df, is.element(place, c("Churchill", "Regina"))), aes (x = daytime, y = precipfda, group = place))
# p0 <- p0 + geom_line (color = "grey", size = 0.5)
# p0 <- p0 + geom_line (aes(y = preciphatfull, color = place))
# p0 <- p0 + labs (x = "Days", y = "Log10 of Precipitation")
# p0 <- p0 + facet_wrap( ~ place)
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)

# Assesing the fit
y2cmap = precSmooth$y2cMap

Errmat = eval.fd(daytime, precfd) - eval.fd(daytime, yhatfullfd) 
SigmaE2 = cov(t(Errmat))

# plot residuals
x11()
plot(precfd-yhatfullfd)

# ggplot
# dffull.res <- data.frame(daytime = daytime, residuals = c(Errmat),
#                      place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365))
# 
# x11()
# p0 <- ggplot (data = dffull.res, aes (x = daytime, y = residuals, colour = place))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = "Residuals") 
# p0 <- p0 + theme_bw ()
# print (p0)

# Functional R^2
nomin <- mean.fd((precfd - yhatfullfd)^2)
denom <- mean.fd((center.fd(precfd))^2)

Rvalues <- 1 - eval.fd(daytime, nomin)/eval.fd(daytime, denom)

plot(daytime, Rvalues, ylim = c(0, 1), type = "l")

# ggplot
# df.r2 <- data.frame(daytime = daytime, r2 = c(Rvalues))
# 
# x11()
# p0 <- ggplot (data = df.r2, aes (x = daytime, y = r2))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = expression(R^2)) 
# p0 <- p0 + theme_bw ()
# print (p0)


# F-Ratio
# Fperm.fd cannot be used!

nomin <- sum((center.fd(precfd))^2)
denom <- sum((precfd - yhatfullfd)^2)

Fvalues <- eval.fd(daytime, nomin)/eval.fd(daytime, denom)*35/(35-1)
kvantil <- qf(0.95, df1 = 35 - 1, df2 = length(precfd$fdnames$reps))

x11()
plot(daytime, Fvalues, ylim = c(0, max(Fvalues)), type = "l")
lines(daytime, rep(kvantil, length(daytime)), lty = "dashed")

# ggplot

# df.Fratio <- data.frame(daytime = daytime, Fratio = c(Fvalues))
# 
# x11()
# p0 <- ggplot (data = df.Fratio, aes (x = daytime, y = Fratio))
# p0 <- p0 + geom_line ()
# p0 <- p0 + geom_hline (yintercept = kvantil, linetype = "dashed")
# p0 <- p0 + theme_bw () + ylim(0, max(Fvalues))
# p0 <- p0 + labs (x = "Days", y = "F-statistic") 
# print (p0)


##### Swedish mortality data

#  Human Mortality Database. University of California, Berkeley (USA),
#  and Max Planck Institute for Demographic Research (Germany).

#   Two data objects are required for these analyses:

#  SwedeMat:  a dataframe object with 81 rows and 144 columns
#             containing the log hazard values for ages 0 through 80
#             and years 1751 through 1884
#  Swede1920: a vector object containing log hazard values for 1914

load("SweMortality.RData")

SwedeMat <- log(Sweden$y[1:81,])

Swede1914 <- SwedeMat[,164]
SwedeLogHazard <- SwedeMat[,1:144]

# SwedeLogHazard = as.matrix(SwedeMat)

dimnames(SwedeLogHazard)[[2]] <- paste('b', 1751:1894, sep='')

# Data for years 1751, 1810, 1860, 1914

select_data <- cbind(SwedeLogHazard[, c('b1751', 'b1810', 'b1860')],
                     Swede1914)
SwedeTime <- 0:80;
matplot(SwedeTime, select_data,
        type = 'l', lwd = 2, xlab = 'age', ylab = 'log Hazard', col=1)

# ggplot
# df.mort <- data.frame(Age = 0:80, Hazard = c(c(SwedeLogHazard),Swede1914), Year = factor(rep(c(1751:1894,1914), each = 81)))
# 
# df.vyb <- subset(df.mort, is.element(Year, c("1751", "1810", "1860", "1914")))
# 
# x11()
# p0 <- ggplot (data = df.vyb, aes (x = Age, y = Hazard, colour = Year))
# p0 <- p0 + geom_line ()
# p0 <- p0 + theme_bw ()
# p0 <- p0 + labs (x = "Age", y = "log Hazard") 
# print (p0)


#  smooth the log hazard observations

SwedeTime <- 0:80
SwedeRng <- c(0,80)

nbasis <- 85
norder <- 6
SwedeBasis <- create.bspline.basis(SwedeRng, nbasis, norder)

D2fdPar <- fdPar(SwedeBasis, lambda = 1e-7)

SwedeLogHazfd <- smooth.basis(SwedeTime, SwedeLogHazard, D2fdPar)$fd
SwedeLogHaz1914fd <- smooth.basis(SwedeTime, Swede1914, D2fdPar)$fd

# The following requires manually clicking on the plot
# for each of 144 birth year cohorts

# plotfit.fd(SwedeLogHazard,SwedeTime,SwedeLogHazfd)
plot(SwedeLogHazfd)

# ggplot

# df.mort.sm <- data.frame(Age = 0:80, Hazard = c(c(eval.fd(SwedeTime, SwedeLogHazfd)), c(eval.fd(SwedeTime, SwedeLogHaz1914fd))), 
#                          Year = factor(rep(c(1751:1894, 1914), each = 81)))
# 
# df.vyb.sm <- subset(df.mort.sm, is.element(Year, c("1751", "1810", "1860", "1914")))
# 
# x11()
# p0 <- ggplot (data = df.vyb.sm, aes (x = Age, y = Hazard, colour = Year))
# p0 <- p0 + geom_line ()
# p0 <- p0 + theme_bw ()
# p0 <- p0 + labs (x = "Age", y = "log Hazard") 
# print (p0)


# Set up for the list of regression coefficient fdPar objects

nbasis <- 23
SwedeRng <- c(0,80)
SwedeBetaBasis <- create.bspline.basis(SwedeRng,nbasis)

SwedeBeta0Par <- fdPar(SwedeBetaBasis, 2, 1e-5)

SwedeBeta1fd  <- bifd(matrix(0,23,23), SwedeBetaBasis, SwedeBetaBasis)

SwedeBeta1Par <- bifdPar(SwedeBeta1fd, 2, 2, 1e3, 1e3)

SwedeBetaList <- list(SwedeBeta0Par, SwedeBeta1Par)

#  Define the dependent and independent variable objects

NextYear <- SwedeLogHazfd[2:144]
LastYear <- SwedeLogHazfd[1:143]

#  Do the regression analysis

Swede.linmod <- linmod(NextYear, LastYear, SwedeBetaList)

# Parameters

# Beta 0
# Swede.beta0 <- eval.fd(SwedeTime, Swede.linmod$beta0estfd)
# plot(Swede.linmod$beta0estfd)

# Beta 1
Swede.ages <- seq(0, 80, 2)
Swede.beta1mat <- eval.bifd(Swede.ages, Swede.ages, Swede.linmod$beta1estbifd)

persp(Swede.ages, Swede.ages, Swede.beta1mat,
      xlab="age", ylab="age",zlab="beta(s,t)",
      cex.lab=1.5,cex.axis=1.5)

# library(rgl)
# open3d()
# persp3d (Swede.ages, Swede.ages, Swede.beta1mat, ylab = "Age [yrs]", xlab = "Age t [yrs]", zlab = "beta(t,s)",
#          col = "lightblue")

# Plot the estimate
yhatfullfd <- Swede.linmod$yhatfdobj

ind1 <- which(yhatfullfd$fdnames$reps == "b1810")
ind2 <- which(yhatfullfd$fdnames$reps == "b1860")

x11()
par(mfrow=c(1, 2))
plot(yhatfullfd[ind1])
lines(SwedeLogHazfd[ind1], col = 2)

plot(yhatfullfd[ind2])
lines(SwedeLogHazfd[ind2], col = 2)

# ggplot

# df.mort.sm$preciphatfull <- c(c(Swede1914), c(eval.fd(yhatfullfd, evalarg = SwedeTime)), c(Swede1914))
# 
# x11()
# p0 <- ggplot (data = subset(df.mort.sm, is.element(Year, c("1810", "1860"))), aes (x = Age, y = Hazard, group = Year))
# p0 <- p0 + geom_line (color = "grey", size = 0.5)
# p0 <- p0 + geom_line (aes(y = preciphatfull, color = Year))
# p0 <- p0 + labs (x = "Age", y = "log Hazard")
# p0 <- p0 + facet_wrap( ~ Year)
# p0 <- p0 + guides (colour = FALSE) 
# p0 <- p0 + theme_bw ()
# print (p0)

# Assesing the fit
y2cmap <- SwedeLogHazfd$y2cMap

Errmat <- eval.fd(SwedeTime, NextYear) - eval.fd(SwedeTime, yhatfullfd) 
SigmaE2 <- cov(t(Errmat))

# plot residuals
matplot(SwedeTime, Errmat[, c(ind1, ind2)], type = "l")

# ggplot
# dffull.res <- data.frame(Age = SwedeTime, residuals = c(Errmat),
#                          Year = factor(rep(1752:1894, each = 81)))
# 
# x11()
# p0 <- ggplot (data = subset(dffull.res, is.element(Year, c("1810", "1860"))), aes (x = Age, y = residuals, colour = Year))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Days", y = "Residuals") 
# p0 <- p0 + theme_bw ()
# print (p0)

# Functional R^2
nomin <- mean.fd((NextYear - yhatfullfd)^2)
denom <- mean.fd((center.fd(NextYear))^2)

Rvalues <- 1 - eval.fd(SwedeTime[-1], nomin)/eval.fd(SwedeTime[-1], denom)

plot(SwedeTime[-1], Rvalues, ylim = c(0, 1), type = "l")

# ggplot
# df.r2 <- data.frame(Age = SwedeTime[-1], r2 = c(Rvalues))
# 
# x11()
# p0 <- ggplot (data = df.r2, aes (x = Age, y = r2))
# p0 <- p0 + geom_hline (yintercept = 0, linetype = "dashed")
# p0 <- p0 + geom_line ()
# p0 <- p0 + labs (x = "Age", y = expression(R^2)) 
# p0 <- p0 + theme_bw ()
# print (p0)

# F-Ratio

nomin <- sum((center.fd(NextYear))^2)
denom <- sum((NextYear - yhatfullfd)^2)

Fvalues <- eval.fd(SwedeTime[-1], nomin)/eval.fd(SwedeTime[-1], denom)*145/(145-1)

kvantil <- qf(0.95, df1 = 145 - 1, df2 = length(NextYear$fdnames$reps))

x11()

plot(SwedeTime[-1], Fvalues, ylim = c(0, max(Fvalues)), type = "l")
lines(SwedeTime[-1], rep(kvantil, length(SwedeTime[-1])), lty = "dashed")

# ggplot
# df.Fratio <- data.frame(Age = SwedeTime[-1], Fratio = c(Fvalues))
# 
# x11()
# p0 <- ggplot (data = df.Fratio, aes (x = Age, y = Fratio))
# p0 <- p0 + geom_line ()
# p0 <- p0 + geom_hline (yintercept = kvantil, linetype = "dashed")
# p0 <- p0 + theme_bw () + ylim(0, max(Fvalues))
# p0 <- p0 + labs (x = "Age", y = "F-statistic") 
# print (p0)

########################
# Problems
######################

##### GAIT data
gaittime <- as.numeric(dimnames(gait)[[1]])*20
gaitrange <- c(0,20)
gaitbasis <- create.fourier.basis(gaitrange, nbasis = 21)
harmaccelLfd <- vec2Lfd(c(0, (2*pi/20)^2, 0), rangeval = gaitrange)
gaitfd <- smooth.basisPar(gaittime, gait,
                          gaitbasis, Lfdobj=harmaccelLfd, lambda=1e-2)

# plot
x11()
par(mfrow = c(2, 1))
plot(gaitfd)


##########################
#  Simulation
##########################

# As the first, generate Xi(t) as the Gaussian process
library(MASS)
N <- 30
C <- function(x, y) 0.01 * exp(-0.001 * (x - y)^2) # covariance function
M <- function(x) sin(x/365) # mean function
t <- seq(1, 365, length.out = 100) # will sample the GP at these points
k <- length(t)
m <- M(t)
S <- matrix(nrow = k, ncol = k)
for (i in 1:k) for (j in 1:k) S[i, j] = C(t[i], t[j])
EX <- t(mvrnorm(N, m, S))
Xi <- EX + matrix(rnorm(k*N, sd = 0.02), nrow = k)


# Regression Model

Beta0 <- 1 + 2*t/30 - (t/30)^2
rozsah <- range(t)
gr <- merge(t, t)
Beta1 <- matrix(sin(2*pi*(gr$x - gr$y)/365), nrow = k)

EX <- t(EX)
EY <- matrix(0, dim(EX)[1], dim(EX)[2])
Yi <- matrix(0, dim(EX)[1], dim(EX)[2])
for (j in 1:N){
  EY[j, ] <- t(Beta0) + (EX[j,] %*% Beta1)*(rozsah[2] - rozsah[1])/k
  Yi[j, ] <- EY[j,] + rnorm(k, 0, 10)
}
Yi <- t(Yi)
EY <- t(EY)
matplot(Yi, type = "l")
matplot(Xi, type = "l")