R/gamm.r

In mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation

##  R routines for mgcv::gamm (c) Simon Wood 2002-2019

### the following two functions are for use in place of log and exp
### in positivity ensuring re-parameterization.... they have `better' 
### over/underflow characteristics, but are still continuous to second
### derivative. 

notExp <- function(x)
# overflow avoiding C2 function for ensuring positivity
{ f <- x
  ind <- x > 1
  f[ind] <- exp(1)*(x[ind]^2+1)/2
  ind <- (x <= 1)&(x > -1)
  f[ind] <- exp(x[ind])
  ind <- (x <= -1)
  x[ind] <- -x[ind] ;f[ind] <-  exp(1)*(x[ind]^2+1)/2; f[ind]<-1/f[ind]
  f
}


notLog <- function(x)
# inverse function of notExp
{ f <- x
  ind <- x> exp(1)
  f[ind] <- sqrt(2*x[ind]/exp(1)-1)
  ind <- !ind & x > exp(-1)
  f[ind] <- log(x[ind])
  ind <- x <= exp(-1)
  x[ind]<- 1/x[ind]; f[ind] <- sqrt(2*x[ind]/exp(1)-1);f[ind] <- -f[ind]
 f
}

## notLog/notExp replacements. 
## around 27/7/05 nlme was modified to use a new optimizer, which fails with 
## indefinite Hessians. This is a problem if smoothing parameters are zero 
## or infinite. The following attempts to make the notLog parameterization 
## non-monotonic, to artificially reduce the likelihood at very large and very
## small parameter values.

## note gamm, pdTens, pdIdnot, notExp and notExp2 .Rd files all modified by
## this change.


notExp2 <- function (x,d=.Options$mgcv.vc.logrange,b=1/d)
## to avoid needing to modify solve.pdIdnot, this transformation must
## maintain the property that 1/notExp2(x) = notExp2(-x)
{ exp(d*sin(x*b))
}

notLog2 <- function(x,d=.Options$mgcv.vc.logrange,b=1/d)
{ x <- log(x)/d
  x <- pmin(1,x)
  x <- pmax(-1,x)
  asin(x)/b
}


#### pdMat class definitions, to enable tensor product smooths to be employed with gamm()
#### Based on various Pinheiro and Bates pdMat classes.

pdTens <- function(value = numeric(0), form = NULL, nam = NULL, data = sys.frame(sys.parent()))
## Constructor for the pdTens pdMat class: 
# the inverse of the scaled random effects covariance matrix for this class
# is given by a weighted sum of the matrices in the list that is the "S" attribute of 
# a pdTens formula. The weights are the exponentials of the class parameters.
# i.e. the inverse of the r.e. covariance matrix is 
#   \sum_i \exp(\theta_i) S_i / \sigma^2 
# The class name relates to the fact that these objects are used with tensor product smooths.  
{
  object <- numeric(0)
  class(object) <- c("pdTens", "pdMat")
  nlme::pdConstruct(object, value, form, nam, data)
}

## Methods for local generics


pdConstruct.pdTens <-
  function(object, value = numeric(0), form = formula(object),
	   nam = nlme::Names(object), data = sys.frame(sys.parent()), ...)
## used to initialize pdTens objects. Note that the initialization matrices supplied
## are (factors of) trial random effects covariance matrices or their inverses.
## Which one is being passed seems to have to be derived from looking at its
##  structure.
## Class tested rather thoroughly with nlme 3.1-52 on R 2.0.0
{
  val <- NextMethod()
  if (length(val) == 0) {               # uninitiliazed object
    class(val) <- c("pdTens","pdMat")
    return(val)
  }
  if (is.matrix(val)) {			# initialize from a positive definite
    S <- attr(form,"S")
    m <- length(S)
    ## codetools gets it wrong about `y'
    y <- as.numeric((crossprod(val)))   # it's a factor that gets returned in val
    lform <- "y ~ as.numeric(S[[1]])"
    if (m>1) for (i in 2:m) lform <- paste(lform," + as.numeric(S[[",i,"]])",sep="")
    lform <- formula(paste(lform,"-1"))
    mod1 <- lm(lform)
    mod1.r2 <- 1-sum(residuals(mod1)^2)/sum((y-mean(y))^2)
    y <- as.numeric(solve(crossprod(val))) ## ignore codetools complaint about this
    mod2 <- lm(lform)
    mod2.r2 <- 1-sum(residuals(mod2)^2)/sum((y-mean(y))^2)
    ## `value' and `val' can relate to the cov matrix or its inverse:
    ## the following seems to be only way to tell which.
    #if (summary(mod2)$r.sq>summary(mod1)$r.sq) mod1<-mod2
    if (mod2.r2 > mod1.r2) mod1 <- mod2
    value <- coef(mod1)  
    value[value <=0] <- .Machine$double.eps * mean(as.numeric(lapply(S,function(x) max(abs(x)))))
    value <- notLog2(value)
    attributes(value) <- attributes(val)[names(attributes(val)) != "dim"]
    class(value) <- c("pdTens", "pdMat")
    return(value)
  }
  m <- length(attr(form,"S"))
  if ((aux <- length(val)) > 0) {
    if (aux && (aux != m)) {
      stop(gettextf("An object of length %d does not match the required parameter size",aux))
    }
  }
  class(val) <- c("pdTens","pdMat")
  val
}


pdFactor.pdTens <- function(object)
## The factor of the inverse of the scaled r.e. covariance matrix is returned here
## it should be returned as a vector. 
{ sp <- as.vector(object)
  m <- length(sp)
  S <- attr(formula(object),"S")
  value <- S[[1]]*notExp2(sp[1])
  if (m>1) for (i in 2:m) value <- value + notExp2(sp[i])*S[[i]] 
  if (sum(is.na(value))>0) warning("NA's in pdTens factor")
  value <- (value+t(value))/2
  c(t(mroot(value,rank=nrow(value))))
}


pdMatrix.pdTens <-
  function(object, factor = FALSE) 
# the inverse of the scaled random effect covariance matrix is returned here, or
# its factor if factor==TRUE. If A is the matrix being factored and B its
# factor, it is required that A=B'B (not the mroot() default!)

{
  if (!nlme::isInitialized(object)) {
    stop("Cannot extract the matrix from an uninitialized object")
  }
  sp <- as.vector(object)
  m <- length(sp)
  S <- attr(formula(object),"S")
  value <- S[[1]]*notExp2(sp[1])   
  if (m>1) for (i in 2:m) value <- value + notExp2(sp[i])*S[[i]]  
 
  value <- (value + t(value))/2 # ensure symmetry
  if (sum(is.na(value))>0) warning("NA's in pdTens matrix")
  if (factor) {
    value <- t(mroot(value,rank=nrow(value)))
  } 
  dimnames(value) <- attr(object, "Dimnames")
  value
}

#### Methods for standard generics

coef.pdTens <-
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) NextMethod()
  else {
    val <- notExp2(as.vector(object))
    names(val) <- paste("sp.",1:length(val), sep ="")
    val
  }
}

summary.pdTens <-
  function(object, structName = "Tensor product smooth term", ...)
{
  NextMethod(object, structName, noCorrelation=TRUE)
}


# .... end of pdMat definitions for tensor product smooths


### pdIdnot: multiple of the identity matrix - the parameter is
### the notLog2 of the multiple. This is directly modified form 
### Pinheiro and Bates pdIdent class.

####* Constructor

pdIdnot <-
  ## Constructor for the pdIdnot class
  function(value = numeric(0), form = NULL, nam = NULL, data = sys.frame(sys.parent()))
{ #cat(" pdIdnot  ")
  object <- numeric(0)
  class(object) <- c("pdIdnot", "pdMat")
  nlme::pdConstruct(object, value, form, nam, data)
}

####* Methods for local generics

corMatrix.pdIdnot <-
  function(object, ...)
{
  if (!nlme::isInitialized(object)) {
    stop("Cannot extract the matrix from an uninitialized pdMat object")
  }
  if (is.null(Ncol <- attr(object, "ncol"))) {
    stop(paste("Cannot extract the matrix with uninitialized dimensions"))
  }
  val <- diag(Ncol)
## REMOVE sqrt() to revert ...
  attr(val, "stdDev") <- rep(sqrt(notExp2(as.vector(object))), Ncol)
  if (length(nm <- nlme::Names(object)) == 0) {
    len <- length(as.vector(object)) 
    nm <- paste("V", 1:len, sep = "")
    dimnames(val) <- list(nm, nm)
  }
  names(attr(val, "stdDev")) <- nm
  val
}

pdConstruct.pdIdnot <-
  function(object, value = numeric(0), form = formula(object),
	   nam = nlme::Names(object), data = sys.frame(sys.parent()), ...)
{ #cat(" pdConstruct.pdIdnot  ")
  val <- NextMethod()
  if (length(val) == 0) {			# uninitialized object
    if ((ncol <- length(nlme::Names(val))) > 0) {
      attr(val, "ncol") <- ncol
    }
    return(val)
  }
  if (is.matrix(val)) {
#    value <- notLog2(sqrt(mean(diag(crossprod(val)))))
    value <- notLog2(mean(diag(crossprod(val)))) ## REPLACE by above to revert
    attributes(value) <- attributes(val)[names(attributes(val)) != "dim"]
    attr(value, "ncol") <- dim(val)[2]
    class(value) <- c("pdIdnot", "pdMat")
    return(value)
  }
  if (length(val) > 1) {
    stop(paste("An object of length", length(val),
	       "does not match the required parameter size"))
  }
  if (((aux <- length(nlme::Names(val))) == 0) && is.null(formula(val))) {
    stop(paste("Must give names when initializing pdIdnot from parameter.",
	       "without a formula"))
  } else {
    attr(val, "ncol") <- aux
  }
  val
}

