R/sme.R

In sme: Smoothing-Splines Mixed-Effects Models

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

# Copyright - Maurice Berk (maurice@mauriceberk.com) (http://www2.imperial.ac.uk/~mab201)
#
# block diagonal matrix code thanks to http://tolstoy.newcastle.edu.au/R/help/04/05/1320.html

sme <- function(
  object,
  tme,
  ind,
  verbose=F,
  lambda.mu=NULL,
  lambda.v=NULL,
  maxIter=500,
  knots=NULL,
  zeroIntercept=F,
  deltaEM=1e-3,
  deltaNM=1e-3,
  criteria="AICc",
  initial.lambda.mu=10000,
  initial.lambda.v=10000,
  normalizeTime=FALSE,
  ...)
{
  UseMethod("sme")
}

sme.default <- function(
  object,
  tme,
  ind,
  verbose=F,
  lambda.mu=NULL,
  lambda.v=NULL,
  maxIter=500,
  knots=NULL,
  zeroIntercept=F,
  deltaEM=1e-3,
  deltaNM=1e-3,
  criteria="AICc",
  initial.lambda.mu=10000,
  initial.lambda.v=10000,
  normalizeTime=FALSE,...)
{
  y <- object
  ind <- as.factor(ind)
  ind <- droplevels(ind)
  ind.unique <- sort(unique(ind))
  y <- y[order(ind)]
  tme <- tme[order(ind)]
  ind <- ind[order(ind)]

  call <- match.call(call=sys.call(-1))

  intercept <- min(tme)
  if(normalizeTime)
  {
    tme.orig <- tme
    tme <- (tme - min(tme)) / max(tme)
  }

  if(!is.null(knots))
  {
    if(normalizeTime)
    {
      knots.orig <- knots
      knots <- (knots - min(tme.orig)) / max(tme.orig)
    }
  
    max.tme <- max(tme)
    min.tme <- min(tme)
    B <- ns(c(min.tme,knots,max.tme),knots=knots,intercept=T)

    G <- roughnessMatrix(incidenceMatrix(c(min.tme,knots,max.tme)))
    G <- t(B) %*% G %*% B
  }
  else
  {
    G <- roughnessMatrix(incidenceMatrix(tme))
  }

  if(zeroIntercept)
  {
    y.orig <- y
    tme.orig2 <- tme
    ind.orig <- ind
    
    y <- y[tme!=min(tme)]
    ind <- ind[tme!=min(tme)]
    tme <- tme[tme!=min(tme)]

    if(normalizeTime)
    {
      tme.orig <- tme.orig[tme.orig!=min(tme.orig)]
    }

    G <- G[-1,-1]
  }
  
  yi <- split(y,ind)
  tmei <- split(tme,ind)
  Ni <- sapply(yi,length)

  if(is.null(knots))
  {
    X <- incidenceMatrix(tme)
    Xi <- split.data.frame(X,ind)
  }
  else
  {
    X <- ns(tme,knots=knots,intercept=T)
    Xi <- split.data.frame(X,ind)
  }

  if(is.null(lambda.mu) || is.null(lambda.v))
  {
    res.em <- .C("SMEOptimization",
                 y=as.double(y),
                 tme=as.double(tme),
                 ind=as.integer(ind),
                 X=as.double(X),
                 N=as.integer(length(y)),
                 n=as.integer(length(yi)),
                 Ni=as.integer(sapply(yi,length)),
                 p=as.integer(ncol(X)),
                 lambdaMu=as.double(initial.lambda.mu),
                 lambdaV=as.double(initial.lambda.v),
                 G=as.double(G),
                 mu=double(ncol(X)),
                 sigmaSquared=double(1),
                 D=double(ncol(X) * ncol(X)),
                 v=double(ncol(X) * length(yi)),
                 likelihood=double(1),
                 dfMu=double(1),
                 dfV=double(1),
                 iterations=integer(1),
                 maxIterations=as.integer(maxIter),
                 deltaEM=as.double(deltaEM),
                 deltaNM=as.double(deltaNM),
                 criteria=as.integer(match(criteria,c("AIC","AICc","BICN","BICn"))),
                 verbose=as.integer(verbose),
                 info=integer(1),
                 PACKAGE="sme")
  }
  else
  {
    res.em <- .C("SME",
                 y=as.double(y),
                 tme=as.double(tme),
                 ind=as.integer(ind),
                 X=as.double(X),
                 N=as.integer(length(y)),
                 n=as.integer(length(yi)),
                 Ni=as.integer(sapply(yi,length)),
                 p=as.integer(ncol(X)),
                 lambdaMu=as.double(lambda.mu),
                 lambdaV=as.double(lambda.v),
                 G=as.double(G),
                 mu=double(ncol(X)),
                 sigmaSquared=double(1),
                 D=double(ncol(X) * ncol(X)),
                 v=double(ncol(X) * length(yi)),
                 likelihood=double(1),
                 dfMu=double(1),
                 dfV=double(1),
                 iterations=integer(1),
                 maxIterations=as.integer(maxIter),
                 deltaEM=as.double(deltaEM),
                 verbose=as.integer(verbose),
                 info=integer(1),
                 PACKAGE="sme")
  }

  return.value <- list()
  if(is.null(knots))
  {
    return.value$coefficients <- rbind(as.vector(res.em$mu),t(matrix(res.em$v,ncol=res.em$n)))
    if(normalizeTime)
    {
      if(zeroIntercept)
      {
        colnames(return.value$coefficients) <- sort(unique(tme.orig[tme.orig!=intercept]))
      }
      else
      {
        colnames(return.value$coefficients) <- sort(unique(tme.orig))
      }
    }
    else
    {
      colnames(return.value$coefficients) <- attr(X,"tau")
    }
  }
  else
  {
    return.value$coefficients <- rbind(as.vector(B %*% res.em$mu),do.call(rbind,lapply(split(res.em$v,rep(1:res.em$n,each=res.em$p)),function(v) as.vector(B %*% v))))
    if(normalizeTime)
    {
      colnames(return.value$coefficients) <- c(min(tme.orig),knots.orig,max(tme.orig))
    }
    else
    {
      colnames(return.value$coefficients) <- c(min.tme,knots,max.tme)
    }
  }
  rownames(return.value$coefficients) <- c("mu",paste("v",ind.unique,sep=""))
  return.value$fitted <- unlist(mapply(Xi,split(res.em$v,rep(1:res.em$n,each=res.em$p)),FUN=function(Xi,vi) Xi %*% (res.em$mu + vi),SIMPLIFY=FALSE))
  return.value$residuals <- y - as.vector(return.value$fitted)
  return.value$data <- data.frame(y=y,tme=if(normalizeTime){ tme.orig }else{ tme },ind=ind)
  #return.value$em <- res.em
  return.value$call <- call
  return.value$dfMu <- res.em$dfMu
  return.value$dfV <- res.em$dfV
  return.value$logLik <- res.em$likelihood
  return.value$nobs <- length(y)
  return.value$df <- c(mu=res.em$dfMu,v=res.em$dfV)
  return.value$smoothingParameters <- c(mu=res.em$lambdaMu,v=res.em$lambdaV)
  return.value$parameters <- list(sigmaSquared=res.em$sigmaSquared,D=matrix(res.em$D,nrow=ncol(X)))
  return.value$iterations <- res.em$iterations
  return.value$info <- res.em$info
  return.value$zeroIntercept <- zeroIntercept
  return.value$normalizeTime <- return.value$normalizeTime
  if(zeroIntercept)
  {
    return.value$intercept <- intercept
  }
  if(!is.null(knots))
  {
    if(normalizeTime)
    {
      return.value$knots <- knots.orig
    }
    else
    {
      return.value$knots <- knots
    }
  }

  class(return.value) <- "sme"

  return(return.value)
}

