# Compute degrees of freedom according to Liang et al. for penalised spline smoothing
liangdf.gamm <- function(y, yhat, x, method="REML", step=0.0001, sigma0sq, sigmasq)
{
n <- length(y)
ldf1 <- 0
ldf2 <- 0
for(i in 1:n)
  {
  newy1 <- y
  mystep1 <- step
  newy1[i] <- y[i]+mystep1
  try(m1 <- gamm(newy1 ~ s(x,bs="ps", m=c(4,2)), method=method))
  while( !exists("m1") || class(m1)=="try-error" ){
    mystep1 <- mystep1 + step
    newy1[i] <- y[i]+mystep1
    try(m1 <- gamm(newy1 ~ s(x,bs="ps", m=c(4,2)), method=method))
    }   
  newyhat1 <- predict(m1$gam)
  sigmasq1 <- m1$lme$sigma^2
  
  newy2 <- y
  mystep2 <- step
  newy2[i] <- y[i]-mystep2
  try(m2 <- gamm(newy2 ~ s(x,bs="ps", m=c(4,2)), method=method))
  while( !exists("m2") || class(m2)=="try-error" ){
    mystep2 <- mystep2 + step
    newy2[i] <- y[i]-mystep2
    try(m2 <- gamm(newy2 ~ s(x,bs="ps", m=c(4,2)), method=method))
    }   
  newyhat2 <- predict(m2$gam)
  sigmasq2 <- m2$lme$sigma^2

  ldf1 <- ldf1 + (newyhat1[i]-yhat[i])/mystep1

  ldf2 <- ldf2 + 2*sigma0sq*(yhat[i]-y[i])*(1/sigmasq1-1/sigmasq)/mystep1 + 
          sigma0sq^2*((1/sigmasq1-1/sigmasq)/mystep1 - (1/sigmasq-1/sigmasq2)/mystep2)/((mystep1+mystep2)/2)
  }
return(c(ldf1,ldf2))
}


# Compute degrees of freedom according to Liang et al. for random effects
liangdf.lme <- function(y, yhat, id, method="REML", step=0.0001, sigma0sq, sigmasq)
{
n <- length(y)
ldf1 <- 0
ldf2 <- 0
for(i in 1:n)
  {
  newy1 <- y
  mystep1 <- step
  newy1[i] <- y[i]+mystep1
  
  try(m1 <- lme(newy1~1, random= ~1|id, method=method))
  while( !exists("m1") || class(m1)=="try-error" ){
    mystep1 <- mystep1 + step
    newy1[i] <- y[i]+mystep1
    try(m1 <- lme(newy1~1, random= ~1|id, method=method))
    }   
  newyhat1 <- predict(m1, level=1)
  sigmasq1 <- m1$sigma^2
  
  newy2 <- y
  mystep2 <- step
  newy2[i] <- y[i]-mystep2
  try(m2 <- lme(newy2~1, random= ~1|id, method=method))
  while( !exists("m2") || class(m2)=="try-error" ){
    mystep2 <- mystep2 + step
    newy2[i] <- y[i]-mystep2
    try(m2 <- lme(newy2~1, random= ~1|id, method=method))
    }   
  newyhat2 <- predict(m2, level=1)
  sigmasq2 <- m2$sigma^2

  ldf1 <- ldf1 + (newyhat1[i]-yhat[i])/mystep1

  ldf2 <- ldf2 + 2*sigma0sq*(yhat[i]-y[i])*(1/sigmasq1-1/sigmasq)/mystep1 + 
          sigma0sq^2*((1/sigmasq1-1/sigmasq)/mystep1 - (1/sigmasq-1/sigmasq2)/mystep2)/((mystep1+mystep2)/2)
  }
return(c(ldf1,ldf2))
}

# adding "p-spline" classes and methods (taken from mgcv help files)

smooth.construct.ps.smooth.spec<-function(object,data,knots)
# a p-spline constructor method function
{ require(splines)
  if (length(object$p.order)==1) m<-rep(object$p.order,2)
  else m<-object$p.order  # m[1] - basis order, m[2] - penalty order
  nk<-object$bs.dim-m[1]  # number of interior knots
  if (nk<=0) stop("basis dimension too small for b-spline order")
  x <- get.var(object$term,data)  # find the data
  xl<-min(x);xu<-max(x);xr<-xu-xl # data limits and range
  xl<-xl-xr*0.001;xu<-xu+xr*0.001;dx<-(xu-xl)/(nk-1)
  if (!is.null(knots)) k <- get.var(object$term,knots)
  else k<-NULL
  if (is.null(k))
  k<-seq(min(x)-dx*(m[1]+1),max(x)+dx*(m[1]+1),length=nk+2*m[1]+2)
  if (length(k)!=nk+2*m[1]+2)
  stop(paste("there should be ",nk+2*m[1]+2," supplied knots"))
  object$X<-spline.des(k,x,m[1]+2,x*0)$design # get model matrix
  if (!object$fixed)
  { S<-diag(object$bs.dim);if (m[2]) for (i in 1:m[2]) S<-diff(S)
    object$S<-list(t(S)%*%S)  # get penalty
    object$S[[1]] <- (object$S[[1]]+t(object$S[[1]]))/2 # exact symmetry
  }
  object$rank<-object$bs.dim-m[2]  # penalty rank
  object$null.space.dim <- m[2]  # dimension of unpenalized space
  object$knots<-k;object$m<-m      # store p-spline specific info.
  object$C<-matrix(colSums(object$X),1,object$bs.dim) #constraint
  object$df<-ncol(object$X)-1      # maximum DoF
  if (object$by!="NA")  # deal with "by" variable
  { by <- get.var(object$by,data) # find by variable
    if (is.null(by)) stop("Can't find by variable")
    object$X<-as.numeric(by)*object$X # form diag(by)%*%X
  }
  class(object)<-"pspline.smooth"  # Give object a class
  object
}