pdFactor.pdIdnot <-
  function(object)
{ ## UNCOMMENT first line, comment 2nd to revert
  # notExp2(as.vector(object)) * diag(attr(object, "ncol"))
  #cat(" pdFactor.pdIdnot  ")
  sqrt(notExp2(as.vector(object))) * diag(attr(object, "ncol"))
}

pdMatrix.pdIdnot <-
  function(object, factor = FALSE)
{ #cat("  pdMatrix.pdIdnot  ")
  if (!nlme::isInitialized(object)) {
    stop("Cannot extract the matrix from an uninitialized pdMat object")
  }
  if (is.null(Ncol <- attr(object, "ncol"))) {
    stop(paste("Cannot extract the matrix with uninitialized dimensions"))
  }
  value <- diag(Ncol)
    
  ## REPLACE by #1,#2,#3 to revert
  if (factor) {
   #1  value <- notExp2(as.vector(object)) * value
   #2  attr(value, "logDet") <- Ncol * log(notExp2(as.vector(object)))
   value <- sqrt(notExp2(as.vector(object))) * value
   attr(value, "logDet") <- Ncol * log(notExp2(as.vector(object)))/2
  } else {
   #3 value <- notExp2(as.vector(object))^2 * value
    value <- notExp2(as.vector(object)) * value
  }
  dimnames(value) <- attr(object, "Dimnames")
  value
}

####* Methods for standard generics

coef.pdIdnot <-
  function(object, unconstrained = TRUE, ...)
{ #cat(" coef.pdIdnot    ")
  if (unconstrained) NextMethod()
  else structure(notExp2(as.vector(object)),
           names = c(paste("sd(", deparse(formula(object)[[2]],backtick=TRUE),")",sep = "")))
}

Dim.pdIdnot <-
  function(object, ...)
{
  if (!is.null(val <- attr(object, "ncol"))) {
    c(val, val)
  } else {
    stop("Cannot extract the dimensions")
  }
}

logDet.pdIdnot <-
  function(object, ...)
{ ## REMOVE /2 to revert ....
  attr(object, "ncol") * log(notExp2(as.vector(object)))/2
  
}

solve.pdIdnot <-
  function(a, b, ...)
{
  if (!nlme::isInitialized(a)) {
    stop("Cannot extract the inverse from an uninitialized object")
  }
  atr <- attributes(a)
  a <- -coef(a, TRUE)
  attributes(a) <- atr
  a
}

summary.pdIdnot <-
  function(object, structName = "Multiple of an Identity", ...)
{ #cat("  summary.pdIdnot  ")
  # summary.pdMat(object, structName, noCorrelation = TRUE)

  ## ... summary.pdMat is not exported in the nlme NAMESPACE file, so....

  NextMethod(object, structName, noCorrelation=TRUE)
}

### end of pdIdnot class

smooth2random <- function(object,vnames,type=1) UseMethod("smooth2random")

smooth2random.fs.interaction <- function(object,vnames,type=1) {
## conversion method for smooth-factor random interactions.
## For use with gamm4, this needs to generate a sparse version of
## each full model matrix, with smooth coefs re-ordered so that the 
## penalties are not interwoven, but blocked (i.e. this ordering is 
## as for gamm case). 
  if (object$fixed) return(list(fixed=TRUE,Xf=object$X))
  ## If smooth constructor was not called with "gamm" attribute set,
  ## then we need to reset model matrix to base model matrix.
  if (!is.null(object$Xb)) {
    object$X <- object$Xb
    object$S <- object$base$S
    if (!is.null(object$S.scale)&&length(object$S)>0) for (i in 1:length(object$S)) object$S[[i]] <- object$S[[i]]/object$S.scale[i] 
  }  
  colx <- ncol(object$X) 
  diagU <- rep(1,colx)
  ind <- 1:colx 
  ## flev <- levels(object$fac)
  n.lev <- length(object$flev)
  if (type==2) {
    ## index which params in fit parameterization are penalized by each penalty.
    ## e.g. pen.ind==1 is TRUE for each param penalized by first penalty and
    ## FALSE otherwise...
    pen.ind <- rep(ind*0,n.lev)
  } else pen.ind <- NULL
  random <- list()
  k <- 1
  rinc <- rind <- rep(0,0)
  for (i in 1:length(object$S)) { ## work through penalties
    indi <- ind[diag(object$S[[i]])!=0] ## index of penalized cols
   
    X <- object$X[,indi,drop=FALSE] ## model matrix for this component
    D <- diag(object$S[[i]])[indi]
    diagU[indi] <- 1/sqrt(D) ## transform that reduces penalty to identity
    X <- X%*%diag(diagU[indi],ncol=length(indi))
    term.name <- new.name("Xr",vnames)
    vnames <- c(vnames,term.name)
    rind <- c(rind,k:(k+ncol(X)-1))
    rinc <- c(rinc,rep(ncol(X),ncol(X)))
    k <- k + n.lev * ncol(X)
    if (type==1) { ## gamm form for use with lme 
      ## env set to avoid 'save' saving whole environment to file...
      form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
      fname <- new.name(object$fterm,vnames)
      vnames <- c(vnames,fname)
      random[[i]] <- pdIdnot(form)
      names(random)[i] <- fname        ## supplied factor name
      attr(random[[i]],"group") <- object$fac  ## factor supplied as part of term 
      attr(random[[i]],"Xr.name") <- term.name
      attr(random[[i]],"Xr") <- X
    } else { ## gamm4 form --- whole sparse matrices
      ## Xr <- as(matrix(0,nrow(X),0),"dgCMatrix") - deprecated, use...
      Xr <- as(as(as(matrix(0,nrow(X),0), "dMatrix"), "generalMatrix"), "CsparseMatrix")
      ii <- 0
      for (j in 1:n.lev) { ## assemble full sparse model matrix
        ## Xr <- cbind2(Xr,as(X*as.numeric(object$fac==object$flev[j]),"dgCMatrix")) - deprecated
	Xr <- cbind2(Xr,
	  as(as(as(X*as.numeric(object$fac==object$flev[j]), "dMatrix"), "generalMatrix"), "CsparseMatrix"))
        pen.ind[indi+ii] <- i;ii <- ii + colx
      }
      random[[i]] <- if (is.null(object$Xb)) Xr else as(Xr,"matrix") 
      names(random)[i] <- term.name
      attr(random[[i]],"s.label") <- object$label
    }
  }
  if (type==2) {
    ## expand the rind (rinc not needed)
    ind <- 1:length(rind)
    ni <- length(ind)
    rind <- rep(rind,n.lev)
    if (n.lev>1) for (k in 2:n.lev) { 
      rind[ind+ni] <- rind[ind]+rinc
      ind <- ind + ni
    }
    D <- rep(diagU,n.lev)
  } else D <- diagU ## b_original = D*b_fit

  Xf <- matrix(0,nrow(object$X),0)
  list(rand=random,trans.D=D,Xf=Xf,fixed=FALSE,rind=rind,rinc=rinc,
       pen.ind=pen.ind) ## pen.ind==i is TRUE for coefs penalized by ith penalty
} ## smooth2random.fs.interaction

smooth2random.t2.smooth <- function(object,vnames,type=1) {
## takes a smooth object and turns it into a form suitable for estimation as a random effect
## vnames is a list of names to avoid when assigning variable names.
## type==1 indicates an lme random effect.
## Returns 1. a list of random effects, including grouping factors, and 
##            a fixed effects matrix. Grouping factors, model matrix and 
##            model matrix name attached as attributes, to each element. 
##         2. rind: and index vector such that if br is the vector of 
##            random coefficients for the term, br[rind] is the coefs in 
##            order for this term. rinc - dummy here.
##         3. A matrix, trans.D, that transforms coefs, in order [rand1, rand2,... fix]
##            back to original parameterization. If null, then not needed.
##         4. A matrix Xf for the fixed effects, if any.
##         5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is
##            not returned.
## This version deals only with t2 smooths conditioned on a whole 
## dataset dummy factor.
## object must contain model matrix for smooth.
  if (object$fixed) return(list(fixed=TRUE,Xf=object$X))

  fixed <- rep(TRUE,ncol(object$X))
  random <- list()
  diagU <- rep(1,ncol(object$X))
  ind <- 1:ncol(object$X)
  pen.ind <- ind*0
  n.para <- 0
  for (i in 1:length(object$S)) { ## work through penalties
    indi <- ind[diag(object$S[[i]])!=0] ## index of penalized cols
    pen.ind[indi] <- i
    X <- object$X[,indi,drop=FALSE] ## model matrix for this component
    D <- diag(object$S[[i]])[indi]
    diagU[indi] <- 1/sqrt(D) ## transform that reduces penalty to identity
    X <- X%*%diag(diagU[indi])
    fixed[indi] <- FALSE
    term.name <- new.name("Xr",vnames)
    group.name <- new.name("g",vnames)
    vnames <- c(vnames,term.name,group.name)
    if (type==1) { ## gamm form for lme
      ## env set to avoid 'save' saving whole environment to file...
      form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
      random[[i]] <- pdIdnot(form)
      names(random)[i] <- group.name
      attr(random[[i]],"group") <- factor(rep(1,nrow(X)))
      attr(random[[i]],"Xr.name") <- term.name
      attr(random[[i]],"Xr") <- X  
    } else { ## lmer form as used by gamm4
      random[[i]] <- X 
      names(random)[i] <- term.name
      attr(random[[i]],"s.label") <- object$label
    }
    n.para <- n.para + ncol(X)
  }
  if (sum(fixed)) { ## then there are fixed effects!
    Xf <- object$X[,fixed,drop=FALSE]
  } else Xf <- matrix(0,nrow(object$X),0)
  list(rand=random,trans.D=diagU,Xf=Xf,fixed=FALSE,
       rind=1:n.para,rinc=rep(n.para,n.para),pen.ind=pen.ind)
} ## smooth2random.t2.smooth