sme.data.frame <- function(object,tme,ind,verbose=F,lambda.mu=NULL,lambda.v=NULL,maxIter=500,knots=NULL,zeroIntercept=F,deltaEM=1e-3,deltaNM=1e-3,criteria="AICc",initial.lambda.mu=10000,initial.lambda.v=10000,normalizeTime=FALSE,...)
{
  if("variable" %in% names(object))
  {
    ys <- split(object$y,object$variable)
    tmes <- split(object$tme, object$variable)
    inds <- split(object$ind, object$variable)
    return(sme.list(ys,tmes,inds,verbose=verbose,lambda.mu=lambda.mu,lambda.v=lambda.v,maxIter=maxIter,knots=knots,zeroIntercept=zeroIntercept,deltaEM=deltaEM,deltaNM=deltaNM,criteria=criteria,initial.lambda.mu=initial.lambda.mu,initial.lambda.v=initial.lambda.v,normalizeTime=normalizeTime,...))
  }
  else
  {
    return(sme(object=object$y,tme=object$tme,ind=object$ind,verbose=verbose,lambda.mu=lambda.mu,lambda.v=lambda.v,maxIter=maxIter,knots=knots,zeroIntercept=zeroIntercept,deltaEM=deltaEM,deltaNM=deltaNM,criteria=criteria,initial.lambda.mu=initial.lambda.mu,initial.lambda.v=initial.lambda.v,normalizeTime=normalizeTime,...))
  }
}