Predict.matrix.pspline.smooth<-function(object,data)
# prediction method function for the p.spline smooth class
{ require(splines)
  x <- get.var(object$term,data)
  X <- spline.des(object$knots,x,object$m[1]+2,x*0)$design
  if (object$by != "NA") { ## handle any by variables
        by <- get.var(object$by, data)
        if (is.null(by))
            stop("Can't find by variable")
        X <- as.numeric(by) * X
  }
  X
}

# extract design matrices, smoothing parameters, etc. from objects returned
# by gamm (taken from package RLRsim)
extract.lmeDesign <-
function(m)
{
    require(mgcv)
    start.level = 1
    data<-m$data
    grps <- getGroups(m)
    n <- length(grps)
    if (is.null(m$modelStruct$varStruct))
        w <- rep(m$sigma, n)
    else {
        w <- 1/varWeights(m$modelStruct$varStruct)
        group.name <- names(m$groups)
        order.txt <- paste("ind<-order(data[[\"", group.name[1],
            "\"]]", sep = "")
        if (length(m$groups) > 1)
            for (i in 2:length(m$groups)) order.txt <- paste(order.txt,
                ",data[[\"", group.name[i], "\"]]", sep = "")
        order.txt <- paste(order.txt, ")")
        eval(parse(text = order.txt))
        w[ind] <- w
        w <- w * m$sigma
    }
    if (is.null(m$modelStruct$corStruct))
        V <- diag(n)
    else {
        c.m <- corMatrix(m$modelStruct$corStruct)
        if (!is.list(c.m))
            V <- c.m
        else {
            V <- matrix(0, n, n)
            gr.name <- names(c.m)
            n.g <- length(c.m)
            j0 <- 1
            ind <- ii <- 1:n
            for (i in 1:n.g) {
                j1 <- j0 + nrow(c.m[[i]]) - 1
                V[j0:j1, j0:j1] <- c.m[[i]]
                ind[j0:j1] <- ii[grps == gr.name[i]]
                j0 <- j1 + 1
            }
            V[ind, ] <- V
            V[, ind] <- V
        }
    }
    V <- as.vector(w) * t(as.vector(w) * V)
    X <- list()
    grp.dims <- m$dims$ncol
    Zt <- model.matrix(m$modelStruct$reStruct, data)
    cov <- as.matrix(m$modelStruct$reStruct)
    i.col <- 1
    n.levels <- length(m$groups)
    Z <- matrix(0, n, 0)
    if (start.level <= n.levels) {
        for (i in 1:(n.levels - start.level + 1)) {
            if(length(levels(m$groups[[n.levels-i+1]]))!=1)
            {
            X[[1]] <- model.matrix(~m$groups[[n.levels - i +
                1]] - 1, contrasts.arg = c("contr.treatment",
                "contr.treatment"))
            }
            else X[[1]]<-matrix(1)
            X[[2]] <- as.matrix(Zt[, i.col:(i.col + grp.dims[i] -
                1)])
            i.col <- i.col + grp.dims[i]
            Z <- cbind(tensor.prod.model.matrix(X),Z)
        }
        Vr <- matrix(0, ncol(Z), ncol(Z))
        start <- 1
        for (i in 1:(n.levels - start.level + 1)) {
            k <- n.levels - i + 1
            for (j in 1:m$dims$ngrps[i]) {
                stop <- start + ncol(cov[[k]]) - 1
                 Vr[ncol(Z)+1-(stop:start),ncol(Z)+1-(stop:start)] <- cov[[k]]
                 start <- stop + 1
            }
        }

        Vlambda <- V/m$sigma^2 + Z %*% Vr %*% t(Z)
    }
X<-model.matrix(formula(m$call$fixed),data)
y<-as.vector(matrix(m$residuals,nc=NCOL(m$residuals))[,NCOL(m$residuals)] +matrix(m$fitted,nc=NCOL(m$fitted))[,NCOL(m$fitted)])
return(list(
             Vlambda=Vlambda, #Cov(y)/Var(Error)
             V=V, #Cov(Error)
             Vr=Vr, #Cov(RanEf)/Var(Error)
             X=X,
             Z=Z,
             sigmasq=m$sigma^2,
             lambda=unique(diag(Vr)),
             y=y
           )
      )
}