smooth2random.mgcv.smooth <- function(object,vnames,type=1) {
## takes a smooth object and turns it into a form suitable for estimation as a random effect
## vnames is a list of names to avoid when assigning variable names.
## type==1 indicates an lme random effect.
## Returns 1. a list of random effects, including grouping factors, and 
##            a fixed effects matrix. Grouping factors, model matrix and 
##            model matrix name attached as attributes, to each element. 
##         2. rind: an index vector such that if br is the vector of 
##            random coefficients for the term, br[rind] is the coefs in 
##            order for this term. rinc - dummy here.
##         3. A matrix, U, + vec D that transforms coefs, in order [rand1, rand2,... fix]
##            back to original parameterization. b_origonal = U%*%(D*b_fit) 
##         4. A matrix Xf for the fixed effects, if any.
##         5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is
##            not returned.
## This version deals only with single penalty smooths conditioned on a whole 
## dataset dummy factor.
## object must contain model matrix for smooth.
  if (object$fixed) return(list(fixed=TRUE,Xf=object$X))

  if (length(object$S)>1) stop("Can not convert this smooth class to a random effect")

  ## reparameterize so that unpenalized basis is separated out and at end...
  ev <- eigen(object$S[[1]],symmetric=TRUE)
  ## following is a hack for developers calling smooth2random for fit and then again
  ## for prediction, rather than as intended (do this on different machines and...)
  if (ev$vectors[1,1]<0) ev$vectors <- -ev$vectors 
  null.rank <- object$df - object$rank
  p.rank <- object$rank
  if (p.rank>ncol(object$X)) p.rank <- ncol(object$X)
  U <- ev$vectors
  D <- c(ev$values[1:p.rank],rep(1,null.rank))
  D <- 1/sqrt(D)
  UD <- t(t(U)*D)      ## U%*%[b,beta] returns coefs in original parameterization 
  X <- object$X%*%UD 
  if (p.rank<object$df) Xf <- X[,(p.rank+1):object$df,drop=FALSE] else 
                        Xf <- matrix(0,nrow(object$X),0) # no fixed terms left
  
  term.name <- new.name("Xr",vnames)
  if (type==1) { ## gamm form for lme
    ## env set to avoid 'save' saving whole environment with fitted model object
    form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
    random <- list(pdIdnot(form))
    group.name <- new.name("g",vnames)
    names(random) <- group.name
    attr(random[[1]],"group") <- factor(rep(1,nrow(X)))
    attr(random[[1]],"Xr.name") <- term.name
    attr(random[[1]],"Xr") <- X[,1:p.rank,drop=FALSE]
  } else { ## lmer form as used in gamm4
    random <- list(X[,1:p.rank,drop=FALSE]) 
    names(random)[1] <- term.name
    attr(random[[1]],"s.label") <- object$label
  }
  rind <- 1:p.rank
  pen.ind <- rep(0,ncol(object$X))
  pen.ind[rind] <- 1
  rinc <- rep(p.rank,p.rank)
  list(rand=random, ## the random effects for this term
         Xf=Xf, ## the fixed effect model matrix for this term
         trans.U=U,   ## orthog matrix mapping coefs back to original parameterization 
         trans.D=D,   ## back trans vector (b_original = U%*%(D*b_fit))
         fixed=FALSE,rind=rind,rinc=rinc,
         pen.ind=pen.ind)
} ## end smooth2random.mgcv.smooth


smooth2random.tensor.smooth <- function(object,vnames,type=1) {
## takes a smooth object and turns it into a form suitable for estimation as a random effect
## vnames is a list of names to avoid when assigning variable names.
## type==1 indicates an lme random effect.
## Returns 1. a list of random effects, including grouping factors, and 
##            a fixed effects matrix. Grouping factors, model matrix and 
##            model matrix name attached as attributes, to each element. 
##         2. rind: an index vector such that if br is the vector of 
##            random coefficients for the term, br[rind] is the coefs in 
##            order for this term. rinc - dummy here.
##         3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix]
##            back to original parameterization. If null, then not needed.
##         4. A matrix Xf for the fixed effects, if any.
##         5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is
##            not returned.
## This version deals only with te smooths conditioned on a whole 
## dataset dummy factor.
## object must contain model matrix for smooth.

  if (type==2) stop("te smooths not useable with gamm4: use t2 instead")

  if (sum(object$fx)==length(object$fx)) return(list(fixed=TRUE,Xf=object$X)) else
  if (sum(object$fx)!=0) warning("gamm can not fix only some margins of tensor product.")
     
  ## first sort out the re-parameterization...
  sum.S <- object$S[[1]]/mean(abs(object$S[[1]]))
  #null.rank <- ncol(object$margin[[1]]$X)-object$margin[[1]]$rank ## null space rank
  # bs.dim <- ncol(object$margin[[1]]$X)
  if (length(object$S)>1) for (l in 2:length(object$S)) { 
    sum.S <- sum.S + object$S[[l]]/mean(abs(object$S[[l]]))
    #dfl <- ncol(object$margin[[l]]$X) ## actual df of term (`df' may not be set by constructor)
    #null.rank <- null.rank * (dfl-object$margin[[l]]$rank) 
    #bs.dim <- bs.dim * dfl
  }
  null.rank <- object$null.space.dim
  #null.rank <- null.rank - bs.dim + object$df
  ##sum.S <- (sum.S+t(sum.S))/2 # ensure symmetry
  ev <- eigen(sum.S,symmetric=TRUE)
  ## following is a hack for developers calling smooth2random for fit and then again
  ## for prediction, rather than as intended...
  if (ev$vectors[1,1]<0) ev$vectors <- -ev$vectors 
  p.rank <- ncol(object$X) - null.rank
  if (p.rank>ncol(object$X)) p.rank <- ncol(object$X)
  U <- ev$vectors
  D <- c(ev$values[1:p.rank],rep(1,null.rank))
  if (sum(D<=0)) stop(
  "Tensor product penalty rank appears to be too low: please email Simon.Wood@R-project.org with details.")
  ## D <- 1/sqrt(D)
  U <- U ## maps coefs back to untransformed versions
  X <- object$X%*%U
  if (p.rank<ncol(X)) Xf <- X[,(p.rank+1):ncol(X),drop=FALSE] else 
                      Xf <- matrix(0,nrow(X),0) # no fixed terms left
 
  for (l in 1:length(object$S)) {   # transform penalty explicitly
    object$S[[l]] <- (t(U)%*%object$S[[l]]%*%U)[1:p.rank,1:p.rank]
    object$S[[l]] <- (object$S[[l]]+t(object$S[[l]]))/2
  }
  
  term.name <- new.name("Xr",vnames)
  form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
  attr(form,"S") <- object$S   ## this class needs penalty matrices to be supplied
  random <- list(pdTens(form)) ## lme random effect class

  group.name <- new.name("g",vnames)
  names(random) <- group.name ## grouping factor name
  attr(random[[1]],"group") <- factor(rep(1,nrow(X))) ## grouping factor
  attr(random[[1]],"Xr.name") <- term.name           ## random effect matrix name 
  attr(random[[1]],"Xr") <- X[,1:p.rank,drop=FALSE] ## random effect model matrix
  rind <- 1:p.rank
  rinc <- rep(p.rank,p.rank)
  list(rand=random, ## the random effects for this term
         Xf=Xf, ## the fixed effect model matrix for this term
         trans.U=U,   ## the matrix mapping coefs back to original parameterization   
         trans.D=rep(1,ncol(U)),
         fixed=FALSE,rind=rind,rinc=rinc)
} ## end smooth2random.tensor.smooth



gamm.setup <- function(formula,pterms,
                      data=stop("No data supplied to gamm.setup"),knots=NULL,
                     parametric.only=FALSE,absorb.cons=FALSE)