sme.list <- function(
  object,
  tme,
  ind,
  verbose=F,
  lambda.mu=NULL,
  lambda.v=NULL,
  maxIter=500,
  knots=NULL,
  zeroIntercept=F,
  deltaEM=1e-3,
  deltaNM=1e-3,
  criteria="AICc",
  initial.lambda.mu=10000,
  initial.lambda.v=10000,
  normalizeTime=FALSE,
  numberOfThreads=-1,
  ...)
{
  ys <- object
  tmes <- tme
  inds <- ind

  for(i in 1:length(ys))
  {
    inds[[i]] <- as.factor(inds[[i]])
    inds[[i]] <- droplevels(inds[[i]])
    ys[[i]] <- ys[[i]][order(inds[[i]])]
    tmes[[i]] <- tmes[[i]][order(inds[[i]])]
    inds[[i]] <- inds[[i]][order(inds[[i]])]
  }

  call <- match.call(call=sys.call(-1))

  intercepts <- lapply(tmes,min)  
  if(normalizeTime)
  {
    tmes.orig <- tmes
    tmes <- lapply(tmes,function(tme) (tme - min(tme)) / max(tme))
  }
  
  if(!is.null(knots))
  {
    if(normalizeTime)
    {
      allKnots <- lapply(tmes.orig,function(tme) c(0,(knots - min(tme)) / max(tme),1))
    }
    else
    {
      max.tmes <- sapply(tmes,max)
      min.tmes <- sapply(tmes,min)
      allKnots <- mapply(max.tmes,min.tmes,FUN=function(max.tme,min.tme) c(min.tme,knots,max.tme),SIMPLIFY=FALSE)
    }

    Bs <- lapply(allKnots,ns,knots=knots,intercept=T)

    Gs <- lapply(allKnots,function(a) roughnessMatrix(incidenceMatrix(a)))
    Gs <- mapply(Bs,Gs,FUN=function(B,G) t(B) %*% G %*% B,SIMPLIFY=FALSE)
  }
  else
  {
    Bs <- rep(NA,length(ys))
    allKnots <- rep(NA,length(ys))
    Gs <- lapply(tmes,function(tme) roughnessMatrix(incidenceMatrix(tme)))
  }

  if(zeroIntercept)
  {
    ys.orig <- ys
    tmes.orig2 <- tmes
    inds.orig <- inds
    
    ys <- mapply(ys,tmes,FUN=function(y,tme) y[tme!=min(tme)],SIMPLIFY=FALSE)
    inds <- mapply(inds,tmes,FUN=function(ind,tme) ind[tme!=min(tme)],SIMPLIFY=FALSE)
    tmes <- lapply(tmes,function(tme) tme[tme!=min(tme)])
    
    if(normalizeTime)
    {
      tmes.orig <- lapply(tmes.orig,function(tme.orig) tme.orig[tme.orig!=min(tme.orig)])
    }
    
    Gs <- lapply(Gs,function(G) G[-1,-1])
  }

  M <- length(ys)
  Ni <- sapply(ys,length)

  yis <- mapply(ys,inds,FUN=function(y,ind) split(y,ind),SIMPLIFY=FALSE)
  ni <- sapply(yis, length)

  tmeis <- mapply(tmes,inds,FUN=function(tme,ind) split(tme,ind),SIMPLIFY=FALSE)
  Nij <- lapply(yis,sapply,length)

  if(is.null(knots))
  {
    Xs <- lapply(tmes, function(tme) incidenceMatrix(tme))
    Xis <- mapply(Xs,inds,FUN=function(X,ind) split.data.frame(X,ind),SIMPLIFY=FALSE)
  }
  else
  {
    Xs <- lapply(tmes,function(tme) ns(tme,knots=knots,intercept=T))
    Xis <- mapply(Xs,inds,FUN=function(X,ind) split.data.frame(X,ind),SIMPLIFY=FALSE)
  }

  ps <- sapply(Xs,ncol)

  if(any(is.null(lambda.mu) | is.null(lambda.v)))
  {
    res.em <- .C("SMEOptimizationMultiple",
                 M=as.integer(M),
                 y=as.double(unlist(ys)),
                 tme=as.double(unlist(tmes)),
                 ind=as.integer(unlist(inds)),
                 X=as.double(unlist(Xs)),
                 Ni=as.integer(Ni),
                 ni=as.integer(ni),
                 Nij=as.integer(unlist(Nij)),
                 p=as.integer(ps),
                 lambdaMu=as.double(rep(initial.lambda.mu,M)),
                 lambdaV=as.double(rep(initial.lambda.v,M)),
                 G=as.double(unlist(Gs)),
                 mu=double(sum(ps)),
                 sigmaSquared=double(M),
                 D=double(sum(sapply(ps,"^",2))),
                 v=double(sum(mapply(ps,ni,FUN="*"))),
                 likelihood=double(M),
                 dfMu=double(M),
                 dfV=double(M),
                 iterations=integer(M),
                 maxIterations=as.integer(maxIter),
                 deltaEM=as.double(deltaEM),
                 deltaNM=as.double(deltaNM),
                 criteria=as.integer(match(criteria,c("AIC","AICc","BICN","BICn"))),
                 verbose=as.integer(verbose),
                 info=integer(M),
                 numberOfThreads=as.integer(numberOfThreads),
                 PACKAGE="sme")
  }
  else
  {
    if(length(lambda.mu)==1)
    {
      lambda.mu <- rep(lambda.mu,M)
    }
    if(length(lambda.v)==1)
    {
      lambda.v <- rep(lambda.v,M)
    }
    res.em <- .C("SMEMultiple",
                 M=as.integer(M),
                 y=as.double(unlist(ys)),
                 tme=as.double(unlist(tmes)),
                 ind=as.integer(unlist(inds)),
                 X=as.double(unlist(Xs)),
                 Ni=as.integer(Ni),
                 ni=as.integer(ni),
                 Nij=as.integer(unlist(Nij)),
                 p=as.integer(ps),
                 lambdaMu=as.double(lambda.mu),
                 lambdaV=as.double(lambda.v),
                 G=as.double(unlist(Gs)),
                 mu=double(sum(ps)),
                 sigmaSquared=double(M),
                 D=double(sum(sapply(ps,"^",2))),
                 v=double(sum(mapply(ps,ni,FUN="*"))),
                 likelihood=double(M),
                 dfMu=double(M),
                 dfV=double(M),
                 iterations=integer(M),
                 maxIterations=as.integer(maxIter),
                 deltaEM=as.double(deltaEM),
                 deltaNM=as.double(deltaNM),
                 verbose=as.integer(verbose),
                 info=integer(M),
                 numberOfThreads=as.integer(numberOfThreads),
                 PACKAGE="sme")

  }
  
  mus <- split(res.em$mu,rep(1:M,ps))
  vs <- split(res.em$v,rep(1:M,mapply(ps,ni,FUN="*")))
  iterations <- res.em$iterations
  dfMus <- res.em$dfMu
  dfVs <- res.em$dfV
  likelihoods <- res.em$likelihood
  Ds <- split(res.em$D,rep(1:M,ps*ps))
  sigmas <- res.em$sigmaSquared
  lambdaMus <- res.em$lambdaMu
  lambdaVs <- res.em$lambdaV
  infos <- res.em$info

  return.value <- mapply(
    mus,
    vs,
    Xs,
    Xis,
    Bs,
    ni,
    ps,
    allKnots,
    ys,
    if(normalizeTime){ tmes.orig }else{ tmes },
    inds,
    iterations,
    dfMus,
    dfVs,
    likelihoods,
    Ds,
    sigmas,
    lambdaMus,
    lambdaVs,
    infos,
    intercepts,
    FUN=function(mu,v,X,Xi,B,n,p,allKnots,y,tme,ind,iterations,dfMu,dfV,likelihood,D,sigma,lambdaMu,lambdaV,info,intercept)
    {
      if(is.null(knots))
      {
        coefficients <- rbind(mu,t(matrix(v,ncol=n)))
        colnames(coefficients) <- sort(unique(tme))
      }
      else
      {
        coefficients <- rbind(as.vector(B %*% mu),do.call(rbind,lapply(split(v,rep(1:n,each=p)),function(v) as.vector(B %*% v))))
        colnames(coefficients) <- c(min(tme),knots,max(tme))
      }
      rownames(coefficients) <- c("mu",paste("v",sort(unique(ind)),sep=""))
      fitted <- unlist(mapply(Xi,split(v,rep(1:n,each=p)),FUN=function(Xi,vi) Xi %*% (mu + vi),SIMPLIFY=FALSE))
      residuals <- y - as.vector(fitted)
      data <- data.frame(y=y,tme=tme,ind=ind)
      em <- list(iterations=iterations,likelihood=likelihood,dfMu=dfMu,dfV=dfV,D=matrix(D,nrow=p),sigmaSquared=sigma,lambdaMu=lambdaMu,lambdaV=lambdaV,info=info)
      return.value <- list(
        coefficients=coefficients,
        fitted=fitted,
        resid=residuals,
        data=data,
        #em=em,
        logLik=likelihood,
        df=c(mu=dfMu,v=dfV),
        call=call,
        nobs=length(y),
        smoothingParameters=c(mu=lambdaMu,v=lambdaV),
        parameters=list(sigmaSquared=sigma,D=matrix(D,nrow=p)),
        iterations=iterations,
        info=info,
        zeroIntercept=zeroIntercept,
        normalizeTime=normalizeTime)
      if(!is.null(knots)) return.value$knots <- knots
      if(zeroIntercept){ return.value$intercept <- intercept }
      class(return.value) <- "sme"
      return(return.value)
    },SIMPLIFY=FALSE)

  class(return.value) <- c("sme.list","list")
  return(return.value)
}

