##############################
# Penalised Spline Smoothing #
##############################

library(mgcv)
library(MASS)
library(nlme)
require(splines)
source("fcts.R")

# sample size
n <- 100

# nonlinearity parameter d
d <- 0.1

# covariate, function and response
x<-seq(0,1,1/(n-1))
f <- 1 + x + 2*d*(0.3-x)^2
y<-f+rnorm(n)

# fit penalised spline based on mixed model representation and REML
m1.reml <- gamm(y~s(x,bs="ps", m=c(4,2)), method="REML")

# extract quantities required for computing the AIC
res <- extract.lmeDesign(m1.reml$lme)
X <- res$X
Z <- res$Z
y <- res$y
lambda <- res$lambda
sigmasq <- res$sigmasq
tausq <- sigmasq*lambda

V <- diag(rep(sigmasq, n)) + Z%*%diag(rep(tausq, ncol(Z)))%*%t(Z)
H <- cbind(X,Z)
H <- t(H)%*%H
H1 <- solve(H + diag(c(rep(0,ncol(X)), rep(1/lambda, ncol(Z)))))

beta <- as.numeric(coef(m1.reml$lme)[1:ncol(X)])
b <- as.numeric(coef(m1.reml$lme)[(ncol(X)+1) : (ncol(X)+ncol(Z))])

# conditional and marginal log-likelihood
cllike1reml <- -0.5*n*log(2*pi)-0.5*n*log(sigmasq) -
                         0.5*sum((y - X%*%beta - Z%*%b)^2)/sigmasq
mllike1reml <- -0.5*(n-ncol(X))*log(2*pi)-0.5*sum(log(eigen(V)$values)) -
                           0.5*t(y-X%*%beta)%*%solve(V)%*%(y-X%*%beta) -
                           0.5*sum(log(eigen(t(X)%*%solve(V)%*%X)$values))

# degrees of freedom
cdf1reml <- sum(diag(H1%*%H))
ldf <- liangdf.gamm(y, predict(m1.reml$gam), x, sigma0sq=sigmasq, sigmasq=sigmasq)
cl0df1reml <- ldf[1]
cl1df1reml <- ldf[1]+ldf[2]/2
mdf1reml <- ncol(X)+1+1

# corresponding penalties for the AICs
cpen1asyreml <- cdf1reml+1
cl0pen1asyreml <- cl0df1reml+1
cl1pen1asyreml <- cl1df1reml
mpen1asyreml <- mdf1reml

# AICs
cAIC1asyreml <- -2*cllike1reml + 2*cpen1asyreml
cl0AIC1asyreml <- -2*cllike1reml + 2*cl0pen1asyreml
cl1AIC1asyreml <- -2*cllike1reml + 2*cl1pen1asyreml
mAIC1asyreml <- -2*mllike1reml + 2*mdf1reml


# fit penalised spline based on mixed model representation and ML
m1.ml <- gamm(y~s(x,bs="ps", m=c(4,2)), method="ML")

# extract quantities required for computing the AIC
res <- extract.lmeDesign(m1.ml$lme)
X <- res$X
Z <- res$Z
y <- res$y
lambda <- res$lambda
sigmasq <- res$sigmasq
tausq <- sigmasq*lambda

V <- diag(rep(sigmasq, n)) + Z%*%diag(rep(tausq, ncol(Z)))%*%t(Z)
H <- cbind(X,Z)
H <- t(H)%*%H
H1 <- solve(H + diag(c(rep(0,ncol(X)), rep(1/lambda, ncol(Z)))))

beta <- as.numeric(coef(m1.ml$lme)[1:ncol(X)])
b <- as.numeric(coef(m1.ml$lme)[(ncol(X)+1) : (ncol(X)+ncol(Z))])

# conditional and marginal log-likelihood
cllike1ml <- -0.5*n*log(2*pi)-0.5*n*log(sigmasq) -
                         0.5*sum((y - X%*%beta - Z%*%b)^2)/sigmasq
mllike1ml <- -0.5*n*log(2*pi)-0.5*sum(log(eigen(V)$values)) -
                         0.5*t(y-X%*%beta)%*%solve(V)%*%(y-X%*%beta)

# degrees of freedom
cdf1ml <- sum(diag(H1%*%H))
ldf <- liangdf.gamm(y, predict(m1.ml$gam), x, method="ML", sigma0sq=sigmasq, sigmasq=sigmasq)
cl0df1ml <- ldf[1]
cl1df1ml <- ldf[1]+ldf[2]/2
mdf1ml <- ncol(X)+1+1

# corresponding penalties for the AICs
cpen1asyml <- cdf1ml+1
cl0pen1asyml <- cl0df1ml+1
cl1pen1asyml <- cl1df1ml
mpen1asyml <- mdf1ml

# AICs
cAIC1asyml <- -2*cllike1ml + 2*cpen1asyml
cl0AIC1asyml <- -2*cllike1ml + 2*cl0pen1asyml
cl1AIC1asyml <- -2*cllike1ml + 2*cl1pen1asyml
mAIC1asyml <- -2*mllike1ml + 2*mdf1ml


# Fit the corresponding linear model
m0 <- lm(y~X[,2], x=TRUE, y=TRUE)

# extract quantities required for computing the AIC
X <- m0$x
y <- m0$y
beta <- as.numeric(coef(m0))
sigmasq.ml <- sum((residuals(m0)^2))/n
sigmasq.reml <- sum((residuals(m0)^2))/(n-ncol(X))

# log-likelihood and restricted log-likelihood
llike0ml <- -0.5*n*log(2*pi)-0.5*n*log(sigmasq.ml) -
                          0.5*sum((y - X%*%beta)^2)/sigmasq.ml
llike0reml <- -0.5*(n-ncol(X))*log(2*pi)-0.5*(n-ncol(X))*log(sigmasq.reml) -
                          0.5*sum((y - X%*%beta)^2)/sigmasq.reml -
                          0.5*sum(log(eigen(t(X)%*%X)$values))
# degrees of freedom

df0 <- ncol(X)+1

# corresponding penalties for the AICs
pen0asyml <- df0
pen0asyreml <- df0

# AICs
AIC0asyml <- -2*llike0ml + 2*pen0asyml
AIC0asyreml <- -2*llike0reml + 2*pen0asyreml