## set up the model matrix, penalty matrices and auxilliary information about the smoothing bases
## needed for a gamm fit.
## NOTE: lme can't deal with offset terms.
## There is an implicit assumption that any rank deficient penalty does not penalize 
## the constant term in a basis. 
## 1. Calls gam.setup, as for a gam to produce object G suitable for estimating a gam.
## 2. Works through smooth list, G$smooth, modifying so that... 
##    i) Smooths are reparameterized to have identity, or tensor, penalty matrices,
##       and separated fixed and random components with transform information 
##       stored in the smooth.
##   ii) The smooth stores the index range of its coefficients in overall
##       fixed and r.e. coef vectors.
##   iii) pdIdnot or pdTens terms are stored in a separate `G$random' list, which
##        refers, by name, to r.e. model matrices stored directly in G.
##   iv) an overall fixed effect matrix X is accumulated.   
{ 
  ## first simply call `gam.setup'....

  G <- gam.setup(formula,pterms,
                 data=data,knots=knots,sp=NULL,
                 min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,gamm.call=TRUE)
 
  if (!is.null(G$L)) stop("gamm can not handle linked smoothing parameters (probably from use of `id' or adaptive smooths)")
  # now perform re-parameterization...

  first.f.para <- G$nsdf+1 
  first.r.para <- 1
 
  G$Xf <- G$X # full GAM model matrix, treating smooths as fixed effects
  random <- list()
  
  if (G$nsdf>0) ind <- 1:G$nsdf else ind <- rep(0,0)  
  X <- G$X[,ind,drop=FALSE] # accumulate fixed effects into here

  xlab <- rep("",0)
  
  ## first have to create a processing order, so that any smooths conditional on 
  ## multi-level factors are processed last, and hence end up at the end of the 
  ## random list (right is nested in left in this list!)
  if (G$m>0) {
    pord <- 1:G$m
    done <- rep(FALSE,length(pord))
    k <- 0
    f.name <- NULL
    for (i in 1:G$m) if (is.null(G$smooth[[i]]$fac)) { 
         k <- k + 1
         pord[k] <- i 
         done[i] <- TRUE
    } else {
      if (is.null(f.name)) f.name <- G$smooth[[i]]$fterm
      else if (f.name!=G$smooth[[i]]$fterm) stop("only one level of smooth nesting is supported by gamm")
      if (!is.null(attr(G$smooth[[i]],"del.index"))) stop("side conditions not allowed for nested smooths")
    }   
    if (k < G$m) pord[(k+1):G$m] <- (1:G$m)[!done] 
    ## .... ordered so that nested smooths are last
  }

  if (G$m) for (i in 1:G$m) { ## work through the smooths
    
    sm <- G$smooth[[pord[i]]]
    sm$X <- G$X[,sm$first.para:sm$last.para,drop=FALSE]
    rasm <- smooth2random(sm,names(data)) ## convert smooth to random effect and fixed effects
    sm$fixed <- rasm$fixed
  
    if (!is.null(sm$fac)) { 
      flev <- levels(sm$fac) ## grouping factor for smooth
      ##n.lev <- length(flev)
    } ##else n.lev <- 1
 
    ## now append constructed variables to model frame and random effects to main list
    n.para <- 0 ## count random coefficients
   ## rinc <- rind <- rep(0,0)
    if (!sm$fixed) {
     # kk <- 1;
      for (k in 1:length(rasm$rand)) {
         group.name <- names(rasm$rand)[k]
         group <- attr(rasm$rand[[k]],"group")
         Xr.name <- attr(rasm$rand[[k]],"Xr.name")
         Xr <- attr(rasm$rand[[k]],"Xr")
         attr(rasm$rand[[k]],"group") <- attr(rasm$rand[[k]],"Xr") <- 
         attr(rasm$rand[[k]],"Xr.name") <- NULL
     #    rind <- c(rind,kk:(kk+ncol(Xr)-1))
     #    rinc <- c(rinc,rep(ncol(Xr),ncol(Xr))) ## increment for rind
     #    kk <- kk + n.lev * ncol(Xr)
         n.para <- n.para + ncol(Xr)
         data[[group.name]] <- group
         data[[Xr.name]] <- Xr
      }
      random <- c(random,rasm$rand)
      sm$trans.U <- rasm$trans.U ## matrix mapping fit coefs back to original
      sm$trans.D <- rasm$trans.D ## so b_original = U%*%(D*b_fit)
    }

    if (ncol(rasm$Xf)) { ## lme requires names
      Xfnames <- rep("",ncol(rasm$Xf)) 
      k <- length(xlab)+1
      for (j in 1:ncol(rasm$Xf)) {
        xlab[k] <- Xfnames[j] <-
        new.name(paste(sm$label,"Fx",j,sep=""),xlab)
        k <- k + 1
      } 
      colnames(rasm$Xf) <- Xfnames
    }

    X <- cbind(X,rasm$Xf) # add fixed model matrix to overall fixed X

    ## update the counters indicating which elements of the whole model 
    ## fixed and random coef vectors contain the coefs for this smooth.
    ## note convention that smooth coefs are [random, fixed]    

    sm$first.f.para <- first.f.para
    first.f.para <- first.f.para + ncol(rasm$Xf)
    sm$last.f.para <- first.f.para - 1 ## note less than sm$first.f.para => no fixed

    sm$rind <- rasm$rind - 1 + first.r.para
    sm$rinc <- rasm$rinc 
 #   sm$first.r.para <- first.r.para
    first.r.para <- first.r.para+n.para
 #   sm$last.r.para <- first.r.para-1

    sm$n.para <- n.para

    ## convention is that random coefs for grouped smooths will be 
    ## packed [coefs for level 1, coefs for level 2, ...]
    ## n.para is number of random paras at each level.
    ## so coefs for ith level will be indexed by 
    ## rind + (i-1)*n.para
    ## first.r.para:last.r.para + (i-1)*n.para

    if (!is.null(sm$fac)) { ## there is a grouping factor for this smooth
      ## have to up this first.r.para to allow a copy of coefs for each level of fac...
      first.r.para <- first.r.para + n.para*(length(flev)-1) 
    }   
 
    sm$X <- NULL ## delete model matrix
  
    if (G$m>0) G$smooth[[pord[i]]] <- sm  ## replace smooth object with extended version 
  }
 
  G$random <- random
  G$X <- X  ## fixed effects model matrix
  G$data <- data
  if (G$m>0) G$pord <- pord ## gamm needs to run through smooths in same order as here

  G
} ## end of gamm.setup



varWeights.dfo <- function(b,data)
## get reciprocal *standard deviations* implied by the estimated variance
## structure of an lme object, b, in *original data frame order*.
{  w <- nlme::varWeights(b$modelStruct$varStruct) 
   # w is not in data.frame order - it's in inner grouping level order
   group.name <- names(b$groups) # b$groups[[i]] doesn't always retain factor ordering
   ind <- NULL
   order.txt <- paste("ind<-order(data[[\"",group.name[1],"\"]]",sep="")
   if (length(b$groups)>1) for (i in 2:length(b$groups)) 
   order.txt <- paste(order.txt,",data[[\"",group.name[i],"\"]]",sep="")
   order.txt <- paste(order.txt,")")
   eval(parse(text=order.txt))
   w[ind] <- w # into data frame order
   w
}