getRoughnessMatrix <- function(object)
{
  if(is.null(object$zeroIntercept))
  {
    object$zeroIntercept <- FALSE
  }
  if(is.null(object$normalizeTime))
  {
    object$normalizeTime <- FALSE
  }

  if(object$normalizeTime)
  {
    if(object$zeroIntercept)
    {
      min.tme <- object$intercept
    }
    else
    {
      min.tme <- min(object$tme)
    }
  
    if(!is.null(object$knots))
    {
      knots <- (knots - min.tme) / max(min.tme)
    }
    object$data$tme <- (object$data$tme - min.tme) / max(object$data$tme)
  }

  if(is.null(object$knots))
  {
    X <- incidenceMatrix(object$data$tme)
    Xi <- split.data.frame(X,object$data$ind)

    if(object$zeroIntercept)
    {
      attr(X,"tau") <- c(if(object$normalizeTime){ 0 }else{ object$intercept },attr(X,"tau"))
      attr(X,"M") <- attr(X,"M")+1
      G <- roughnessMatrix(X)
      G <- G[-1,-1]
    }
    else
    {
      G <- roughnessMatrix(X)
    }
  }
  else
  {
    X <- ns(object$data$tme,knots=object$knots,intercept=T)
    Xi <- split.data.frame(X,object$data$ind)
  
    max.tme <- max(object$data$tme)
    min.tme <- min(object$data$tme)
    B <- ns(c(min.tme,object$knots,max.tme),knots=object$knots,intercept=T)
    G <- roughnessMatrix(incidenceMatrix(c(min.tme,object$knots,max.tme)))
    G <- t(B) %*% G %*% B
  }

  G
}

