mgcv.r 180 KB
 Dirk Eddelbuettel committed Apr 10, 2018 1 ## R routines for the package mgcv (c) Simon Wood 2000-2013  Dirk Eddelbuettel committed Apr 10, 2018 2 ## With contributions from Henric Nilsson  Dirk Eddelbuettel committed Apr 10, 2018 3   Dirk Eddelbuettel committed Apr 10, 2018 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Rrank <- function(R,tol=.Machine$double.eps^.9) { ## Finds rank of upper triangular matrix R, by estimating condition ## number of upper rank by rank block, and reducing rank until this is ## acceptably low... assumes R pivoted rank <- m <- ncol(R) ok <- TRUE while (ok) { Rcond <- .C(C_R_cond,R=as.double(R),r=as.integer(m),c=as.integer(rank), work=as.double(rep(0,4*m)),Rcond=as.double(1))$Rcond if (Rcond*tol<1) ok <- FALSE else rank <- rank - 1 } rank }  Dirk Eddelbuettel committed Apr 10, 2018 18 slanczos <- function(A,k=10,kl=-1,tol=.Machine$double.eps^.5,nt=1) {  Dirk Eddelbuettel committed Apr 10, 2018 19 20 21 22 23 24 25 26 27 28 29 ## Computes truncated eigen decomposition of symmetric A by ## Lanczos iteration. If kl < 0 then k largest magnitude ## eigenvalues returned, otherwise k highest and kl lowest. ## Eigenvectors are always returned too. ## set.seed(1);n <- 1000;A <- matrix(runif(n*n),n,n);A <- A+t(A);er <- slanczos(A,10) ## um <- eigen(A,symmetric=TRUE);ind <- c(1:5,(n-5+1):n) ## range(er$values-um$values[ind]);range(abs(er$vectors)-abs(um$vectors[,ind])) ## It seems that when k (or k+kl) is beyond 10-15% of n ## then you might as well use eigen(A,symmetric=TRUE), but the ## extra cost is the expensive accumulation of eigenvectors. ## Should re-write whole thing using LAPACK routines for eigenvectors.  Dirk Eddelbuettel committed Apr 10, 2018 30  if (tol<=0||tol>.01) stop("silly tolerance supplied")  Dirk Eddelbuettel committed Apr 10, 2018 31  k <- round(k);kl <- round(kl)  Dirk Eddelbuettel committed Apr 10, 2018 32 33  if (k<0) stop("argument k must be positive.") m <- k + max(0,kl)  Dirk Eddelbuettel committed Apr 10, 2018 34  n <- nrow(A)  Dirk Eddelbuettel committed Apr 10, 2018 35  if (m<1) return(list(values=rep(0,0),vectors=matrix(0,n,0),iter=0))  Dirk Eddelbuettel committed Apr 10, 2018 36  if (n != ncol(A)) stop("A not square")  Dirk Eddelbuettel committed Apr 10, 2018 37  if (m>n) stop("Can not have more eigenvalues than nrow(A)")  Dirk Eddelbuettel committed Apr 10, 2018 38  oo <- .C(C_Rlanczos,A=as.double(A),U=as.double(rep(0,n*m)),D=as.double(rep(0,m)),  Dirk Eddelbuettel committed Apr 10, 2018 39  n=as.integer(n),m=as.integer(k),ml=as.integer(kl),tol=as.double(tol),nt=as.integer(nt))  Dirk Eddelbuettel committed Apr 10, 2018 40 41 42  list(values = oo$D,vectors = matrix(oo$U,n,m),iter=oo$n) }  Dirk Eddelbuettel committed Apr 10, 2018 43 44 45 46 47 48 49 50 51 52 53 54 rig <- function(n,mean,scale) { ## inverse guassian deviates generated by algorithm 5.7 of ## Gentle, 2003. scale = 1/lambda. if (length(n)>1) n <- length(n) y <- rnorm(n)^2 mu2 <- 0*y + mean^2 ## y is there to ensure mu2 is a vector x <- mean + 0.5*scale*(mu2*y - mean*sqrt(4*mean*y/scale + mu2*y^2)) ind <- runif(n) > mean/(mean+x) x[ind] <- mu2[ind]/x[ind] x ## E(x) = mean; var(x) = scale*mean^3 }  Dirk Eddelbuettel committed Apr 10, 2018 55 56 57 58 59 60 61 62 63 64 65 66 strip.offset <- function(x) # sole purpose is to take a model frame and rename any "offset(a.name)" # columns "a.name" { na <- names(x) for (i in 1:length(na)) { if (substr(na[i],1,7)=="offset(") na[i] <- substr(na[i],8,nchar(na[i])-1) } names(x) <- na x }  Dirk Eddelbuettel committed Apr 10, 2018 67   Dirk Eddelbuettel committed Apr 10, 2018 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 pcls <- function(M) # Function to perform penalized constrained least squares. # Problem to be solved is: # # minimise ||W^0.5 (y - Xp)||^2 + p'Bp # subject to Ain p >= b & C p = "constant" # # where B = \sum_{i=1}^m \theta_i S_i and W=diag(w) # on entry this routine requires a list M, with the following elements: # M$X - the design matrix for this problem. # M$p - a feasible initial parameter vector - note that this should not be set up to # lie exactly on all the inequality constraints - which can easily happen if M$p=0! # M$y - response variable # M$w - weight vector: W= diag(M$w) # M$Ain - matrix of inequality constraints # M$bin - b above # M$C - fixed constraint matrix # M$S - List of (minimal) penalty matrices # M$off - used for unpacking M$S # M$sp - array of theta_i's # Ain, bin and p are not in the object needed to call mgcv.... # { nar<-c(length(M$y),length(M$p),dim(M$Ain)[1],dim(M$C)[1],0) H<-0 ## sanity checking ... if (nrow(M$X)!=nar[1]) stop("nrow(M$X) != length(M$y)") if (ncol(M$X)!=nar[2]) stop("ncol(M$X) != length(M$p)") if (length(M$w)!=nar[1]) stop("length(M$w) != length(M$y)") if (nar[3]!=length(M$bin)) stop("nrow(M$Ain) != length(M$bin)")  Dirk Eddelbuettel committed Apr 10, 2018 97 98  if (nrow(M$Ain)>0) { if (ncol(M$Ain)!=nar[2]) stop("nrow(M$Ain) != length(M$p)")  Dirk Eddelbuettel committed Apr 10, 2018 99  res <- as.numeric(M$Ain%*%M$p) - as.numeric(M$bin)  Dirk Eddelbuettel committed Apr 10, 2018 100 101 102 103 104  if (sum(res<0)>0) stop("initial parameters not feasible") res <- abs(res) if (sum(res<.Machine$double.eps^.5)>0) warning("initial point very close to some inequality constraints") res <- mean(res)  Dirk Eddelbuettel committed Apr 10, 2018 105  if (res<.Machine$double.eps^.5)  Dirk Eddelbuettel committed Apr 10, 2018 106  warning("initial parameters very close to inequality constraints")  Dirk Eddelbuettel committed Apr 10, 2018 107 108 109 110 111 112 113 114 115 116 117 118 119  } if (nrow(M$C)>0) if (ncol(M$C)!=nar[2]) stop("ncol(M$C) != length(M$p)") if (length(M$S)!=length(M$off)) stop("M$S and M$off have different lengths") if (length(M$S)!=length(M$sp)) stop("M$sp has different length to M$S and M$off") # pack the S array for mgcv call m<-length(M$S) Sa<-array(0,0);df<-0 if (m>0) for (i in 1:m) { Sa<-c(Sa,M$S[[i]]) df[i]<-nrow(M$S[[i]])  Dirk Eddelbuettel committed Apr 10, 2018 120  if (M$off[i]+df[i]-1>nar[2]) stop(gettextf("M$S[%d] is too large given M$off[%d]", i, i))  Dirk Eddelbuettel committed Apr 10, 2018 121  }  Dirk Eddelbuettel committed Apr 10, 2018 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137  qra.exist <- FALSE if (ncol(M$X)>nrow(M$X)) { if (m>0) stop("Penalized model matrix must have no more columns than rows") else { ## absorb M$C constraints qra <- qr(t(M$C)) j <- nrow(M$C);k <- ncol(M$X) M$X <- t(qr.qty(qra,t(M$X))[(j+1):k,]) M$Ain <- t(qr.qty(qra,t(M$Ain))[(j+1):k,]) M$C <- matrix(0,0,0) M$p <- rep(0,ncol(M$X)) nar[2] <- length(M$p) nar[4] <- 0 qra.exist <- TRUE if (ncol(M$X)>nrow(M$X)) stop("Model matrix not full column rank") } }  Dirk Eddelbuettel committed Apr 10, 2018 138  o<-.C(C_RPCLS,as.double(M$X),as.double(M$p),as.double(M$y),as.double(M$w),as.double(M$Ain),as.double(M$bin)  Dirk Eddelbuettel committed Apr 10, 2018 139  ,as.double(M$C),as.double(H),as.double(Sa),as.integer(M$off),as.integer(df),as.double(M$sp),  Dirk Eddelbuettel committed Apr 10, 2018 140  as.integer(length(M$off)),as.integer(nar))  Dirk Eddelbuettel committed Apr 10, 2018 141 142 143  p <- array(o[[2]],length(M$p)) if (qra.exist) p <- qr.qy(qra,c(rep(0,j),p)) p  Dirk Eddelbuettel committed Apr 10, 2018 144 } ## pcls  Dirk Eddelbuettel committed Apr 10, 2018 145 146   Dirk Eddelbuettel committed Apr 10, 2018 147 interpret.gam0 <- function(gf,textra=NULL)  Dirk Eddelbuettel committed Apr 10, 2018 148 149 150 151 152 # interprets a gam formula of the generic form: # y~x0+x1+x3*x4 + s(x5)+ s(x6,x7) .... # and returns: # 1. a model formula for the parametric part: pf (and pfok indicating whether it has terms) # 2. a list of descriptors for the smooths: smooth.spec  Dirk Eddelbuettel committed Apr 10, 2018 153 154 155 # this is function does the work, and is called by in interpret.gam { p.env <- environment(gf) # environment of formula tf <- terms.formula(gf,specials=c("s","te","ti","t2")) # specials attribute indicates which terms are smooth  Dirk Eddelbuettel committed Apr 10, 2018 156   Dirk Eddelbuettel committed Apr 10, 2018 157 158  terms <- attr(tf,"term.labels") # labels of the model terms nt <- length(terms) # how many terms?  Dirk Eddelbuettel committed Apr 10, 2018 159   Dirk Eddelbuettel committed Apr 10, 2018 160 161 162 163  if (attr(tf,"response") > 0) { # start the replacement formulae response <- as.character(attr(tf,"variables")[2]) } else { response <- NULL  Dirk Eddelbuettel committed Apr 10, 2018 164  }  Dirk Eddelbuettel committed Apr 10, 2018 165 166  sp <- attr(tf,"specials")$s # array of indices of smooth terms tp <- attr(tf,"specials")$te # indices of tensor product terms  Dirk Eddelbuettel committed Apr 10, 2018 167  tip <- attr(tf,"specials")$ti # indices of tensor product pure interaction terms  Dirk Eddelbuettel committed Apr 10, 2018 168  t2p <- attr(tf,"specials")$t2 # indices of type 2 tensor product terms  Dirk Eddelbuettel committed Apr 10, 2018 169 170 171 172  off <- attr(tf,"offset") # location of offset in formula ## have to translate sp, tp, tip, t2p so that they relate to terms, ## rather than elements of the formula...  Dirk Eddelbuettel committed Apr 10, 2018 173 174 175 176 177 178 179 180  vtab <- attr(tf,"factors") # cross tabulation of vars to terms if (length(sp)>0) for (i in 1:length(sp)) { ind <- (1:nt)[as.logical(vtab[sp[i],])] sp[i] <- ind # the term that smooth relates to } if (length(tp)>0) for (i in 1:length(tp)) { ind <- (1:nt)[as.logical(vtab[tp[i],])] tp[i] <- ind # the term that smooth relates to  Dirk Eddelbuettel committed Apr 10, 2018 181  }  Dirk Eddelbuettel committed Apr 10, 2018 182 183 184 185  if (length(tip)>0) for (i in 1:length(tip)) { ind <- (1:nt)[as.logical(vtab[tip[i],])] tip[i] <- ind # the term that smooth relates to }  Dirk Eddelbuettel committed Apr 10, 2018 186  if (length(t2p)>0) for (i in 1:length(t2p)) {  Dirk Eddelbuettel committed Apr 10, 2018 187 188  ind <- (1:nt)[as.logical(vtab[t2p[i],])] t2p[i] <- ind # the term that smooth relates to  Dirk Eddelbuettel committed Apr 10, 2018 189  } ## re-referencing is complete  Dirk Eddelbuettel committed Apr 10, 2018 190   Dirk Eddelbuettel committed Apr 10, 2018 191  k <- kt <- kti <- kt2 <- ks <- kp <- 1 # counters for terms in the 2 formulae  Dirk Eddelbuettel committed Apr 10, 2018 192 193  len.sp <- length(sp) len.tp <- length(tp)  Dirk Eddelbuettel committed Apr 10, 2018 194  len.tip <- length(tip)  Dirk Eddelbuettel committed Apr 10, 2018 195  len.t2p <- length(t2p)  Dirk Eddelbuettel committed Apr 10, 2018 196  ns <- len.sp + len.tp + len.tip + len.t2p # number of smooths  Dirk Eddelbuettel committed Apr 10, 2018 197  pav <- av <- rep("",0)  Dirk Eddelbuettel committed Apr 10, 2018 198  smooth.spec <- list()  Dirk Eddelbuettel committed Apr 10, 2018 199  mgcvat <- "package:mgcv" %in% search() ## is mgcv in search path?  Dirk Eddelbuettel committed Apr 10, 2018 200 201 202  if (nt) for (i in 1:nt) { # work through all terms if (k <= ns&&((ks<=len.sp&&sp[ks]==i)||(kt<=len.tp&&tp[kt]==i)|| (kti<=len.tip&&tip[kti]==i)||(kt2<=len.t2p&&t2p[kt2]==i))) { # it's a smooth  Dirk Eddelbuettel committed Apr 10, 2018 203 204 205 206 207 208 209 210 211 212 213  ## have to evaluate in the environment of the formula or you can't find variables ## supplied as smooth arguments, e.g. k <- 5;gam(y~s(x,k=k)), fails, ## but if you don't specify namespace of mgcv then stuff like ## loadNamespace('mgcv'); k <- 10; mgcv::interpret.gam(y~s(x,k=k)) fails (can't find s) ## eval(parse(text=terms[i]),envir=p.env,enclos=loadNamespace('mgcv')) fails?? ## following may supply namespace of mgcv explicitly if not on search path... if (mgcvat) st <- eval(parse(text=terms[i]),envir=p.env) else { st <- try(eval(parse(text=terms[i]),envir=p.env),silent=TRUE) if (inherits(st,"try-error")) st <- eval(parse(text=terms[i]),enclos=p.env,envir=loadNamespace('mgcv')) }  Dirk Eddelbuettel committed Apr 10, 2018 214 215 216 217 218  if (!is.null(textra)) { ## modify the labels on smooths with textra pos <- regexpr("(",st$lab,fixed=TRUE)[1] st$label <- paste(substr(st$label,start=1,stop=pos-1),textra, substr(st$label,start=pos,stop=nchar(st$label)),sep="") }  Dirk Eddelbuettel committed Apr 10, 2018 219  smooth.spec[[k]] <- st  Dirk Eddelbuettel committed Apr 10, 2018 220 221  if (ks<=len.sp&&sp[ks]==i) ks <- ks + 1 else # counts s() terms if (kt<=len.tp&&tp[kt]==i) kt <- kt + 1 else # counts te() terms  Dirk Eddelbuettel committed Apr 10, 2018 222  if (kti<=len.tip&&tip[kti]==i) kti <- kti + 1 else # counts ti() terms  Dirk Eddelbuettel committed Apr 10, 2018 223  kt2 <- kt2 + 1 # counts t2() terms  Dirk Eddelbuettel committed Apr 10, 2018 224  k <- k + 1 # counts smooth terms  Dirk Eddelbuettel committed Apr 10, 2018 225 226  } else { # parametric av[kp] <- terms[i] ## element kp on rhs of parametric  Dirk Eddelbuettel committed Apr 10, 2018 227  kp <- kp+1 # counts parametric terms  Dirk Eddelbuettel committed Apr 10, 2018 228 229  } }  Dirk Eddelbuettel committed Apr 10, 2018 230  if (!is.null(off)) { ## deal with offset  Dirk Eddelbuettel committed Apr 10, 2018 231  av[kp] <- as.character(attr(tf,"variables")[1+off])  Dirk Eddelbuettel committed Apr 10, 2018 232 233  kp <- kp+1 }  Dirk Eddelbuettel committed Apr 10, 2018 234 235  pf <- paste(response,"~",paste(av,collapse=" + "))  Dirk Eddelbuettel committed Apr 10, 2018 236 237 238  if (attr(tf,"intercept")==0) { pf <- paste(pf,"-1",sep="") if (kp>1) pfok <- 1 else pfok <- 0  Dirk Eddelbuettel committed Apr 10, 2018 239  } else {  Dirk Eddelbuettel committed Apr 10, 2018 240 241  pfok <- 1;if (kp==1) { pf <- paste(pf,"1");  Dirk Eddelbuettel committed Apr 10, 2018 242 243  } }  Dirk Eddelbuettel committed Apr 10, 2018 244   Dirk Eddelbuettel committed Apr 10, 2018 245  fake.formula <- pf  Dirk Eddelbuettel committed Apr 10, 2018 246   Dirk Eddelbuettel committed Apr 10, 2018 247  if (length(smooth.spec)>0)  Dirk Eddelbuettel committed Apr 10, 2018 248 249  for (i in 1:length(smooth.spec)) { nt <- length(smooth.spec[[i]]$term)  Dirk Eddelbuettel committed Apr 10, 2018 250 251  ff1 <- paste(smooth.spec[[i]]$term[1:nt],collapse="+") fake.formula <- paste(fake.formula,"+",ff1)  Dirk Eddelbuettel committed Apr 10, 2018 252 253 254 255  if (smooth.spec[[i]]$by!="NA") { fake.formula <- paste(fake.formula,"+",smooth.spec[[i]]$by) av <- c(av,smooth.spec[[i]]$term,smooth.spec[[i]]$by) } else av <- c(av,smooth.spec[[i]]$term)  Dirk Eddelbuettel committed Apr 10, 2018 256  }  Dirk Eddelbuettel committed Apr 10, 2018 257  fake.formula <- as.formula(fake.formula,p.env)  Dirk Eddelbuettel committed Apr 10, 2018 258 259 260 261 262  if (length(av)) { pred.formula <- as.formula(paste("~",paste(av,collapse="+"))) pav <- all.vars(pred.formula) ## trick to strip out 'offset(x)' etc... pred.formula <- reformulate(pav) } else pred.formula <- ~1  Dirk Eddelbuettel committed Apr 10, 2018 263  ret <- list(pf=as.formula(pf,p.env),pfok=pfok,smooth.spec=smooth.spec,  Dirk Eddelbuettel committed Apr 10, 2018 264 265  fake.formula=fake.formula,response=response,fake.names=av, pred.names=pav,pred.formula=pred.formula)  Dirk Eddelbuettel committed Apr 10, 2018 266  class(ret) <- "split.gam.formula"  Dirk Eddelbuettel committed Apr 10, 2018 267  ret  Dirk Eddelbuettel committed Apr 10, 2018 268 269 270 271 272 273 274 } ## interpret.gam0 interpret.gam <- function(gf) { ## wrapper to allow gf to be a list of formulae or ## a single formula. This facilitates general penalized ## likelihood models in which several linear predictors ## may be involved...  Dirk Eddelbuettel committed Apr 10, 2018 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 ## ## The list syntax is as follows. The first formula must have a response on ## the lhs, rather than labels. For m linear predictors, there ## must be m 'base formulae' in linear predictor order. lhs labels will ## be ignored in a base formula. Empty base formulae have '-1' on rhs. ## Further formulae have labels up to m labels 1,...,m on the lhs, in a ## syntax like this: 3 + 5 ~ s(x), which indicates that the same s(x) ## should be added to both linear predictors 3 and 5. ## e.g. A bivariate normal model with common expected values might be ## list(y1~-1,y2~-1,1+2~s(x)), whereas if the second component was contaminated ## by something else we might have list(y1~-1,y2~s(v)-1,1+2~s(x)) ## ## For a list argument, this routine returns a list of slit.formula objects ## with an extra field "lpi" indicating the linear predictors to which each ## contributes...  Dirk Eddelbuettel committed Apr 10, 2018 290 291  if (is.list(gf)) { d <- length(gf)  Dirk Eddelbuettel committed Apr 10, 2018 292 293 294 295 296 297  ## make sure all formulae have a response, to avoid ## problems with parametric sub formulae of the form ~1 #if (length(gf[[1]])<3) stop("first formula must specify a response") resp <- gf[[1]][2]  Dirk Eddelbuettel committed Apr 10, 2018 298  ret <- list()  Dirk Eddelbuettel committed Apr 10, 2018 299  pav <- av <- rep("",0)  Dirk Eddelbuettel committed Apr 10, 2018 300  nlp <- 0 ## count number of linear predictors (may be different from number of formulae)  Dirk Eddelbuettel committed Apr 10, 2018 301 302  for (i in 1:d) { textra <- if (i==1) NULL else paste(".",i-1,sep="") ## modify smooth labels to identify to predictor  Dirk Eddelbuettel committed Apr 10, 2018 303 304 305 306 307 308  lpi <- getNumericResponse(gf[[i]]) ## get linear predictors to which this applies, if explicit if (length(lpi)==1) warning("single linear predictor indices are ignored") if (length(lpi)>0) gf[[i]][[2]] <- NULL else { ## delete l.p. labels from formula response nlp <- nlp + 1;lpi <- nlp ## this is base formula for l.p. number nlp }  Dirk Eddelbuettel committed Apr 10, 2018 309  ret[[i]] <- interpret.gam0(gf[[i]],textra)  Dirk Eddelbuettel committed Apr 10, 2018 310 311  ret[[i]]$lpi <- lpi ## record of the linear predictors to which this applies  Dirk Eddelbuettel committed Apr 10, 2018 312 313 314 315 316 317 318 319 320  ## make sure all parametric formulae have a response, to avoid ## problems with parametric sub formulae of the form ~1 respi <- rep("",0) ## no extra response terms if (length(ret[[i]]$pf)==2) { ret[[i]]$pf[3] <- ret[[i]]$pf[2];ret[[i]]$pf[2] <- resp respi <- rep("",0) } else if (i>1) respi <- ret[[i]]$response ## extra response terms av <- c(av,ret[[i]]$fake.names,respi) ## accumulate all required variable names pav <- c(pav,ret[[i]]$pred.names) ## predictors only  Dirk Eddelbuettel committed Apr 10, 2018 321 322  } av <- unique(av) ## strip out duplicate variable names  Dirk Eddelbuettel committed Apr 10, 2018 323  pav <- unique(pav)  Dirk Eddelbuettel committed Apr 10, 2018 324 325  ret$fake.formula <- if (length(av)>0) reformulate(av,response=ret[[1]]$response) else ret[[1]]$fake.formula ## create fake formula containing all variables  Dirk Eddelbuettel committed Apr 10, 2018 326  ret$pred.formula <- if (length(pav)>0) reformulate(pav) else ~1 ## predictor only formula  Dirk Eddelbuettel committed Apr 10, 2018 327  ret$response <- ret[[1]]$response  Dirk Eddelbuettel committed Apr 10, 2018 328 329  ret$nlp <- nlp ## number of linear predictors for (i in 1:d) if (max(ret[[i]]$lpi)>nlp||min(ret[[i]]$lpi)<1) stop("linear predictor labels out of range")  Dirk Eddelbuettel committed Apr 10, 2018 330 331 332 333  class(ret) <- "split.gam.formula" return(ret) } else interpret.gam0(gf) } ## interpret.gam  Dirk Eddelbuettel committed Apr 10, 2018 334   Dirk Eddelbuettel committed Apr 10, 2018 335   Dirk Eddelbuettel committed Apr 10, 2018 336 fixDependence <- function(X1,X2,tol=.Machine$double.eps^.5,rank.def=0,strict=FALSE)  Dirk Eddelbuettel committed Apr 10, 2018 337 338 # model matrix X2 may be linearly dependent on X1. This # routine finds which columns of X2 should be zeroed to  Dirk Eddelbuettel committed Apr 10, 2018 339 340 # fix this. If rank.def>0 then it is taken as the known degree # of dependence of X2 on X1 and tol is ignored.  Dirk Eddelbuettel committed Apr 10, 2018 341 { qr1 <- qr(X1,LAPACK=TRUE)  Dirk Eddelbuettel committed Apr 10, 2018 342  R11 <- abs(qr.R(qr1)[1,1])  Dirk Eddelbuettel committed Apr 10, 2018 343  r<-ncol(X1);n<-nrow(X1)  Dirk Eddelbuettel committed Apr 10, 2018 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360  if (strict) { ## only delete columns of X2 individually dependent on X1 ## Project columns of X2 into space of X1 and look at difference ## to orignal X2 to check for deficiency... QtX2 <- qr.qty(qr1,X2) QtX2[-(1:r),] <- 0 mdiff <- colMeans(abs(X2 - qr.qy(qr1,QtX2))) if (rank.def>0) ind <- (1:ncol(X2))[rank(mdiff) <= rank.def] else ind <- (1:ncol(X2))[mdiff < R11*tol] if (length(ind)<1) ind <- NULL } else { ## make X2 full rank given X1 QtX2 <- qr.qty(qr1,X2)[(r+1):n,] # Q'X2 qr2 <- qr(QtX2,LAPACK=TRUE) R <- qr.R(qr2) # now final diagonal block of R may be zero, indicating rank # deficiency. r0 <- r <- nrow(R) if (rank.def > 0 && rank.def <= nrow(R)) r0 <- r - rank.def else ## degree of rank def known  Dirk Eddelbuettel committed Apr 10, 2018 361  while (r0>0 && mean(abs(R[r0:r,r0:r]))< R11*tol) r0 <- r0 -1 ## compute rank def  Dirk Eddelbuettel committed Apr 10, 2018 362 363 364 365 366  r0 <- r0 + 1 if (r0>r) return(NULL) else ind <- qr2$pivot[r0:r] # the columns of X2 to zero in order to get independence } ind  Dirk Eddelbuettel committed Apr 10, 2018 367 } ## fixDependence  Dirk Eddelbuettel committed Apr 10, 2018 368 369   Dirk Eddelbuettel committed Apr 10, 2018 370 371 372 373 374 augment.smX <- function(sm,nobs,np) { ## augments a smooth model matrix with a square root penalty matrix for ## identifiability constraint purposes. ns <- length(sm$S) ## number of penalty matrices if (ns==0) { ## nothing to do  Dirk Eddelbuettel committed Apr 10, 2018 375  return(rbind(sm$X,matrix(0,np,ncol(sm$X))))  Dirk Eddelbuettel committed Apr 10, 2018 376  }  Dirk Eddelbuettel committed Apr 10, 2018 377 378 379 380 381 382 383 384 385  ind <- colMeans(abs(sm$S[[1]]))!=0 sqrmaX <- mean(abs(sm$X[,ind]))^2 alpha <- sqrmaX/mean(abs(sm$S[[1]][ind,ind])) St <- sm$S[[1]]*alpha if (ns>1) for (i in 2:ns) { ind <- colMeans(abs(sm$S[[i]]))!=0 alpha <- sqrmaX/mean(abs(sm$S[[i]][ind,ind])) St <- St + sm$S[[i]]*alpha }  Dirk Eddelbuettel committed Apr 10, 2018 386  rS <- mroot(St,rank=ncol(St)) ## get sqrt of penalty  Dirk Eddelbuettel committed Apr 10, 2018 387  X <- rbind(sm$X,matrix(0,np,ncol(sm$X))) ## create augmented model matrix  Dirk Eddelbuettel committed Apr 10, 2018 388 389  X[nobs+sm$p.ind,] <- t(rS) ## add in X ## scaled augmented model matrix  Dirk Eddelbuettel committed Apr 10, 2018 390 } ## augment.smX  Dirk Eddelbuettel committed Apr 10, 2018 391 392  gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)  Dirk Eddelbuettel committed Apr 10, 2018 393 394 395 # works through a list of smooths, sm, aiming to identify nested or partially # nested terms, and impose identifiability constraints on them. # Xp is the parametric model matrix. It is needed in order to check whether  Dirk Eddelbuettel committed Apr 10, 2018 396 # there is a constant (or equivalent) in the model. If there is, then this needs  Dirk Eddelbuettel committed Apr 10, 2018 397 398 # to be included when working out side constraints, otherwise dependencies can be # missed.  Dirk Eddelbuettel committed Apr 10, 2018 399 400 # Note that with.pen is quite extreme, since you then pretty much only pick # up dependencies in the null spaces  Dirk Eddelbuettel committed Apr 10, 2018 401 402 403 404 { if (!with.pen) { ## check that's possible and reset if not! with.pen <- nrow(Xp) < ncol(Xp) + sum(unlist(lapply(sm,function(x) ncol(x$X)))) } m <- length(sm)  Dirk Eddelbuettel committed Apr 10, 2018 405 406 407 408 409 410  if (m==0) return(sm) v.names<-array("",0);maxDim<-1 for (i in 1:m) { ## collect all term names and max smooth dim' vn <- sm[[i]]$term ## need to include by variables in names if (sm[[i]]$by!="NA") vn <- paste(vn,sm[[i]]$by,sep="")  Dirk Eddelbuettel committed Apr 10, 2018 411 412  ## need to distinguish levels of factor by variables... if (!is.null(sm[[i]]$by.level)) vn <- paste(vn,sm[[i]]$by.level,sep="")  Dirk Eddelbuettel committed Apr 10, 2018 413 414 415 416 417 418 419  sm[[i]]$vn <- vn ## use this record to identify variables from now v.names <- c(v.names,vn) if (sm[[i]]$dim > maxDim) maxDim <- sm[[i]]$dim } lv <- length(v.names) v.names <- unique(v.names) if (lv == length(v.names)) return(sm) ## no repeats => no nesting  Dirk Eddelbuettel committed Apr 10, 2018 420 421 422 423 424  ## Only get this far if there is nesting. ## Need to test for intercept or equivalent in Xp intercept <- FALSE if (ncol(Xp)) {  Dirk Eddelbuettel committed Apr 10, 2018 425  ## first check columns directly...  Dirk Eddelbuettel committed Apr 10, 2018 426  if (sum(apply(Xp,2,sd)<.Machine$double.eps^.75)>0) intercept <- TRUE else {  Dirk Eddelbuettel committed Apr 10, 2018 427  ## no constant column, so need to check span of Xp...  Dirk Eddelbuettel committed Apr 10, 2018 428 429 430 431  f <- rep(1,nrow(Xp)) ff <- qr.fitted(qr(Xp),f) if (max(abs(ff-f))<.Machine$double.eps^.75) intercept <- TRUE }  Dirk Eddelbuettel committed Apr 10, 2018 432 433  }  Dirk Eddelbuettel committed Apr 10, 2018 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451  sm.id <- as.list(v.names) names(sm.id) <- v.names for (i in 1:length(sm.id)) sm.id[[i]]<-array(0,0) sm.dim <- sm.id for (d in 1:maxDim) { for (i in 1:m) { if (sm[[i]]$dim==d) for (j in 1:d) { ## work through terms term<-sm[[i]]$vn[j] a <- sm.id[[term]] la <- length(a)+1 sm.id[[term]][la] <- i ## record smooth i.d. for this variable sm.dim[[term]][la] <- d ## ... and smooth dim. } } } ## so now each unique variable name has an associated array of ## the smooths of which it is an argument, arranged in ascending ## order of dimension.  Dirk Eddelbuettel committed Apr 10, 2018 452  if (maxDim==1) warning("model has repeated 1-d smooths of same variable.")  Dirk Eddelbuettel committed Apr 10, 2018 453 454 455 456 457 458 459 460 461 462 463 464 465 466  ## Now set things up to enable term specific model matrices to be ## augmented with square root penalties, on the fly... if (with.pen) { k <- 1 for (i in 1:m) { ## create parameter indices for each term k1 <- k + ncol(sm[[i]]$X) - 1 sm[[i]]$p.ind <- k:k1 k <- k1 + 1 } np <- k-1 ## number of penalized parameters } nobs <- nrow(sm[[1]]$X) ## number of observations  Dirk Eddelbuettel committed Apr 10, 2018 467  for (d in 1:maxDim) { ## work up through dimensions  Dirk Eddelbuettel committed Apr 10, 2018 468  for (i in 1:m) { ## work through smooths  Dirk Eddelbuettel committed Apr 10, 2018 469 470 471  if (sm[[i]]$dim == d&&sm[[i]]$side.constrain) { ## check for nesting if (with.pen) X1 <- matrix(c(rep(1,nobs),rep(0,np)),nobs+np,as.integer(intercept)) else X1 <- matrix(1,nobs,as.integer(intercept))  Dirk Eddelbuettel committed Apr 10, 2018 472 473 474 475  for (j in 1:d) { ## work through variables b <- sm.id[[sm[[i]]$vn[j]]] # list of smooths dependent on this variable k <- (1:length(b))[b==i] ## locate current smooth in list if (k>1) for (l in 1:(k-1)) { ## collect X columns  Dirk Eddelbuettel committed Apr 10, 2018 476 477 478 479  if (with.pen) { ## need to use augmented model matrix in testing if (is.null(sm[[b[l]]]$Xa)) sm[[b[l]]]$Xa <- augment.smX(sm[[b[l]]],nobs,np) X1 <- cbind(X1,sm[[b[l]]]$Xa) } else X1 <- cbind(X1,sm[[b[l]]]$X) ## penalties not considered  Dirk Eddelbuettel committed Apr 10, 2018 480 481  } } ## Now X1 contains columns for all lower dimensional terms  Dirk Eddelbuettel committed Apr 10, 2018 482 483 484 485 486 487  if (ncol(X1)==as.integer(intercept)) ind <- NULL else { if (with.pen) { if (is.null(sm[[i]]$Xa)) sm[[i]]$Xa <- augment.smX(sm[[i]],nobs,np) ind <- fixDependence(X1,sm[[i]]$Xa,tol=tol) } else ind <- fixDependence(X1,sm[[i]]$X,tol=tol) }  Dirk Eddelbuettel committed Apr 10, 2018 488 489 490  ## ... the columns to zero to ensure independence if (!is.null(ind)) { sm[[i]]$X <- sm[[i]]$X[,-ind]  Dirk Eddelbuettel committed Apr 10, 2018 491  ## work through list of penalty matrices, applying constraints...  Dirk Eddelbuettel committed Apr 10, 2018 492 493  nsmS <- length(sm[[i]]$S) if (nsmS>0) for (j in nsmS:1) { ## working down so that dropping is painless  Dirk Eddelbuettel committed Apr 10, 2018 494  sm[[i]]$S[[j]] <- sm[[i]]$S[[j]][-ind,-ind]  Dirk Eddelbuettel committed Apr 10, 2018 495 496  if (sum(sm[[i]]$S[[j]]!=0)==0) rank <- 0 else rank <- qr(sm[[i]]$S[[j]],tol=tol,LAPACK=FALSE)$rank  Dirk Eddelbuettel committed Apr 10, 2018 497  sm[[i]]$rank[j] <- rank ## replace previous rank with new rank  Dirk Eddelbuettel committed Apr 10, 2018 498 499 500  if (rank == 0) { ## drop the penalty sm[[i]]$rank <- sm[[i]]$rank[-j] sm[[i]]$S[[j]] <- NULL  Dirk Eddelbuettel committed Apr 10, 2018 501  sm[[i]]$S.scale <- sm[[i]]$S.scale[-j]  Dirk Eddelbuettel committed Apr 10, 2018 502 503 504  if (!is.null(sm[[i]]$L)) sm[[i]]$L <- sm[[i]]$L[-j,,drop=FALSE] } } ## penalty matrices finished  Dirk Eddelbuettel committed Apr 10, 2018 505  ## Now we need to establish null space rank for the term  Dirk Eddelbuettel committed Apr 10, 2018 506 507  mi <- length(sm[[i]]$S) if (mi>0) {  Dirk Eddelbuettel committed Apr 10, 2018 508  St <- sm[[i]]$S[[1]]/norm(sm[[i]]$S[[1]],type="F")  Dirk Eddelbuettel committed Apr 10, 2018 509  if (mi>1) for (j in 1:mi) St <- St +  Dirk Eddelbuettel committed Apr 10, 2018 510 511 512 513 514  sm[[i]]$S[[j]]/norm(sm[[i]]$S[[j]],type="F") es <- eigen(St,symmetric=TRUE,only.values=TRUE) sm[[i]]$null.space.dim <- sum(es$values don't clone } specb ## return clone  Dirk Eddelbuettel committed Apr 10, 2018 582 } ## clone.smooth.spec  Dirk Eddelbuettel committed Apr 10, 2018 583 584   Dirk Eddelbuettel committed Apr 10, 2018 585 586 587 588 589 590 591 592 593 parametricPenalty <- function(pterms,assign,paraPen,sp0) { ## routine to process any penalties on the parametric part of the model. ## paraPen is a list whose items have names corresponding to the ## term.labels in pterms. Each list item may have named elements ## L, rank and sp. All other elements should be penalty coefficient matrices. S <- list() ## penalty matrix list off <- rep(0,0) ## offset array rank <- rep(0,0) ## rank array sp <- rep(0,0) ## smoothing param array  Dirk Eddelbuettel committed Apr 10, 2018 594  full.sp.names <- rep("",0) ## names for sp's multiplying penalties (not underlying)  Dirk Eddelbuettel committed Apr 10, 2018 595 596 597 598  L <- matrix(0,0,0) k <- 0 tind <- unique(assign) ## unique term indices n.t <- length(tind)  Dirk Eddelbuettel committed Apr 10, 2018 599 600  if (n.t>0) for (j in 1:n.t) if (tind[j]>0) { term.label <- attr(pterms[tind[j]],"term.label")  Dirk Eddelbuettel committed Apr 10, 2018 601 602  P <- paraPen[[term.label]] ## get any penalty information for this term if (!is.null(P)) { ## then there is information  Dirk Eddelbuettel committed Apr 10, 2018 603  ind <- (1:length(assign))[assign==tind[j]] ## index of coefs involved here  Dirk Eddelbuettel committed Apr 10, 2018 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626  Li <- P$L;P$L <- NULL spi <- P$sp;P$sp <- NULL ranki <- P$rank;P$rank <- NULL ## remaining terms should be penalty matrices... np <- length(P) if (!is.null(ranki)&&length(ranki)!=np) stop("rank' has wrong length in paraPen'") if (np) for (i in 1:np) { ## unpack penalty matrices, offsets and ranks k <- k + 1 S[[k]] <- P[[i]] off[k] <- min(ind) ## index of first coef penalized by this term if ( ncol(P[[i]])!=nrow(P[[i]])||nrow(P[[i]])!=length(ind)) stop(" a parametric penalty has wrong dimension") if (is.null(ranki)) { ev <- eigen(S[[k]],symmetric=TRUE,only.values=TRUE)$values rank[k] <- sum(ev>max(ev)*.Machine$double.eps*10) ## estimate rank } else rank[k] <- ranki[i] } ## now deal with L matrices if (np) { ## only do this stuff if there are any penalties! if (is.null(Li)) Li <- diag(np) if (nrow(Li)!=np) stop("L has wrong dimension in paraPen'") L <- rbind(cbind(L,matrix(0,nrow(L),ncol(Li))), cbind(matrix(0,nrow(Li),ncol(L)),Li))  Dirk Eddelbuettel committed Apr 10, 2018 627  ind <- (length(sp)+1):(length(sp)+ncol(Li))  Dirk Eddelbuettel committed Apr 10, 2018 628  ind2 <- (length(sp)+1):(length(sp)+nrow(Li)) ## used to produce names for full sp array  Dirk Eddelbuettel committed Apr 10, 2018 629  if (is.null(spi)) {  Dirk Eddelbuettel committed Apr 10, 2018 630  sp[ind] <- -1 ## auto-initialize  Dirk Eddelbuettel committed Apr 10, 2018 631 632  } else { if (length(spi)!=ncol(Li)) stop("sp' dimension wrong in paraPen'")  Dirk Eddelbuettel committed Apr 10, 2018 633  sp[ind] <- spi  Dirk Eddelbuettel committed Apr 10, 2018 634  }  Dirk Eddelbuettel committed Apr 10, 2018 635 636 637  ## add smoothing parameter names.... if (length(ind)>1) names(sp)[ind] <- paste(term.label,ind-ind[1]+1,sep="") else names(sp)[ind] <- term.label  Dirk Eddelbuettel committed Apr 10, 2018 638 639 640  if (length(ind2)>1) full.sp.names[ind2] <- paste(term.label,ind2-ind2[1]+1,sep="") else full.sp.names[ind2] <- term.label  Dirk Eddelbuettel committed Apr 10, 2018 641 642 643 644 645 646 647 648 649 650 651 652  } } ## end !is.null(P) } ## looped through all terms if (k==0) return(NULL) if (!is.null(sp0)) { if (length(sp0)0) return(rep(0,0)) ## deparse turns lhs into a string; strsplit extracts the characters ## corresponding to numbers; unlist deals with the fact that deparse ## will split long lines resulting in multiple list items from ## strsplit; as.numeric converts the numbers; na.omit drops NAs ## resulting from "" elements; unique & round are obvious... round(unique(na.omit(as.numeric(unlist(strsplit(deparse(form[[2]]), "[^0-9]+")))))) } ## getNumericResponse olid <- function(X,nsdf,pstart,flpi,lpi) { ## X is a model matrix, made up of nf=length(nsdf) column blocks. ## The ith block starts at column pstart[i] and its first nsdf[i] ## columns are unpenalized. X is used to define nlp=length(lpi) ## linear predictors. lpi[[i]] gives the columns of X used in the ## ith linear predictor. flpi[j] gives the linear predictor(s) ## to which the jth block of X belongs. The problem is that the ## unpenalized blocks need not be identifiable when used in combination. ## This function returns a vector dind of columns of X to drop for ## identifiability, along with modified lpi, pstart and nsdf vectors. nlp <- length(lpi) ## number of linear predictors n <- nrow(X) nf <- length(nsdf) ## number of formulae blocks Xp <- matrix(0,n*nlp,sum(nsdf)) start <- 1 ii <- 1:n tind <- rep(0,0) ## complete index of all parametric columns in X ## create a block matrix, Xp, with the same identifiability properties as ## unpenalized part of model... for (i in 1:nf) { stop <- start - 1 + nsdf[i] if (stop>=start) { ind <- pstart[i] + 1:nsdf[i] - 1 for (k in flpi[[i]]) { Xp[ii+(k-1)*n,start:stop] <- X[,ind] } tind <- c(tind,ind) start <- start + nsdf[i] } } ## rank deficiency of Xp will reveal number of redundant parametric ## terms, and a pivoted QR will reveal which to drop to restore ## full rank... qrx <- qr(Xp,LAPACK=TRUE,tol=0.0) ## unidentifiable columns get pivoted to final cols r <- Rrank(qr.R(qrx)) ## get rank from R factor of pivoted QR if (r==ncol(Xp)) { ## full rank, all fine, drop nothing dind <- rep(0,0) } else { ## reduced rank, drop some columns dind <- tind[sort(qrx$pivot[(r+1):ncol(X)],decreasing=TRUE)] ## columns to drop ## now we need to adjust nsdf, pstart and lpi for (d in dind) { ## working down through drop indices ## following commented out code is useful should it ever prove necessary to ## adjust pstart and nsdf, but at present these are only used in prediction, ## and it is cleaner to leave them unchanged, and simply drop using dind during prediction. #k <- if (d>=pstart[nf]) nlp else which(d >= pstart[1:(nf-1)] & d < pstart[2:nf]) #nsdf[k] <- nsdf[k] - 1 ## one less unpenalized column in this block #if (k0) lpi[[i]] <- lpi[[i]][-k] ## drop row k <- which(lpi[[i]]>d) if (length(k)>0) lpi[[i]][k] <- lpi[[i]][k] - 1 ## close up } } ## end of drop index loop } list(dind=dind,lpi=lpi) ##,pstart=pstart,nsdf=nsdf) } ## olid  Dirk Eddelbuettel committed Apr 10, 2018 741 gam.setup.list <- function(formula,pterms,  Dirk Eddelbuettel committed Apr 10, 2018 742  data=stop("No data supplied to gam.setup"),knots=NULL,sp=NULL,  Dirk Eddelbuettel committed Apr 10, 2018 743 744  min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,select=FALSE,idLinksBases=TRUE, scale.penalty=TRUE,paraPen=NULL,gamm.call=FALSE,drop.intercept=FALSE) {  Dirk Eddelbuettel committed Apr 10, 2018 745 746 747 748 ## version of gam.setup for when gam is called with a list of formulae, ## specifying several linear predictors... ## key difference to gam.setup is an attribute to the model matrix, "lpi", which is a list ## of column indices for each linear predictor  Dirk Eddelbuettel committed Apr 10, 2018 749 750  if (!is.null(paraPen)) stop("paraPen not supported for multi-formula models") if (!absorb.cons) stop("absorb.cons must be TRUE for multi-formula models")  Dirk Eddelbuettel committed Apr 10, 2018 751 752 753 754  d <- length(pterms) ## number of formulae lp.overlap <- if (formula$nlp0) um$term.names[1:um$nsdf] <- paste(um$term.names[1:um$nsdf],i-1,sep=".") G$term.names <- c(G$term.names,um$term.names) G$lsp0 <- c(G$lsp0,um$lsp0)  Dirk Eddelbuettel committed Apr 10, 2018 835 836  G$sp <- c(G$sp,um$sp) pof <- ncol(G$X)  Dirk Eddelbuettel committed Apr 10, 2018 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861  } ## formula loop end ## If there is overlap then there is a danger of lack of identifiability of the ## parameteric terms, especially if there are factors present in shared components. ## The following code deals with this possibility... if (lp.overlap) { rt <- olid(G$X,G$nsdf,pstart,flpi,lpi) if (length(rt$dind)>0) { ## then columns have to be dropped warning("dropping unidentifiable parameteric terms from model",call.=FALSE) G$X <- G$X[,-rt$dind] ## drop cols G$cmX <- G$cmX[-rt$dind] G$term.names <- G$term.names[-rt$dind] ## adjust indexing in smooth list, noting that coefs of smooths ## are never dropped by dind for (i in 1:length(G$smooth)) { k <- sum(rt$dind < G$smooth[[i]]$first.para) G$smooth[[i]]$first.para <- G$smooth[[i]]$first.para - k G$smooth[[i]]$last.para <- G$smooth[[i]]$last.para - k } for (i in 1:length(G$off)) G$off[i] <- G$off[i] - sum(rt$dind < G$off[i]) ## replace various indices with updated versions... # pstart <- rt$pstart; G$nsdf <- rt$nsdf ## these two only needed in predict.gam - cleaner to leave unchanged lpi <- rt$lpi attr(G$nsdf,"drop.ind") <- rt$dind ## store drop index }  Dirk Eddelbuettel committed Apr 10, 2018 862  }  Dirk Eddelbuettel committed Apr 10, 2018 863  attr(lpi,"overlap") <- lp.overlap  Dirk Eddelbuettel committed Apr 10, 2018 864  attr(G$X,"lpi") <- lpi  Dirk Eddelbuettel committed Apr 10, 2018 865  attr(G$nsdf,"pstart") <- pstart ##unlist(lapply(lpi,min))  Dirk Eddelbuettel committed Apr 10, 2018 866 867 868  G } ## gam.setup.list  Dirk Eddelbuettel committed Apr 10, 2018 869 870   Dirk Eddelbuettel committed Apr 10, 2018 871 872 gam.setup <- function(formula,pterms, data=stop("No data supplied to gam.setup"),knots=NULL,sp=NULL,  Dirk Eddelbuettel committed Apr 10, 2018 873  min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,select=FALSE,idLinksBases=TRUE,  Dirk Eddelbuettel committed Apr 10, 2018 874  scale.penalty=TRUE,paraPen=NULL,gamm.call=FALSE,drop.intercept=FALSE,  Dirk Eddelbuettel committed Apr 10, 2018 875  diagonal.penalty=FALSE,apply.by=TRUE)  Dirk Eddelbuettel committed Apr 10, 2018 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 ## set up the model matrix, penalty matrices and auxilliary information about the smoothing bases ## needed for a gam fit. ## elements of returned object: ## * m - number of smooths ## * min.sp - minimum smoothing parameters ## * H supplied H matrix ## * pearson.extra, dev.extra, n.true --- entries to hold these quantities ## * pterms - terms object for parametric terms ## * intercept TRUE if intercept present ## * offset - the model offset ## * nsdf - number of strictly parameteric coefs ## * contrasts ## * xlevels - records levels of factors ## * assign - indexes which parametric model matrix columns map to which term in pterms ## * smooth - list of smooths ## * S - penalties (non-zero block only) ## * off - first coef penalized by each element of S ## * cmX - col mean of X ## * P - maps parameters in fit constraint parameterization to those in prediction parameterization ## * X - model matrix ## * sp ## * rank ## * n.paraPen ## * L ## * lsp0 ## * y - response ## * C - constraint matrix - only if absorb.cons==FALSE ## * n - dim(y) ## * w - weights ## * term.names ## * nP  Dirk Eddelbuettel committed Apr 10, 2018 907 { # split the formula if the object being passed is a formula, otherwise it's already split  Dirk Eddelbuettel committed Apr 10, 2018 908   Dirk Eddelbuettel committed Apr 10, 2018 909 910  if (inherits(formula,"split.gam.formula")) split <- formula else if (inherits(formula,"formula")) split <- interpret.gam(formula)  Dirk Eddelbuettel committed Apr 10, 2018 911 912  else stop("First argument is no sort of formula!")  Dirk Eddelbuettel committed Apr 10, 2018 913 914 915 916 917 918 919  if (length(split$smooth.spec)==0) { if (split$pfok==0) stop("You've got no model....") m <- 0 } else m <- length(split$smooth.spec) # number of smooth terms G <- list(m=m,min.sp=min.sp,H=H,pearson.extra=0, dev.extra=0,n.true=-1,pterms=pterms) ## dev.extra gets added to deviance if REML/ML used in gam.fit3  Dirk Eddelbuettel committed Apr 10, 2018 920 921  if (is.null(attr(data,"terms"))) # then data is not a model frame  Dirk Eddelbuettel committed Apr 10, 2018 922 923  mf <- model.frame(split$pf,data,drop.unused.levels=FALSE) # must be false or can end up with wrong prediction matrix! else mf <- data # data is already a model frame  Dirk Eddelbuettel committed Apr 10, 2018 924 925 926  G$intercept <- attr(attr(mf,"terms"),"intercept")>0 G$offset <- model.offset(mf) # get the model offset (if any)  Dirk Eddelbuettel committed Apr 10, 2018 927  if (!is.null(G$offset)) G$offset <- as.numeric(G$offset)  Dirk Eddelbuettel committed Apr 10, 2018 928   Dirk Eddelbuettel committed Apr 10, 2018 929  # construct strictly parametric model matrix....  Dirk Eddelbuettel committed Apr 10, 2018 930  if (drop.intercept) attr(pterms,"intercept") <- 1 ## ensure there is an intercept to drop  Dirk Eddelbuettel committed Apr 10, 2018 931  X <- model.matrix(pterms,mf)  Dirk Eddelbuettel committed Apr 10, 2018 932 933 934 935 936 937 938  if (drop.intercept) { ## some extended families require intercept to be dropped xat <- attributes(X);ind <- xat$assign>0 X <- X[,xat$assign>0,drop=FALSE] ## some extended families need to drop intercept xat$assign <- xat$assign[ind];xat$dimnames[[2]]<-xat$dimnames[[2]][ind]; xat$dim[2] <- xat$dim[2]-1;attributes(X) <- xat G$intercept <- FALSE }  Dirk Eddelbuettel committed Apr 10, 2018 939  rownames(X) <- NULL ## save memory  Dirk Eddelbuettel committed Apr 10, 2018 940   Dirk Eddelbuettel committed Apr 10, 2018 941 942 943 944 945  G$nsdf <- ncol(X) G$contrasts <- attr(X,"contrasts") G$xlevels <- .getXlevels(pterms,mf) G$assign <- attr(X,"assign") # used to tell which coeffs relate to which pterms  Dirk Eddelbuettel committed Apr 10, 2018 946 947  ## now deal with any user supplied penalties on the parametric part of the model... PP <- parametricPenalty(pterms,G$assign,paraPen,sp)  Dirk Eddelbuettel committed Apr 10, 2018 948 949 950 951 952 953 954 955  if (!is.null(PP)) { ## strip out supplied sps already used ind <- 1:length(PP$sp) if (!is.null(sp)) sp <- sp[-ind] if (!is.null(min.sp)) { PP$min.sp <- min.sp[ind] min.sp <- min.sp[-ind] } }  Dirk Eddelbuettel committed Apr 10, 2018 956 957 958  # next work through smooth terms (if any) extending model matrix.....  Dirk Eddelbuettel committed Apr 10, 2018 959 960  G$smooth <- list() G$S <- list()  Dirk Eddelbuettel committed Apr 10, 2018 961   Dirk Eddelbuettel committed Apr 10, 2018 962  if (gamm.call) { ## flag that this is a call from gamm --- some smoothers need to know!  Dirk Eddelbuettel committed Apr 10, 2018 963 964  if (m>0) for (i in 1:m) attr(split$smooth.spec[[i]],"gamm") <- TRUE }  Dirk Eddelbuettel committed Apr 10, 2018 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980  if (m>0 && idLinksBases) { ## search smooth.spec[[]] for terms linked by common id's id.list <- list() ## id information list for (i in 1:m) if (!is.null(split$smooth.spec[[i]]$id)) { id <- as.character(split$smooth.spec[[i]]$id) if (length(id.list)&&id%in%names(id.list)) { ## it's an existing id ni <- length(id.list[[id]]$sm.i) ## number of terms so far with this id id.list[[id]]$sm.i[ni+1] <- i ## adding smooth.spec index to this id's list ## clone smooth.spec from base smooth spec.... base.i <- id.list[[id]]$sm.i[1] split$smooth.spec[[i]] <- clone.smooth.spec(split$smooth.spec[[base.i]], split$smooth.spec[[i]]) ## add data for this term to the data list for basis setup... temp.term <- split$smooth.spec[[i]]$term  Dirk Eddelbuettel committed Apr 10, 2018 981 982 983  ## note cbind deliberate in next line, as construction will handle matrix argument ## correctly...  Dirk Eddelbuettel committed Apr 10, 2018 984 985  for (j in 1:length(temp.term)) id.list[[id]]$data[[j]] <- cbind(id.list[[id]]$data[[j]], get.var(temp.term[j],data,vecMat=FALSE))  Dirk Eddelbuettel committed Apr 10, 2018 986 987  } else { ## new id  Dirk Eddelbuettel committed Apr 10, 2018 988 989 990 991 992 993 994 995 996 997  id.list[[id]] <- list(sm.i=i) ## start the array of indices of smooths with this id id.list[[id]]$data <- list() ## need to collect together all data for which this basis will be used, ## for basis setup... term <- split$smooth.spec[[i]]$term for (j in 1:length(term)) id.list[[id]]$data[[j]] <- get.var(term[j],data,vecMat=FALSE) } ## new id finished } } ## id.list complete  Dirk Eddelbuettel committed Apr 10, 2018 998 999  G$off<-array(0,0) first.para<-G$nsdf+1  Dirk Eddelbuettel committed Apr 10, 2018 1000  sm <- list()  Dirk Eddelbuettel committed Apr 10, 2018 1001  newm <- 0  Dirk Eddelbuettel committed Apr 10, 2018 1002 1003  if (m>0) for (i in 1:m) { # idea here is that terms are set up in accordance with information given in split$smooth.spec  Dirk Eddelbuettel committed Apr 10, 2018 1004 1005  # appropriate basis constructor is called depending on the class of the smooth # constructor returns penalty matrices model matrix and basis specific information  Dirk Eddelbuettel committed Apr 10, 2018 1006  ## sm[[i]] <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,sparse.cons=sparse.cons) ## old code  Dirk Eddelbuettel committed Apr 10, 2018 1007 1008  id <- split$smooth.spec[[i]]$id if (is.null(id)||!idLinksBases) { ## regular evaluation  Dirk Eddelbuettel committed Apr 10, 2018 1009  sml <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,  Dirk Eddelbuettel committed Apr 10, 2018 1010 1011  null.space.penalty=select,sparse.cons=sparse.cons, diagonal.penalty=diagonal.penalty,apply.by=apply.by)  Dirk Eddelbuettel committed Apr 10, 2018 1012 1013 1014  } else { ## it's a smooth with an id, so basis setup data differs from model matrix data names(id.list[[id]]$data) <- split$smooth.spec[[i]]$term ## give basis data suitable names sml <- smoothCon(split$smooth.spec[[i]],id.list[[id]]$data,knots,  Dirk Eddelbuettel committed Apr 10, 2018 1015  absorb.cons,n=nrow(data),dataX=data,scale.penalty=scale.penalty,  Dirk Eddelbuettel committed Apr 10, 2018 1016 1017  null.space.penalty=select,sparse.cons=sparse.cons, diagonal.penalty=diagonal.penalty,apply.by=apply.by)  Dirk Eddelbuettel committed Apr 10, 2018 1018 1019 1020 1021 1022  } for (j in 1:length(sml)) { newm <- newm + 1 sm[[newm]] <- sml[[j]] }  Dirk Eddelbuettel committed Apr 10, 2018 1023 1024  }  Dirk Eddelbuettel committed Apr 10, 2018 1025 1026  G$m <- m <- newm ## number of actual smooths  Dirk Eddelbuettel committed Apr 10, 2018 1027 1028  ## at this stage, it is neccessary to impose any side conditions required ## for identifiability  Dirk Eddelbuettel committed Apr 10, 2018 1029 1030 1031 1032 1033 1034 1035 1036 1037  if (m>0) { sm <- gam.side(sm,X,tol=.Machine$double.eps^.5) if (!apply.by) for (i in 1:length(sm)) { ## restore any by-free model matrices if (!is.null(sm[[i]]$X0)) { ## there is a by-free matrix to restore ind <- attr(sm[[i]],"del.index") ## columns, if any to delete sm[[i]]$X <- if (is.null(ind)) sm[[i]]$X0 else sm[[i]]$X0[,-ind,drop=FALSE] } } }  Dirk Eddelbuettel committed Apr 10, 2018 1038   Dirk Eddelbuettel committed Apr 10, 2018 1039 1040 1041 1042 1043  ## The matrix, L, mapping the underlying log smoothing parameters to the ## log of the smoothing parameter multiplying the S[[i]] must be ## worked out... idx <- list() ## idx[[id]]$c contains index of first col in L relating to id L <- matrix(0,0,0)  Dirk Eddelbuettel committed Apr 10, 2018 1044  lsp.names <- sp.names <- rep("",0) ## need a list of names to identify sps in global sp array  Dirk Eddelbuettel committed Apr 10, 2018 1045 1046 1047  if (m>0) for (i in 1:m) { id <- sm[[i]]$id ## get the L matrix for this smooth...  Dirk Eddelbuettel committed Apr 10, 2018 1048 1049  length.S <- length(sm[[i]]$S) if (is.null(sm[[i]]$L)) Li <- diag(length.S) else Li <- sm[[i]]$L  Dirk Eddelbuettel committed Apr 10, 2018 1050 1051 1052 1053 1054 1055 1056 1057 1058  if (length.S > 0) { ## there are smoothing parameters to name if (length.S == 1) spn <- sm[[i]]$label else { Sname <- names(sm[[i]]$S) if (is.null(Sname)) spn <- paste(sm[[i]]$label,1:length.S,sep="") else spn <- paste(sm[[i]]$label,Sname,sep="") } }  Dirk Eddelbuettel committed Apr 10, 2018 1059 1060 1061 1062 1063 1064 1065 1066  ## extend the global L matrix... if (is.null(id)||is.null(idx[[id]])) { ## new id' if (!is.null(id)) { ## create record in idx' idx[[id]]$c <- ncol(L)+1 ## starting column in L for this id' idx[[id]]$nc <- ncol(Li) ## number of columns relating to this id' } L <- rbind(cbind(L,matrix(0,nrow(L),ncol(Li))), cbind(matrix(0,nrow(Li),ncol(L)),Li))  Dirk Eddelbuettel committed Apr 10, 2018 1067 1068  if (length.S > 0) { ## there are smoothing parameters to name sp.names <- c(sp.names,spn) ## extend the sp name vector  Dirk Eddelbuettel committed Apr 10, 2018 1069  lsp.names <- c(lsp.names,spn) ## extend full.sp name vector  Dirk Eddelbuettel committed Apr 10, 2018 1070  }  Dirk Eddelbuettel committed Apr 10, 2018 1071 1072 1073 1074 1075 1076 1077  } else { ## it's a repeat id => shares existing sp's L0 <- matrix(0,nrow(Li),ncol(L)) if (ncol(Li)>idx[[id]]$nc) { stop("Later terms sharing an id' can not have more smoothing parameters than the first such term") } L0[,idx[[id]]$c:(idx[[id]]$c+ncol(Li)-1)] <- Li L <- rbind(L,L0)  Dirk Eddelbuettel committed Apr 10, 2018 1078 1079 1080  if (length.S > 0) { ## there are smoothing parameters to name lsp.names <- c(lsp.names,spn) ## extend full.sp name vector }  Dirk Eddelbuettel committed Apr 10, 2018 1081 1082 1083  } }  Dirk Eddelbuettel committed Apr 10, 2018 1084 1085 1086  ## create the model matrix... Xp <- NULL ## model matrix under prediction constraints, if given  Dirk Eddelbuettel committed Apr 10, 2018 1087 1088  if (m>0) for (i in 1:m) { n.para<-ncol(sm[[i]]$X)  Dirk Eddelbuettel committed Apr 10, 2018 1089  # define which elements in the parameter vector this smooth relates to....  Dirk Eddelbuettel committed Apr 10, 2018 1090  sm[[i]]$first.para<-first.para  Dirk Eddelbuettel committed Apr 10, 2018 1091  first.para<-first.para+n.para  Dirk Eddelbuettel committed Apr 10, 2018 1092  sm[[i]]$last.para<-first.para-1  Dirk Eddelbuettel committed Apr 10, 2018 1093 1094 1095 1096 1097 1098 1099  ## termwise offset handling ... Xoff <- attr(sm[[i]]$X,"offset") if (!is.null(Xoff)) { if (is.null(G$offset)) G$offset <- Xoff else G$offset <- G$offset + Xoff } ## model matrix accumulation ...  Dirk Eddelbuettel committed Apr 10, 2018 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109  ## alternative version under alternative constraint first (prediction only) if (is.null(sm[[i]]$Xp)) { if (!is.null(Xp)) Xp <- cbind2(Xp,sm[[i]]$X) } else { if (is.null(Xp)) Xp <- X Xp <- cbind2(Xp,sm[[i]]$Xp);sm[[i]]$Xp <- NULL } ## now version to use for fitting ... X <- cbind2(X,sm[[i]]$X);sm[[i]]$X<-NULL  Dirk Eddelbuettel committed Apr 10, 2018 1110   Dirk Eddelbuettel committed Apr 10, 2018 1111  G$smooth[[i]] <- sm[[i]]  Dirk Eddelbuettel committed Apr 10, 2018 1112  }  Dirk Eddelbuettel committed Apr 10, 2018 1113   Dirk Eddelbuettel committed Apr 10, 2018 1114 1115 1116 1117  if (is.null(Xp)) { G$cmX <- colMeans(X) ## useful for componentwise CI construction } else { G$cmX <- colMeans(Xp)  Dirk Eddelbuettel committed Apr 10, 2018 1118  ## transform from fit params to prediction params...  Dirk Eddelbuettel committed Apr 10, 2018 1119  ## G$P <- qr.coef(qr(Xp),X) ## old code assumes always full rank!!  Dirk Eddelbuettel committed Apr 10, 2018 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137  qrx <- qr(Xp,LAPACK=TRUE) R <- qr.R(qrx) p <- ncol(R) rank <- Rrank(R) ## rank of Xp/R QtX <- qr.qty(qrx,X)[1:rank,] if (rank0) G$cmX[-(1:G$nsdf)] <- 0 ## zero the smooth parts here else G$cmX <- G$cmX * 0  Dirk Eddelbuettel committed Apr 10, 2018 1141 1142  G$X <- X;rm(X) n.p <- ncol(G$X)  Dirk Eddelbuettel committed Apr 10, 2018 1143  # deal with penalties  Dirk Eddelbuettel committed Apr 10, 2018 1144   Dirk Eddelbuettel committed Apr 10, 2018 1145   Dirk Eddelbuettel committed Apr 10, 2018 1146 1147 1148 1149 1150 1151  ## min.sp must be length nrow(L) to make sense ## sp must be length ncol(L) --- need to partition ## L into columns relating to free log smoothing parameters, ## and columns, L0, corresponding to values supplied in sp. ## lsp0 = L0%*%log(sp[sp>=0]) [need to fudge sp==0 case by ## setting log(0) to log(effective zero) computed case-by-case]  Dirk Eddelbuettel committed Apr 10, 2018 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163  ## following deals with supplied and estimated smoothing parameters... ## first process the sp' array supplied to gam'... if (!is.null(sp)) # then user has supplied fixed smoothing parameters { if (length(sp)!=ncol(L)) { warning("Supplied smoothing parameter vector is too short - ignored.")} if (sum(is.na(sp))) { warning("NA's in supplied smoothing parameter vector - ignoring.")} G$sp <- sp } else { # set up for auto-initialization G$sp<-rep(-1,ncol(L)) }  Dirk Eddelbuettel committed Apr 10, 2018 1164 1165  names(G$sp) <- sp.names  Dirk Eddelbuettel committed Apr 10, 2018 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193  ## now work through the smooths searching for any sp' elements ## supplied in s' or te' terms.... This relies on idx' created ## above... k <- 1 ## current location in sp' array if (m>0) for (i in 1:m) { id <- sm[[i]]$id if (is.null(sm[[i]]$L)) Li <- diag(length(sm[[i]]$S)) else Li <- sm[[i]]$L if (is.null(id)) { ## it's a smooth without an id spi <- sm[[i]]$sp if (!is.null(spi)) { ## sp supplied in s' or te' if (length(spi)!=ncol(Li)) stop("incorrect number of smoothing parameters supplied for a smooth term") G$sp[k:(k+ncol(Li)-1)] <- spi } k <- k + ncol(Li) } else { ## smooth has an id spi <- sm[[i]]$sp if (is.null(idx[[id]]$sp.done)) { ## not already dealt with these sp's if (!is.null(spi)) { ## sp supplied in s' or te' if (length(spi)!=ncol(Li)) stop("incorrect number of smoothing parameters supplied for a smooth term") G$sp[idx[[id]]$c:(idx[[id]]$c+idx[[id]]$nc-1)] <- spi } idx[[id]]$sp.done <- TRUE ## only makes sense to use supplied sp' from defining term k <- k + idx[[id]]$nc } } } ## finished processing sp' vectors supplied in s' or te' terms  Dirk Eddelbuettel committed Apr 10, 2018 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215  ## copy initial sp's back into smooth objects, so there is a record of ## fixed and free... k <- 1 if (length(idx)) for (i in 1:length(idx)) idx[[i]]$sp.done <- FALSE if (m>0) for (i in 1:m) { ## work through all smooths id <- sm[[i]]$id if (!is.null(id)) { ## smooth with id if (idx[[id]]$nc>0) { ## only copy if there are sp's G$smooth[[i]]$sp <- G$sp[idx[[id]]$c:(idx[[id]]$c+idx[[id]]$nc-1)] } if (!idx[[id]]$sp.done) { ## only update k on first encounter with this smooth idx[[id]]$sp.done <- TRUE k <- k + idx[[id]]$nc } } else { ## no id, just work through sp if (is.null(sm[[i]]$L)) nc <- length(sm[[i]]$S) else nc <- ncol(sm[[i]]$L) if (nc>0) G$smooth[[i]]$sp <- G$sp[k:(k+nc-1)] k <- k + nc } } ## now all elements of G$smooth have a record of initial sp.  Dirk Eddelbuettel committed Apr 10, 2018 1216 1217 1218 1219 1220 1221 1222  if (!is.null(min.sp)) # then minimum s.p.'s supplied { if (length(min.sp)!=nrow(L)) stop("length of min.sp is wrong.") if (sum(is.na(min.sp))) stop("NA's in min.sp.") if (sum(min.sp<0)) stop("elements of min.sp must be non negative.") }  Dirk Eddelbuettel committed Apr 10, 2018 1223 1224 1225 1226  k.sp <- 0 # count through sp and S G$rank <- array(0,0) if (m>0) for (i in 1:m) { sm<-G$smooth[[i]]  Dirk Eddelbuettel committed Apr 10, 2018 1227  if (length(sm$S)>0)  Dirk Eddelbuettel committed Apr 10, 2018 1228 1229 1230 1231  for (j in 1:length(sm$S)) { # work through penalty matrices k.sp <- k.sp+1 G$off[k.sp] <- sm$first.para G$S[[k.sp]] <- sm$S[[j]]  Dirk Eddelbuettel committed Apr 10, 2018 1232  G$rank[k.sp]<-sm$rank[j]  Dirk Eddelbuettel committed Apr 10, 2018 1233 1234 1235  if (!is.null(min.sp)) { if (is.null(H)) H<-matrix(0,n.p,n.p) H[sm$first.para:sm$last.para,sm$first.para:sm$last.para] <-  Dirk Eddelbuettel committed Apr 10, 2018 1236 1237 1238 1239  H[sm$first.para:sm$last.para,sm$first.para:sm$last.para]+min.sp[k.sp]*sm$S[[j]] } } }  Dirk Eddelbuettel committed Apr 10, 2018 1240   Dirk Eddelbuettel committed Apr 10, 2018 1241  ## need to modify L, lsp.names, G$S, G$sp, G$rank and G$off to include any penalties  Dirk Eddelbuettel committed Apr 10, 2018 1242 1243 1244 1245 1246 1247 1248 1249  ## on parametric stuff, at this point.... if (!is.null(PP)) { ## deal with penalties on parametric terms L <- rbind(cbind(L,matrix(0,nrow(L),ncol(PP$L))), cbind(matrix(0,nrow(PP$L),ncol(L)),PP$L)) G$off <- c(PP$off,G$off) G$S <- c(PP$S,G$S) G$rank <- c(PP$rank,G$rank) G$sp <- c(PP$sp,G$sp)  Dirk Eddelbuettel committed Apr 10, 2018 1250  lsp.names <- c(PP$full.sp.names,lsp.names)  Dirk Eddelbuettel committed Apr 10, 2018 1251  G$n.paraPen <- length(PP$off)  Dirk Eddelbuettel committed Apr 10, 2018 1252 1253 1254 1255 1256 1257 1258  if (!is.null(PP$min.sp)) { ## deal with minimum sps if (is.null(H)) H <- matrix(0,n.p,n.p) for (i in 1:length(PP$S)) { ind <- PP$off[i]:(PP$off[i]+ncol(PP$S[[i]])-1) H[ind,ind] <- H[ind,ind] + PP$min.sp[i] * PP$S[[i]] } } ## min.sp stuff finished  Dirk Eddelbuettel committed Apr 10, 2018 1259  } else G$n.paraPen <- 0  Dirk Eddelbuettel committed Apr 10, 2018 1260 1261 1262 1263 1264 1265 1266 1267 1268  ## Now remove columns of L and rows of sp relating to fixed ## smoothing parameters, and use removed elements to create lsp0 fix.ind <- G$sp>=0 if (sum(fix.ind)) { lsp0 <- G$sp[fix.ind]  Dirk Eddelbuettel committed Apr 10, 2018 1269 1270 1271 1272 1273 1274 1275  ind <- lsp0==0 ## find the zero s.p.s ef0 <- indi <- (1:length(ind))[ind] if (length(indi)>0) for (i in 1:length(indi)) { ## find "effective zero" to replace each zero s.p. with ii <- G$off[i]:(G$off[i]+ncol(G$S[[i]])-1) ef0[i] <- norm(G$X[,ii],type="F")^2/norm(G$S[[i]],type="F")*.Machine$double.eps*.1 }  Dirk Eddelbuettel committed Apr 10, 2018 1276  lsp0[!ind] <- log