extract.lme.cov2<-function(b,data=NULL,start.level=1)
# function to extract the response data covariance matrix from an lme fitted
# model object b, fitted to the data in data. "inner" == "finest" grouping 
# start.level is the r.e. grouping level at which to start the construction, 
# levels outer to this will not be included in the calculation - this is useful
# for gamm calculations
# 
# This version aims to be efficient, by not forming the complete matrix if it
# is diagonal or block diagonal. To this end the matrix is returned in a form
# that relates to the data re-ordered according to the coarsest applicable 
# grouping factor. ind[i] gives the row in the original data frame
# corresponding to the ith row/column of V. 
# V is either returned as an array, if it's diagonal, a matrix if it is
# a full matrix or a list of matrices if it is block diagonal.
{ if (!inherits(b,"lme")) stop("object does not appear to be of class lme")
  if (is.null(data)) {
    na.act <- na.action(b)
    data <- if (is.null(na.act)) b$data else b$data[-na.act,] 
  }
  grps <- nlme::getGroups(b) # labels of the innermost groupings - in data frame order
  n <- length(grps)    # number of data
  n.levels <- length(b$groups) # number of nested grouping factors (random effects only)
#  if (n.levels<start.level) { ## then examine correlation groups
    n.corlevels <- if (is.null(b$modelStruct$corStruct)) 0 else
                   length(all.vars(nlme::getGroupsFormula(b$modelStruct$corStruct)))
#  } else n.corlevels <- 0 ## used only to signal irrelevance
  ## so at this stage n.corlevels > 0 iff it determines the coarsest grouping
  ## level if > start.level. 
  if (n.levels<n.corlevels) { ## then cor groups are finest
    ## correlation groups are awkward. The groups returned by 
    ## getGroups(b$modelStruct$corStruct) are not in data frame 
    ## order, but in a sorted order. This is odd because if you
    ## take the original uninitialized corStruct and Initialize 
    ## using what is stored in b$data, and then apply getGroups,
    ## the groups are in data frame order (suggesting one route
    ## for getting the groups into data frame order). The approach
    ## taken here simply re-constructs the groups in data frame
    ## order, thereby avoiding the need to store the uninitialized 
    ## corStruct. Note that covariates of corStructs do not 
    ## cause re-ordering: only grouping variables do that. 
    ## The above has been tested at nlme_3.1-109
    ## grps <- nlme::getGroups(b$modelStruct$corStruct) # replace grps (but not n.levels)
    getGroupsFormula(b$modelStruct$corStruct)
    vnames <- all.vars(nlme::getGroupsFormula(b$modelStruct$corStruct))
    ## attr(b$modelStruct$corStruct,"formula") # would include covariates
    lab <- paste(eval(parse(text=vnames[1]),envir=b$data))
    if (length(vnames)>1) for (i in 2:length(vnames)) {
      lab <- paste(lab,"/",eval(parse(text=vnames[i]),envir=b$data),sep="")
    }
    grps <- factor(lab)
  }
  if (n.levels >= start.level||n.corlevels >= start.level) {
    if (n.levels >= start.level)
    Cgrps <- nlme::getGroups(b,level=start.level) # outer grouping labels (dforder) 
    else Cgrps <- grps
    #Cgrps <- nlme::getGroups(b$modelStruct$corStruct) # ditto
    Cind <- sort(as.numeric(Cgrps),index.return=TRUE)$ix
    # Cind[i] is where row i of sorted Cgrps is in original data frame order 
    rCind <- 1:n; rCind[Cind] <- 1:n
    # rCind[i] is location of ith original datum in the coarse ordering
    ## CFgrps <- grps[Cind] # fine group levels in coarse group order (unused!!)
    Clevel <- levels(Cgrps) # levels of coarse grouping factor
    n.cg <- length(Clevel)  # number of outer groups
    size.cg <- array(0,n.cg)  
    for (i in 1:n.cg) size.cg[i] <- sum(Cgrps==Clevel[i]) # size of each coarse group
    ## Cgrps[Cind] is sorted by coarsest grouping factor level
    ## so e.g. y[Cind] would be data in c.g.f. order
  } else { n.cg <- 1; Cind<-1:n }
  if (is.null(b$modelStruct$varStruct)) w<-rep(b$sigma,n) ### 
  else {
    w <- 1/nlme::varWeights(b$modelStruct$varStruct) 
    # w is not in data.frame order - it's in inner grouping level order
    group.name<-names(b$groups) # b$groups[[i]] doesn't always retain factor ordering
    order.txt <- paste("ind<-order(data[[\"",group.name[1],"\"]]",sep="")
    if (length(b$groups)>1) for (i in 2:length(b$groups)) 
    order.txt <- paste(order.txt,",data[[\"",group.name[i],"\"]]",sep="")
    order.txt <- paste(order.txt,")")
    eval(parse(text=order.txt))
    w[ind] <- w # into data frame order
    w <- w*b$sigma
  }
  w <- w[Cind] # re-order in coarse group order
  if (is.null(b$modelStruct$corStruct)) V <- array(1,n) 
  else {
    c.m <- nlme::corMatrix(b$modelStruct$corStruct) # correlation matrices for each innermost group
    if (!is.list(c.m)) { # copy and re-order into coarse group order
      V <- c.m;V[Cind,] -> V;V[,Cind] -> V 
    } else { 
      V<-list()   # V[[i]] is cor matrix for ith coarse group
      ind <- list() # ind[[i]] is order index for V[[i]] 
      for (i in 1:n.cg) { 
        V[[i]] <- matrix(0,size.cg[i],size.cg[i]) 
        ind[[i]] <- 1:size.cg[i]
      }
      # Voff[i] is where, in coarse order data, first element of V[[i]]
      # relates to ... 
      Voff <- cumsum(c(1,size.cg)) 
      gr.name <- names(c.m) # the names of the innermost groups
      n.g<-length(c.m)   # number of innermost groups
      j0<-rep(1,n.cg) # place holders in V[[i]]'s
      ii <- 1:n
      for (i in 1:n.g) # work through innermost groups
      { # first identify coarse grouping
        Clev <- unique(Cgrps[grps==gr.name[i]])  # level for coarse grouping factor
        if (length(Clev)>1) stop("inner groupings not nested in outer!!")
        k <- (1:n.cg)[Clevel==Clev] # index of coarse group - i.e. update V[[k]] 
        # now need to get c.m into right place within V[[k]]
        j1<-j0[k]+nrow(c.m[[i]])-1
        V[[k]][j0[k]:j1,j0[k]:j1]<-c.m[[i]]
        ind1 <- ii[grps==gr.name[i]] 
        # ind1 is the rows of original data.frame to which c.m[[i]] applies 
        # assuming that data frame order is preserved at the inner grouping
        ind2 <- rCind[ind1] 
        # ind2 contains the rows of the coarse ordering to which c.m[[i]] applies
        ind[[k]][j0[k]:j1] <- ind2 - Voff[k] + 1
        # ind[k] accumulates rows within coarse group k to which V[[k]] applies
        j0[k]<-j1+1  
      }
      for (k in 1:n.cg) { # pasting correlations into right place in each matrix
        V[[k]][ind[[k]],]<-V[[k]];V[[k]][,ind[[k]]]<-V[[k]] 
      }
    }
  } 
  # now form diag(w)%*%V%*%diag(w), depending on class of V
  if (is.list(V)) # it's  a block diagonal structure
  { for (i in 1:n.cg)
    { wi <- w[Voff[i]:(Voff[i]+size.cg[i]-1)] 
      V[[i]] <- as.vector(wi)*t(as.vector(wi)*V[[i]])
    }
  } else
  if (is.matrix(V))
  { V <- as.vector(w)*t(as.vector(w)*V)
  } else # it's a diagonal matrix
  { V <- w^2*V
  }
  # ... covariance matrix according to fitted correlation structure in coarse
  # group order
  
  ## Now work on the random effects ..... 
  X <- list()
  grp.dims <- b$dims$ncol # number of Zt columns for each grouping level (inner levels first)
  # inner levels are first in Zt
  Zt <- model.matrix(b$modelStruct$reStruct,data)  # a sort of proto - Z matrix
  # b$groups and cov (defined below have the inner levels last)
  cov <- as.matrix(b$modelStruct$reStruct) # list of estimated covariance matrices (inner level last)
  i.col<-1
  Z <- matrix(0,n,0) # Z matrix
  if (start.level<=n.levels)
  { for (i in 1:(n.levels-start.level+1)) # work through the r.e. groupings inner to outer
    { # get matrix with columns that are indicator variables for ith set of groups...
      # groups has outer levels first 
      if(length(levels(b$groups[[n.levels-i+1]]))==1) { ## model.matrix needs >1 level 
        X[[1]] <- matrix(rep(1,nrow(b$groups))) } else {
	clist <- list('b$groups[[n.levels - i + 1]]'=c("contr.treatment","contr.treatment"))
        X[[1]] <- model.matrix(~b$groups[[n.levels - i + 1]]-1,
                  contrasts.arg=clist) }
      # Get `model matrix' columns relevant to current grouping level...
      X[[2]] <- Zt[,i.col:(i.col+grp.dims[i]-1),drop=FALSE]
      i.col <- i.col+grp.dims[i]
      # tensor product the X[[1]] and X[[2]] rows...
      Z <- cbind(Z,tensor.prod.model.matrix(X))
    } # so Z assembled from inner to outer levels
    # Now construct overall ranef covariance matrix
    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:b$dims$ngrps[i]) 
      { stop <- start+ncol(cov[[k]])-1
        Vr[start:stop,start:stop]<-cov[[k]]
        start <- stop+1
      }
    }
    Vr <- Vr*b$sigma^2
    ## Now re-order Z into coarse group order
    Z <- Z[Cind,]
    ## Now Z %*% Vr %*% t(Z) is block diagonal: if Z' = [Z1':Z2':Z3': ... ]
    ## where Zi contains th rows of Z for the ith level of the coarsest
    ## grouping factor, then the ith block of (Z Vr Z') is (Zi Vr Zi')
    if (n.cg == 1) { 
      if (is.matrix(V)) { 
        V <- V+Z%*%Vr%*%t(Z)
      } else V <- diag(V) + Z%*%Vr%*%t(Z) 
    } else { # V has a block - diagonal structure
      j0 <- 1
      Vz <- list()
      for (i in 1:n.cg) {
        j1 <- size.cg[i] + j0 -1
        Zi <- Z[j0:j1,,drop=FALSE]
        Vz[[i]] <- Zi %*% Vr %*% t(Zi) 
        j0 <- j1+1
      }
      if (is.list(V)) {
        for (i in 1:n.cg) V[[i]] <- V[[i]]+Vz[[i]] 
      } else { 
        j0 <-1
        for (i in 1:n.cg) {
          kk <- size.cg[i]
          j1 <- kk + j0 -1
          Vz[[i]] <- Vz[[i]] + diag(x=V[j0:j1],nrow=kk,ncol=kk)
          j0 <- j1+1
        }
        V <- Vz
      }
    }
  }
  list(V=V,ind=Cind)
} ## extract.lme.cov2

extract.lme.cov<-function(b,data=NULL,start.level=1)
# function to extract the response data covariance matrix from an lme fitted
# model object b, fitted to the data in data. "inner" == "finest" grouping 
# start.level is the r.e. grouping level at which to start the construction, 
# levels outer to this will not be included in the calculation - this is useful
# for gamm calculations
{ if (!inherits(b,"lme")) stop("object does not appear to be of class lme")
  if (is.null(data)) {
    na.act <- na.action(b)
    data <- if (is.null(na.act)) b$data else b$data[-na.act,] 
  }
  grps<-nlme::getGroups(b) # labels of the innermost groupings - in data frame order
  n<-length(grps)    # number of data
  if (is.null(b$modelStruct$varStruct)) w<-rep(b$sigma,n) ### 
  else 
  { w<-1/nlme::varWeights(b$modelStruct$varStruct) 
    # w is not in data.frame order - it's in inner grouping level order
    group.name<-names(b$groups) # b$groups[[i]] doesn't always retain factor ordering
    order.txt <- paste("ind<-order(data[[\"",group.name[1],"\"]]",sep="")
    if (length(b$groups)>1) for (i in 2:length(b$groups)) 
    order.txt <- paste(order.txt,",data[[\"",group.name[i],"\"]]",sep="")
    order.txt <- paste(order.txt,")")
    eval(parse(text=order.txt))
    w[ind] <- w # into data frame order
    w<-w*b$sigma
  }
  if (is.null(b$modelStruct$corStruct)) V<-diag(n) #*b$sigma^2
  else
  { c.m<-nlme::corMatrix(b$modelStruct$corStruct) # correlation matrices for each group
    if (!is.list(c.m)) V<-c.m
    else
    { V<-matrix(0,n,n)   # data cor matrix
      gr.name <- names(c.m) # the names of the groups
      n.g<-length(c.m)   # number of innermost groups
      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 # pasting correlations into right place in overall matrix
      # V<-V*b$sigma^2
    }
  }  
  V <- as.vector(w)*t(as.vector(w)*V) # diag(w)%*%V%*%diag(w)
  # ... covariance matrix according to fitted correlation structure
  X<-list()
  grp.dims<-b$dims$ncol # number of Zt columns for each grouping level (inner levels first)
  # inner levels are first in Zt
  Zt<-model.matrix(b$modelStruct$reStruct,data)  # a sort of proto - Z matrix
  # b$groups and cov (defined below have the inner levels last)
  cov<-as.matrix(b$modelStruct$reStruct) # list of estimated covariance matrices (inner level last)
  i.col<-1
  n.levels<-length(b$groups)
  Z<-matrix(0,n,0) # Z matrix
  if (start.level<=n.levels)
  { for (i in 1:(n.levels-start.level+1)) # work through the r.e. groupings inner to outer
    { # get matrix with columns that are indicator variables for ith set of groups...
      # groups has outer levels first 
      if(length(levels(b$groups[[n.levels-i+1]]))==1) { ## model.matrix needs >1 level 
        X[[1]] <- matrix(rep(1,nrow(b$groups))) } else {
	clist <- list('b$groups[[n.levels - i + 1]]'=c("contr.treatment","contr.treatment"))
        X[[1]] <- model.matrix(~b$groups[[n.levels - i + 1]]-1,
                  contrasts.arg=clist) }
      # Get `model matrix' columns relevant to current grouping level...
      X[[2]] <- Zt[,i.col:(i.col+grp.dims[i]-1),drop=FALSE]
      i.col <- i.col+grp.dims[i]
      # tensor product the X[[1]] and X[[2]] rows...
      Z <- cbind(Z,tensor.prod.model.matrix(X))
    } # so Z assembled from inner to outer levels
    # Now construct overall ranef covariance matrix
    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:b$dims$ngrps[i]) 
      { stop <- start+ncol(cov[[k]])-1
        Vr[start:stop,start:stop]<-cov[[k]]
        start <- stop+1
      }
    }
    Vr <- Vr*b$sigma^2
    V <- V+Z%*%Vr%*%t(Z)
  }
  V
} ## extract.lme.cov