vcov.sme <- function(object,...)
{
  if(is.null(object$knots))
  {
    X <- incidenceMatrix(object$data$tme)
    Xi <- split.data.frame(X,object$data$ind)
  }
  else
  {
    X <- ns(object$data$tme,knots=object$knots,intercept=T)
    Xi <- split.data.frame(X,object$data$ind)
  }

  G <- getRoughnessMatrix(object)
  Dv <- solve(solve(object$parameters$D) + object$smoothingParameters["v"] * G)
  Vi <- lapply(Xi,function(Xi) Xi %*% Dv %*% t(Xi) + diag(object$parameters$sigmaSquared,nrow=nrow(Xi)))
  inverseVi <- lapply(Vi,solve)
  inverseV <- blockDiagMat(inverseVi)
  tXinverseVX <- t(X) %*% inverseV %*% X
  conditionalCovarianceY <- solve(tXinverseVX + object$smoothingParameters["mu"] * G)
  
  vcov <- conditionalCovarianceY %*% tXinverseVX %*% conditionalCovarianceY
  rownames(vcov) <- colnames(coef(object))
  colnames(vcov) <- rownames(vcov)
  return(vcov)
}

rstandard.sme <- function(model,...)
{
  resid(model) / sqrt(model$parameters$sigmaSquared)
}