formXtViX <- function(V,X)
## forms X'V^{-1}X as efficiently as possible given the structure of
## V (diagonal, block-diagonal, full)
## Actually returns R where R'R =  X'V^{-1}X
{ X <- X[V$ind,,drop=FALSE] # have to re-order X according to V ordering
  if (is.list(V$V)) {     ### block diagonal case
    Z <- X
    j0 <- 1
    for (i in 1:length(V$V))
    { Cv <- chol(V$V[[i]])
      j1 <- j0+nrow(V$V[[i]])-1
      Z[j0:j1,] <- backsolve(Cv,X[j0:j1,,drop=FALSE],transpose=TRUE)
      j0 <- j1 + 1
    }
    #res <- t(Z)%*%Z
  } else if (is.matrix(V$V)) { ### full matrix case
    Cv <- chol(V$V)
    Z <- backsolve(Cv,X,transpose=TRUE)
    #res <- t(Z)%*%Z
  } else {                ### diagonal matrix case
    Z <- X/sqrt(as.numeric(V$V))
    #res <- t(X)%*%(X/as.numeric(V$V))
  }
  qrz <- qr(Z,LAPACK=TRUE)
  R <- qr.R(qrz);R[,qrz$pivot] <- R
  colnames(R) <- colnames(X)
  #res <- crossprod(R)
  #res
  R
}


new.name <- function(proposed,old.names)
# finds a name based on proposed, that is not in old.names
# if the proposed name is in old.names then ".xx" is added to it 
# where xx is a number chosen to ensure the a unique name
{ prop <- proposed
  k <- 0
  while (sum(old.names==prop))
  { prop<-paste(proposed,".",k,sep="")
    k <- k + 1
  }
  prop
}

gammPQL <- function (fixed, random, family, data, correlation, weights,
    control, niter = 30, verbose = TRUE, mustart=NULL, etastart=NULL, ...) {
## service routine for `gamm' to do PQL fitting. Based on glmmPQL
## from the MASS library (Venables & Ripley). Because `gamm' already
## does some of the preliminary stuff that glmmPQL does, gammPQL can
## be simpler. It also deals with the possibility of the original
## data frame containing variables called `zz' `wts' or `invwt'.
## Modified 2019 to use standard GLM initialization to improve convergence.
## Old glmmPQL style initialization commented out (single # first col) for
## eventual removal.
  off <- model.offset(data)
  if (is.null(off)) off <- 0
  
  y <- model.response(data) 
  nobs <- nrow(data) 
  if (is.null(weights)) weights <- rep(1, nrow(data))

  start <-  NULL ## never used
  if (is.null(mustart)) {
    eval(family$initialize)
  } else {
    mukeep <- mustart
    eval(family$initialize)
    mustart <- mukeep
  }
  #eval(family$initialize) ## NEW

  wts <- weights
# if (is.null(wts)) wts <- rep(1, nrow(data))
  wts.name <- new.name("wts",names(data)) ## avoid overwriting what's already in `data'
  data[[wts.name]] <- wts 
 
#  fit0 <- NULL ## keep checking tools happy 
  ## initial fit (might be better replaced with `gam' call)
#  eval(parse(text=paste("fit0 <- glm(formula = fixed, family = family, data = data,",
#                        "weights =",wts.name,",...)")))
#  w <- fit0$prior.weights
#  eta <- fit0$linear.predictors
  fam <- family
  eta <- if (!is.null(etastart)) 
            etastart
	 else fam$linkfun(mustart)
  mu <- fam$linkinv(eta)	 
  w <- wts;  
#  zz <- eta + fit0$residuals - off
  mu.eta.val <- fam$mu.eta(eta)
  zz <- eta + (y - mu)/mu.eta.val - off
#  wz <- fit0$weights
  wz <- w * mu.eta.val^2/fam$variance(mustart)
  
  ## find non clashing name for pseudodata and insert in formula
  zz.name <- new.name("zz",names(data))
  eval(parse(text=paste("fixed[[2]] <- quote(",zz.name,")")))
 
  data[[zz.name]] <- zz ## pseudodata to `data' 
  
  ## find non-clashing name for inverse weights, and make 
  ## varFixed formula using it...
  
  invwt.name <- new.name("invwt",names(data))
  data[[invwt.name]] <- 1/wz
  w.formula <- as.formula(paste("~",invwt.name,sep=""))

  converged <- FALSE 

  if (family$family %in% c("poisson","binomial")) {
    control$sigma <- 1 ## set scale parameter to 1
    control$apVar <- FALSE ## not available
  }

  for (i in 1:niter) {
    if (verbose) message(gettextf("iteration %d", i))
    fit <- lme(fixed=fixed,random=random,data=data,correlation=correlation,
             control=control,weights=varFixed(w.formula),method="ML",...)
    etaold <- eta
    eta <- fitted(fit) + off
    if (sum((eta - etaold)^2) < 1e-06 * sum(eta^2)) {
      converged <- TRUE
      break
    }
    mu <- fam$linkinv(eta)
    mu.eta.val <- fam$mu.eta(eta)
    ## get pseudodata and insert in `data' 
#    data[[zz.name]] <- eta + (fit0$y - mu)/mu.eta.val - off
    data[[zz.name]] <- eta + (y - mu)/mu.eta.val - off
    wz <- w * mu.eta.val^2/fam$variance(mu)
    data[[invwt.name]] <- 1/wz
  } ## end i in 1:niter
  if (!converged) warning("gamm not converged, try increasing niterPQL")
# fit$y <- fit0$y
  fit$y <- y 
  fit$w <- w ## prior weights
  fit
} ## gammPQL



gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=list(),weights=NULL,
      subset=NULL,na.action,knots=NULL,control=list(niterEM=0,optimMethod="L-BFGS-B",returnObject=TRUE),
      niterPQL=20,verbosePQL=TRUE,method="ML",drop.unused.levels=TRUE,mustart=NULL, etastart=NULL,...)
# Routine to fit a GAMM to some data. Fixed and smooth terms are defined in the formula, but the wiggly 
# parts of the smooth terms are treated as random effects. The onesided formula random defines additional 
# random terms. correlation describes the correlation structure. This routine is basically an interface
# between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model.
{ if (inherits(family,"extended.family")) warning("family are not designed for use with gamm!")
  ## lmeControl turns sigma=NULL into sigma=0, but if you supply sigma=0 rejects it,
  ## which will break the standard handling of the control list. Following line fixes.
  ## but actually Martin M has now fixed lmeControl...
  ##if (!is.null(control$sigma)&&control$sigma==0) control$sigma <- NULL
  if (inherits(family,"extended.family")) warning("gamm is not designed to use extended families")
  control <- do.call("lmeControl",control) 
    # check that random is a named list
    if (!is.null(random))
    { if (is.list(random)) 
      { r.names<-names(random)
        if (is.null(r.names)) stop("random argument must be a *named* list.")
        else if (sum(r.names=="")) stop("all elements of random list must be named")
      }
      else stop("gamm() can only handle random effects defined as named lists")
      random.vars<-c(unlist(lapply(random, function(x) all.vars(formula(x)))),r.names)
    } else random.vars<-NULL

    if (!is.null(correlation)) { 
      cor.for<-attr(correlation,"formula")
      if (!is.null(cor.for)) cor.vars<-all.vars(cor.for)
    } else cor.vars<-NULL

    ## now establish whether weights is varFunc or not...  
    wisvf <- try(inherits(weights,"varFunc"),silent=TRUE)
    if (inherits(wisvf,"try-error")) wisvf <- FALSE
    if (wisvf) { ## collect its variables
      if (inherits(weights,"varComb")) { ## actually a list of varFuncs
        vf.vars <- rep("",0)
        for (i in 1:length(weights)) {
          vf.vars <- c(vf.vars,all.vars(attr(weights[[i]],"formula")))
        }
        vf.vars <- unique(vf.vars)
      } else { ## single varFunc
        vf.for<-attr(weights,"formula")
        if (!is.null(vf.for)) vf.vars<-all.vars(vf.for)
      }
    } else vf.vars <- NULL

    # create model frame.....
    gp <- interpret.gam(formula) # interpret the formula 
    ##cl<-match.call() # call needed in gamm object for update to work
    mf <- match.call(expand.dots=FALSE)
    mf$formula <- gp$fake.formula
    if (wisvf) {
      mf$correlation <- mf$random <- mf$family <- mf$control <- mf$scale <- mf$knots <- mf$sp <- mf$weights <-
      mf$min.sp <- mf$H <- mf$gamma <- mf$fit <- mf$niterPQL <- mf$verbosePQL <- mf$G <- mf$method <- mf$... <- NULL
    } else {
      mf$correlation <- mf$random <- mf$family <- mf$control <- mf$scale <- mf$knots <- mf$sp <-
      mf$min.sp <- mf$H <- mf$gamma <- mf$fit <- mf$niterPQL <- mf$verbosePQL <- mf$G <- mf$method <- mf$... <- NULL
    }
    mf$drop.unused.levels <- drop.unused.levels
    mf[[1]] <- quote(stats::model.frame) ## as.name("model.frame")
    pmf <- mf
    gmf <- eval(mf, parent.frame()) # the model frame now contains all the data, for the gam part only 
    gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly

    if (!wisvf) weights <- gmf[["(weights)"]]

    allvars <- c(cor.vars,random.vars,vf.vars)
    if (length(allvars)) {
      mf$formula <- as.formula(paste(paste(deparse(gp$fake.formula,backtick=TRUE),collapse=""),
                           "+",paste(allvars,collapse="+")))
      mf <- eval(mf, parent.frame()) # the model frame now contains *all* the data
    }
    else mf <- gmf
    rm(gmf)
    if (nrow(mf)<2) stop("Not enough (non-NA) data to do anything meaningful")
    ##Terms <- attr(mf,"terms")    
  
    ## summarize the *raw* input variables
    ## note can't use get_all_vars here -- buggy with matrices
    vars <- all_vars1(gp$fake.formula[-2]) ## drop response here
    inp <- parse(text = paste("list(", paste(vars, collapse = ","),")"))
    dl <- eval(inp, data, parent.frame())
    names(dl) <- vars ## list of all variables needed
    var.summary <- variable.summary(gp$pf,dl,nrow(mf)) ## summarize the input data
    rm(dl) ## save space 

    pmf$formula <- gp$pf
    pmf <- eval(pmf, parent.frame()) # pmf contains all data for parametric part 
    pTerms <- attr(pmf,"terms")

    if (is.character(family)) family<-eval(parse(text=family))
    if (is.function(family)) family <- family()
    if (is.null(family$family)) stop("family not recognized")
  
    # now call gamm.setup 

    G <- gamm.setup(gp,pterms=pTerms,
                    data=mf,knots=knots,parametric.only=FALSE,absorb.cons=TRUE)
    #G$pterms <- pTerms
    G$var.summary <- var.summary    
    mf <- G$data

    n.sr <- length(G$random) # number of random smooths (i.e. s(...,fx=FALSE,...) terms)

    if (is.null(random)&&n.sr==0) 
    stop("gamm models must have at least 1 smooth with unknown smoothing parameter or at least one other random effect")

    offset.name <- attr(mf,"names")[attr(attr(mf,"terms"),"offset")]

    yname <- new.name("y",names(mf))
    eval(parse(text=paste("mf$",yname,"<-G$y",sep="")))
    Xname <- new.name("X",names(mf))
    eval(parse(text=paste("mf$",Xname,"<-G$X",sep="")))
    
    fixed.formula <- paste(yname,"~",Xname,"-1")
    ## following appears to serve no purpose except confusing later lme versions
    #if (length(offset.name)) {
    #  fixed.formula <- paste(fixed.formula,"+",offset.name) 
    #}
    fixed.formula <- as.formula(fixed.formula)
    
    ## Add any explicit random effects to the smooth induced r.e.s 
    rand <- G$random  
    if (!is.null(random))
    { r.m <- length(random)
      r.names <- c(names(rand),names(random))
      for (i in 1:r.m) rand[[n.sr+i]]<-random[[i]]   
      names(rand) <- r.names
    }

    ## need to modify the correlation structure formula, in order that any
    ## grouping factors for correlation get nested within at least the 
    ## constructed dummy grouping factors.

    if (length(formula(correlation))) # then modify the correlation formula
    { # first get the existing grouping structure ....
      corGroup <- paste(names(rand),collapse="/")
      groupForm<-nlme::getGroupsFormula(correlation)
      if (!is.null(groupForm)) {
        groupFormNames <- all.vars(groupForm)
        exind <- groupFormNames %in% names(rand)
        groupFormNames <- groupFormNames[!exind] ## dumping duplicates 
        if (length(groupFormNames)) corGroup <- 
             paste(corGroup,paste(groupFormNames,collapse="/"),sep="/")
      }
      # now make a new formula for the correlation structure including these groups
      corForm <- as.formula(paste(deparse(nlme::getCovariateFormula(correlation)),"|",corGroup))
      attr(correlation,"formula") <- corForm
    }

    ### Actually do fitting ....
 
    ret <- list()

    if (family$family=="gaussian"&&family$link=="identity"&&
    length(offset.name)==0) lme.used <- TRUE else lme.used <- FALSE
    if (lme.used&&!is.null(weights)&&!wisvf) lme.used <- FALSE   

    if (lme.used) { ## following construction is a work-around for problem in nlme 3-1.52 
      eval(parse(text=paste("ret$lme<-lme(",deparse(fixed.formula),
          ",random=rand,data=strip.offset(mf),correlation=correlation,",
          "control=control,weights=weights,method=method)"
            ,sep=""    )))
      ## need to be able to work out full edf for following to work...	    
      # if (is.null(correlation)) { ## compute conditional aic precursor
      #  dev <- sum(family$dev.resids(G$y,fitted(ret$lme),weights))
      #	ret$lme$aic <- family$aic(G$y,1,fitted(ret$lme),weights,dev)
      # }
      ##ret$lme<-lme(fixed.formula,random=rand,data=mf,correlation=correlation,control=control)
    } else { ## Again, construction is a work around for nlme 3-1.52
      if (wisvf) stop("weights must be like glm weights for generalized case")
      if (verbosePQL) cat("\n Maximum number of PQL iterations: ",niterPQL,"\n")
      eval(parse(text=paste("ret$lme<-gammPQL(",deparse(fixed.formula),
          ",random=rand,data=strip.offset(mf),family=family,",
          "correlation=correlation,control=control,",
            "weights=weights,niter=niterPQL,verbose=verbosePQL,mustart=mustart,etastart=etastart,...)",sep=""))) 
      G$y <- ret$lme$y ## makes sure that binomial response is returned as a vector!
     
    }

    ### .... fitting finished

    # now fake a gam object 
    
    object <- list(model=mf,formula=formula,smooth=G$smooth,nsdf=G$nsdf,family=family,
                 df.null=nrow(G$X),y=G$y,terms=gam.terms,pterms=G$pterms,xlevels=G$xlevels,
                 contrasts=G$contrasts,assign=G$assign,na.action=attr(mf,"na.action"),
                 cmX=G$cmX,var.summary=G$var.summary,scale.estimated=TRUE)
    
    pvars <- all.vars(delete.response(object$terms))
    object$pred.formula <- if (length(pvars)>0) reformulate(pvars) else NULL

    #######################################################
    ## Transform  parameters back to the original space....
    #######################################################

    bf <- as.numeric(ret$lme$coefficients$fixed)

#    br <- as.numeric(unlist(ret$lme$coefficients$random))

    ## Grouped random coefs are returned in matrices. Each row relates to one 
    ## level of the grouping factor. So to get coefs in order, with all coefs
    ## for each grouping factor level contiguous, requires the following...
 
    br <- as.numeric(unlist(lapply(ret$lme$coefficients$random,t)))

    fs.present <- FALSE

    if (G$nsdf) p <- bf[1:G$nsdf] else p <- array(0,0)
    if (G$m>0) for (i in 1:G$m)
    { fx <- G$smooth[[i]]$fixed 
      first <- G$smooth[[i]]$first.f.para;last <- G$smooth[[i]]$last.f.para
      if (first <=last) beta <- bf[first:last] else beta <- array(0,0) ## fixed coefs
      if (fx) p <- c(p, beta) else { ## not fixed so need to undo transform of random effects etc. 
        ind <- G$smooth[[i]]$rind ##G$smooth[[i]]$first.r.para:G$smooth[[i]]$last.r.para
        if (!is.null(G$smooth[[i]]$fac)) { ## nested term, need to unpack coefs at each level of fac
          fs.present <- TRUE
          if (first<=last) stop("Nested smooths must be fully random")
          flev <- levels(G$smooth[[i]]$fac)
          for (j in 1:length(flev)) {
            b <- br[ind] 
            b <- G$smooth[[i]]$trans.D*b
            if (!is.null(G$smooth[[i]]$trans.U)) b <- G$smooth[[i]]$trans.U%*%b
            ind <- ind + G$smooth[[i]]$rinc
            p <- c(p,b)
          }
        } else { ## single level
          b <- c(br[ind],beta)
          b <- G$smooth[[i]]$trans.D*b
          if (!is.null(G$smooth[[i]]$trans.U)) b <- G$smooth[[i]]$trans.U%*%b 
          p <- c(p,b)
        }
      }
    }
 
    var.param <- coef(ret$lme$modelStruct$reStruct)
    n.v <- length(var.param) 
    # k <- 1
    spl <- list()
    if (G$m>0) for (i in 1:G$m) # var.param in reverse term order, but forward order within terms!!
    { ii <- G$pord[i]
      n.sp <- length(object$smooth[[ii]]$S) # number of s.p.s for this term 
      if (n.sp>0) {
        if (inherits(object$smooth[[ii]],"tensor.smooth")) ## ... really mean pdTens used here...
        ## object$sp[k:(k+n.sp-1)] 
          spl[[ii]] <- notExp2(var.param[(n.v-n.sp+1):n.v])
        else ## object$sp[k:(k+n.sp-1)] 
          spl[[ii]] <- 1/notExp2(var.param[n.v:(n.v-n.sp+1)])  
      }
      # k <- k + n.sp
      n.v <- n.v - n.sp
    }
    
    object$sp <- rep(0,0)
    if (length(spl)) for (i in 1:length(spl)) if (!is.null(spl[[i]])) object$sp <- c(object$sp,spl[[i]]) 
    if (length(object$sp)==0) object$sp <- NULL  

    object$coefficients <- p
    
    V <- extract.lme.cov2(ret$lme,mf,n.sr+1) # the data covariance matrix, excluding smooths
 

    ## obtain XVX and S....
    first.para <- last.para <- rep(0,G$m) ## collect first and last para info relevant to expanded Xf
    if (fs.present) { ## First create expanded Xf...
      Xf <- G$Xf[,1:G$nsdf,drop=FALSE] 
      if (G$m>0) for (i in 1:G$m) {# Accumulate the total model matrix
        ind <- object$smooth[[i]]$first.para:object$smooth[[i]]$last.para
        if (is.null(object$smooth[[i]]$fac)) { ## normal smooth
           first.para[i] <- ncol(Xf)+1
           Xf <- cbind(Xf,G$Xf[,ind])
           last.para[i] <- ncol(Xf)
        } else { ## smooth conditioned on multilevel factor. Expand Xf
          flev <- levels(object$smooth[[i]]$fac)
          first.para[i] <- ncol(Xf)+1
          for (k in 1:length(flev)) {
            Xf <- cbind(Xf,G$Xf[,ind]*as.numeric(object$smooth[[i]]$fac==flev[k]))
          }
          last.para[i] <- ncol(Xf)
        }
      }

      object$R <- formXtViX(V,Xf) ## inefficient, if there are smooths conditioned on factors
      XVX <- crossprod(object$R)
      nxf <- ncol(Xf)
    } else {
      if (G$m>0) for (i in 1:G$m) {
        first.para[i] <- object$smooth[[i]]$first.para
        last.para[i] <- object$smooth[[i]]$last.para
      }
      object$R <- formXtViX(V,G$Xf)
      XVX <- crossprod(object$R)
      nxf <- ncol(G$Xf)
    }

    object$R <- object$R*ret$lme$sigma ## correction to what is required by summary.gam

    ## Now S...
    S <- matrix(0,nxf,nxf) ## penalty matrix
    first <- G$nsdf+1
    k <- 1
    if (G$m>0) for (i in 1:G$m) {# Accumulate the total penalty matrix 
      if (is.null(object$smooth[[i]]$fac)) { ## simple regular smooth...
        ind <- first.para[i]:last.para[i]
        ns <-length(object$smooth[[i]]$S) 
        if (ns) for (l in 1:ns) { 
           S[ind,ind] <- S[ind,ind] + 
                  object$smooth[[i]]$S[[l]]*object$sp[k]
           k <- k+1
        }
      } else { ## smooth conditioned on factor
        flev <- levels(object$smooth[[i]]$fac)
        ind <- first.para[i]:(first.para[i]+object$smooth[[i]]$n.para-1)
        ns <- length(object$smooth[[i]]$S)
        for (j in 1:length(flev)) {
           if (ns) for (l in 1:ns) { 
             S[ind,ind] <- S[ind,ind] + 
                    object$smooth[[i]]$S[[l]]*object$sp[k]
             k <- k+1
           }
           k <- k - ns ## same smoothing parameters repeated for each factor level
           ind <- ind + object$smooth[[i]]$n.para
        }
        k <- k + ns
      }
    }
    S <- S/ret$lme$sigma^2 # X'V^{-1}X divided by \sigma^2, so should S be
    
    ## now replace original first.para and last.para with expanded versions...
    if (G$m) for (i in 1:G$m) { 
      object$smooth[[i]]$first.para <- first.para[i]
      object$smooth[[i]]$last.para <- last.para[i] 
    }
    ## stable computation of coefficient covariance matrix...
    ev <- eigen(XVX+S,symmetric=TRUE)
    ind <- ev$values != 0
    iv <- ev$values;iv[ind] <- 1/ev$values[ind]
    Vb <- ev$vectors%*%(iv*t(ev$vectors))    

    object$edf<-rowSums(Vb*t(XVX))   
    object$df.residual <- length(object$y) - sum(object$edf)    

    #if (!is.null(ret$lme$aic)) { ## finish the conditional AIC (only happens if no correlation) 
    #  object$aic <- ret$lme$aic + sum(object$edf) ## requires edf for r.e. as well!
    #  ret$lme$aic<- NULL
    #}

    object$sig2 <- ret$lme$sigma^2
    if (lme.used) { object$method <- paste("lme.",method,sep="")} 
    else { object$method <- "PQL"}

    if (!lme.used||method=="ML") Vb <- Vb*length(G$y)/(length(G$y)-G$nsdf)
    object$Vp <- Vb
    object$Ve <- Vb%*%XVX%*%Vb
    
    object$prior.weights <- weights
    class(object) <- "gam"

    ## If prediction parameterization differs from fit parameterization, transform now...
    ## (important for t2 smooths, where fit constraint is not good for component wise 
    ##  prediction s.e.s)

    if (!is.null(G$P)) {
      object$coefficients <- G$P %*% object$coefficients
      object$Vp <- G$P %*% object$Vp %*% t(G$P) 
      object$Ve <- G$P %*% object$Ve %*% t(G$P) 
    }

    object$linear.predictors <- predict.gam(object,type="link")
    object$fitted.values <- object$family$linkinv(object$linear.predictors)  
 
    object$residuals <- residuals(ret$lme) #as.numeric(G$y) - object$fitted.values

    if (G$nsdf>0) term.names<-colnames(G$X)[1:G$nsdf] else term.names<-array("",0)
    n.smooth <- length(G$smooth) 
    if (n.smooth) {
      for (i in 1:n.smooth)
      { k <- 1
        for (j in object$smooth[[i]]$first.para:object$smooth[[i]]$last.para)
        { term.names[j] <- paste(object$smooth[[i]]$label,".",as.character(k),sep="")
          k <- k+1
        }
      }
      if (!is.null(object$sp)) names(object$sp) <- names(G$sp)
    }

    names(object$coefficients) <- term.names  # note - won't work on matrices!!
    names(object$edf) <- term.names
    if (is.null(weights)) object$prior.weights <- object$y*0+1
    else if (wisvf) object$prior.weights <- varWeights.dfo(ret$lme,mf)^2
    else object$prior.weights <- ret$lme$w
    
    object$weights <- object$prior.weights   

    if (!is.null(G$Xcentre)) object$Xcentre <- G$Xcentre ## column centering values
    ## set environments to global to avoid enormous saved object files
    environment(attr(object$model,"terms")) <- 
    environment(object$terms) <- environment(object$pterms) <- 
    environment(object$formula) <- .GlobalEnv
    if (!is.null(object$pred.formula)) environment(object$pred.formula) <-  .GlobalEnv
    ret$gam <- object
    environment(attr(ret$lme$data,"terms")) <- environment(ret$lme$terms) <- .GlobalEnv
    if (!is.null(ret$lme$modelStruct$varStruct)) {
      environment(attr(ret$lme$modelStruct$varStruct,"formula")) <- .GlobalEnv
    }  
    if (!is.null(ret$lme$modelStruct$corStruct)) {
      environment(attr(ret$lme$modelStruct$corStruct,"formula")) <- .GlobalEnv
    }
    class(ret) <- c("gamm","list")
    ret

} ## end gamm