logLik.sme <- function(object,...)
{
  logLikelihood <- object$logLik
  attr(logLikelihood,"df") <- unname(object$df[1] + object$df[2])
  attr(logLikelihood,"nobs") <- length(resid(object))
  logLikelihood
}

AICc <- function(object)
{
  logLikelihood <- logLik(object)
  df <- attr(logLikelihood,"df")
  nobs <- nobs(object)
  AIC(object) + 2 * df * (df + 1) / (nobs - df - 1)
}

BICn <- function(object,...)
{
  UseMethod("BICn")
}

BICn.sme <- function(object,...)
{
  n <- nrow(coef(object))-1
  logLikelihood <- logLik(object,...)
  -2 * logLikelihood + log(n) * attr(logLikelihood,"df")
}

plot.sme <- function(x,type="model",...)
{
  if(type=="model")
  {
    plotSmeModel(x,...)
  }
  else if(type=="raw")
  {
    plotSmeRaw(x,...)
  }
  else if(type=="diagnostic")
  {
    plotSmeDiagnostic(x,...)
  }
  else
  {
    stop(paste("invalid plot type '",type,"'",sep=""))
  }
}

plotSmeModel <- function(x,xlab="Time",ylab="Y",showIndividuals=T,showConfidenceBands=F,col.meanCurve="red",...)
{
  if(is.null(x$zeroIntercept))
  {
    x$zeroIntercept <- FALSE
  }
  if(x$zeroIntercept)
  {
    x.orig <- x
  
    x$coefficients <- cbind(Intercept=x$intercept,x$coefficients)
    colnames(x$coefficients)[1] <- as.character(x$intercept)
    
    x$data <- rbind(x$data,data.frame(y=0,tme=x$intercept,ind="Intercept"))
  }

  mu <- spline(x=as.numeric(colnames(x$coefficients)),y=x$coefficients[1,],n=500,method="natural")
  fs <- lapply(2:nrow(x$coefficients),function(i){ spline(x=as.numeric(colnames(x$coefficients)),y=x$coefficients[1,] + x$coefficients[i,],method="natural",n=500) })
  ylim <- range(x$data$y, mu$y, sapply(fs,"[[","y"))
  xlim <- range(as.numeric(colnames(x$coefficients)))

  if(showConfidenceBands)
  {
    if(x$zeroIntercept)
    {
      mu.variance <- c(0,diag(vcov(x.orig)))
    }
    else
    {
      mu.variance <- diag(vcov(x))
    }

    upper.band <- spline(x=as.numeric(colnames(x$coefficients)),y=x$coefficients[1,] + 1.96 * sqrt(mu.variance),method="natural",n=500)
    lower.band <- spline(x=as.numeric(colnames(x$coefficients)),y=x$coefficients[1,] - 1.96 * sqrt(mu.variance),method="natural",n=500)
    ylim <- range(ylim, upper.band$y, lower.band$y)
  }

  plot(x=x$data$tme,y=x$data$y,ylim=ylim,xlim=xlim,xlab=xlab,ylab=ylab,...)
  if(showIndividuals)
  {
    for(i in 1:length(fs)){ lines(fs[[i]],lty="dashed") }
  }
  lines(mu,lwd=2,col=col.meanCurve)

  if(showConfidenceBands)
  {
    col.meanCurve.rgb <- col2rgb(col.meanCurve)
    polygon(x=c(upper.band$x,rev(lower.band$x)),y=c(upper.band$y,rev(lower.band$y)),col=rgb(col.meanCurve.rgb[1],col.meanCurve.rgb[2],col.meanCurve.rgb[3],alpha=125,maxColorValue=255),border=NA)
  }
}

plotSmeRaw <- function(x,xlab="Time",ylab="Y",mainTitle="",showModelFits=TRUE,showRawLines=FALSE,...)
{
  n <- nrow(coef(x))-1

  if(is.null(x$zeroIntercept))
  {
    x$zeroIntercept <- FALSE
  }
  if(x$zeroIntercept)
  {
    x$coefficients <- cbind(Intercept=x$intercept,x$coefficients)
    colnames(x$coefficients)[1] <- as.character(x$intercept)
    
    x$data <- rbind(x$data,data.frame(y=0,tme=x$intercept,ind=unique(x$data$ind)))
  }

  xyplot(y ~ tme | ind,
         data=x$data,
         xlab=xlab,
         ylab=ylab,
         main=mainTitle,
         subscripts=T,
         panel=
          function(...)
          {
            panel.xyplot(...)
            if(showRawLines)
            {
              panel.lines(...)
            }
            if(showModelFits)
            {
              ind <- x$data$ind[list(...)$subscripts[1]]
              f <- spline(x=as.numeric(colnames(coef(x))),y=coef(x)[1,] + coef(x)[paste("v",ind,sep=""),],method="natural",n=100)
              panel.lines(x=f$x,y=f$y,lty="dashed")
            }
          },
         ...)
}

plotSmeDiagnostic <- function(x)
{
  layout(matrix(c(1,3,2,4),nrow=2))

  sres <- rstandard(x)

  plot(x=fitted(x),
       y=sres,
       xlab="Fitted values",
       ylab="Standardized residuals",
       main="Standardized residuals\nagainst fitted values")

  plot(x=x$data$tme,
       y=sres,
       xlab="Time",
       ylab="Standardized residuals",
       main="Standardized residuals\nagainst time")

  plot(x=x$data$y,
       y=sres,
       xlab="Response",
       ylab="Standardized residuals",
       main="Standardized residuals\nagainst response")

  qqnorm(sres,main="Q-Q plot for\nstandardized residuals")
  qqline(sres)
  layout(1)
}

blockDiagMat <- function(x)
{
     if(!is.list(x)) stop("x not a list")

     x <- x[as.logical(sapply(x, length))]

     n <- length(x)
     if(n==0) return(matrix(nrow=0,ncol=0))

     d <- array(unlist(lapply(x, dim)), c(2, n))
     rr <- d[1,]
     cc <- d[2,]
     rsum <- sum(rr)
     csum <- sum(cc)
     out <- array(0, c(rsum, csum))
     ind <- array(0, c(4, n))
     rcum <- cumsum(rr)
     ccum <- cumsum(cc)
     ind[1,-1] <- rcum[-n]
     ind[2,] <- rcum
     ind[3,-1] <- ccum[-n]
     ind[4,] <- ccum
     imat <- array(1:(rsum * csum), c(rsum, csum))
     iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2],
(y[3]+1):y[4]], imat=imat)
     iuse <- as.vector(unlist(iuse))
     out[iuse] <- unlist(x)
     return(out)
}

incidenceMatrix <- function(x)
{
  tau <- unique(x)
  tau <- sort(tau)
  M <- length(tau)

  X <- t(sapply(x,function(t){ as.numeric(t == tau) }))
 
  attr(X,"tau") <- tau
  attr(X,"M") <- M

  X
}

roughnessMatrix <- function(x)
{
  tau <- attr(x,"tau")
  K <- attr(x,"M")

  h <- vector(mode="numeric",length=K-1)

  for(r in 1:(K-1))
  {
    h[r] <- tau[r+1] - tau[r] 
  }

  A <- matrix(0, nrow=K, ncol=K-2)

  for(r in 1:(K-2))
  {
    A[r,r] <- 1 / h[r]

    A[r+1,r] <- -( (1 / h[r]) + (1 / h[r+1]) )

    A[r+2,r] <- 1 / h[r+1]
  }

  B <- matrix(0, nrow=K-2, ncol=K-2)

  B[1,1] <- (h[1] + h[2]) / 3

  B[2,1] <- h[2] / 6

  if(K > 4)
  {
    for(r in 1:(K-4))
    {
      B[r,r+1] <- h[r+1] / 6

      B[r+1,r+1] <- (h[r+1] + h[r+2]) / 3

      B[r+2,r+1] <- h[r+2] / 6
    }
  }

  B[K-3,K-2] <- h[K-2] / 6

  B[K-2,K-2] <- (h[K-2] + h[K-1]) / 3

  G <- A %*% solve(B) %*% t(A)
  attr(G,"A") <- A
  attr(G,"B") <- B
  G
}