test.gamm <- function(control=nlme::lmeControl(niterEM=3,tolerance=1e-11,msTol=1e-11))
## this function is a response to repeated problems with nlme/R updates breaking
## the pdTens class. It tests for obvious breakages!
{ test1<-function(x,z,sx=0.3,sz=0.4)
  { x<-x*20
    (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
    0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
  }
  compare <- function(b,b1,edf.tol=.001) 
  { edf.diff <- abs(sum(b$edf)-sum(b1$edf))
    fit.cor <- cor(fitted(b),fitted(b1))
    if (fit.cor<.999) { cat("FAILED: fit.cor = ",fit.cor,"\n");return()}
    if (edf.diff>edf.tol) { cat("FAILED: edf.diff = ",edf.diff,"\n");return()}
    cat("PASSED \n")
  }
  n<-500
  x<-runif(n)/20;z<-runif(n);
  f <- test1(x,z)
  y <- f + rnorm(n)*0.2
  control$sigma <- NULL ## avoid failure on silly test
  cat("testing covariate scale invariance ... ")  
  b <- gamm(y~te(x,z), control=control )
  x1 <- x*100 
  b1 <- gamm(y~te(x1,z),control=control)
  res <- compare(b$gam,b1$gam)
 
  cat("testing invariance w.r.t. response ... ")
  y1 <- y*100 
  b1 <- gamm(y1~te(x,z),control=control)
  res <- compare(b$gam,b1$gam)
  
  cat("testing equivalence of te(x) and s(x) ... ")
  b2 <- gamm(y~te(x,k=10,bs="cr"),control=control)
  b1 <- gamm(y~s(x,bs="cr",k=10),control=control)
  res <- compare(b2$gam,b1$gam,edf.tol=.1)

  cat("testing equivalence of gam and gamm with same sp ... ")
  b1 <- gam(y~te(x,z),sp=b$gam$sp)
  res <- compare(b$gam,b1)  
  
  if (FALSE) cat(res,x1,y1) ## fool codetools
}