1## Fit ssanova model 2ssanova0 <- function(formula,type=NULL,data=list(),weights,subset, 3 offset,na.action=na.omit,partial=NULL, 4 method="v",varht=1,prec=1e-7,maxiter=30) 5{ 6 ## Obtain model frame and model terms 7 mf <- match.call() 8 mf$type <- mf$method <- mf$varht <- mf$partial <- NULL 9 mf$prec <- mf$maxiter <- NULL 10 mf[[1]] <- as.name("model.frame") 11 mf <- eval(mf,parent.frame()) 12 ## Generate terms 13 term <- mkterm(mf,type) 14 ## Generate s, q, and y 15 nobs <- dim(mf)[1] 16 s <- q <- NULL 17 nq <- 0 18 for (label in term$labels) { 19 if (label=="1") { 20 s <- cbind(s,rep(1,len=nobs)) 21 next 22 } 23 x <- mf[,term[[label]]$vlist] 24 nphi <- term[[label]]$nphi 25 nrk <- term[[label]]$nrk 26 if (nphi) { 27 phi <- term[[label]]$phi 28 for (i in 1:nphi) 29 s <- cbind(s,phi$fun(x,nu=i,env=phi$env)) 30 } 31 if (nrk) { 32 rk <- term[[label]]$rk 33 for (i in 1:nrk) { 34 nq <- nq+1 35 q <- array(c(q,rk$fun(x,x,nu=i,env=rk$env,out=TRUE)),c(nobs,nobs,nq)) 36 } 37 } 38 } 39 ## Add the partial term 40 if (!is.null(partial)) { 41 mf.p <- model.frame(partial,data) 42 for (lab in colnames(mf.p)) mf[,lab] <- mf.p[,lab] 43 mt.p <- attr(mf.p,"terms") 44 lab.p <- labels(mt.p) 45 matx.p <- model.matrix(mt.p,data)[,-1,drop=FALSE] 46 if (dim(matx.p)[1]!=dim(mf)[1]) 47 stop("gss error in ssanova: partial data are of wrong size") 48 matx.p <- scale(matx.p) 49 center.p <- attr(matx.p,"scaled:center") 50 scale.p <- attr(matx.p,"scaled:scale") 51 s <- cbind(s,matx.p) 52 part <- list(mt=mt.p,center=center.p,scale=scale.p) 53 } 54 else part <- lab.p <- NULL 55 ## Prepare the data 56 y <- model.response(mf,"numeric") 57 w <- model.weights(mf) 58 offset <- model.offset(mf) 59 if (!is.null(offset)) { 60 term$labels <- c(term$labels,"offset") 61 term$offset <- list(nphi=0,nrk=0) 62 y <- y - offset 63 } 64 if (!is.null(w)) { 65 w <- sqrt(w) 66 y <- w*y 67 s <- w*s 68 for (i in 1:nq) q[,,i] <- w*t(w*q[,,i]) 69 } 70 if (qr(s)$rank<dim(s)[2]) 71 stop("gss error in ssanova0: unpenalized terms are linearly dependent") 72 if (!nq) stop("gss error in ssanova0: use lm for models with only unpenalized terms") 73 ## Fit the model 74 if (nq==1) { 75 q <- q[,,1] 76 z <- sspreg0(s,q,y,method,varht) 77 } 78 else z <- mspreg0(s,q,y,method,varht,prec,maxiter) 79 ## Brief description of model terms 80 desc <- NULL 81 for (label in term$labels) 82 desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")]))) 83 if (!is.null(partial)) { 84 desc <- rbind(desc,matrix(c(1,0),length(lab.p),2,byrow=TRUE)) 85 } 86 desc <- rbind(desc,apply(desc,2,sum)) 87 if (is.null(partial)) rownames(desc) <- c(term$labels,"total") 88 else rownames(desc) <- c(term$labels,lab.p,"total") 89 colnames(desc) <- c("Unpenalized","Penalized") 90 ## Return the results 91 obj <- c(list(call=match.call(),mf=mf,terms=term,partial=part,lab.p=lab.p, 92 desc=desc),z) 93 class(obj) <- c("ssanova0","ssanova") 94 obj 95} 96 97## Fit Single Smoothing Parameter REGression 98sspreg0 <- function(s,q,y,method="v",varht=1) 99{ 100 ## Check inputs 101 if (is.vector(s)) s <- as.matrix(s) 102 if (!(is.matrix(s)&is.matrix(q)&is.vector(y)&is.character(method))) { 103 stop("gss error in sspreg: inputs are of wrong types") 104 } 105 nobs <- length(y) 106 nnull <- dim(s)[2] 107 if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs) 108 &(nobs>=nnull)&(nnull>0))) { 109 stop("gss error in sspreg: inputs have wrong dimensions") 110 } 111 ## Set method for smoothing parameter selection 112 code <- (1:3)[c("v","m","u")==method] 113 if (!length(code)) { 114 stop("gss error: unsupported method for smoothing parameter selection") 115 } 116 ## Call RKPACK driver DSIDR 117 z <- .Fortran("dsidr0", 118 as.integer(code), 119 swk=as.double(s), as.integer(nobs), 120 as.integer(nobs), as.integer(nnull), 121 as.double(y), 122 qwk=as.double(q), as.integer(nobs), 123 as.double(0), as.integer(0), double(2), 124 nlambda=double(1), score=double(1), varht=as.double(varht), 125 c=double(nobs), d=double(nnull), 126 qraux=double(nnull), jpvt=integer(nnull), 127 double(3*nobs), 128 info=integer(1),PACKAGE="gss") 129 ## Check info for error 130 if (info<-z$info) { 131 if (info>0) 132 stop("gss error in sspreg: matrix s is rank deficient") 133 if (info==-2) 134 stop("gss error in sspreg: matrix q is indefinite") 135 if (info==-1) 136 stop("gss error in sspreg: input data have wrong dimensions") 137 if (info==-3) 138 stop("gss error in sspreg: unknown method for smoothing parameter selection.") 139 } 140 ## Return the fit 141 c(list(method=method,theta=0), 142 z[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")]) 143} 144 145## Fit Multiple Smoothing Parameter REGression 146mspreg0 <- function(s,q,y,method="v",varht=1,prec=1e-7,maxiter=30) 147{ 148 ## Check inputs 149 if (is.vector(s)) s <- as.matrix(s) 150 if (!(is.matrix(s)&is.array(q)&(length(dim(q))==3) 151 &is.vector(y)&is.character(method))) { 152 stop("gss error in mspreg: inputs are of wrong types") 153 } 154 nobs <- length(y) 155 nnull <- dim(s)[2] 156 nq <- dim(q)[3] 157 if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs) 158 &(nobs>=nnull)&(nnull>0)&(nq>1))) { 159 stop("gss error in mspreg: inputs have wrong dimensions") 160 } 161 ## Set method for smoothing parameter selection 162 code <- (1:3)[c("v","m","u")==method] 163 if (!length(code)) { 164 stop("gss error: unsupported method for smoothing parameter selection") 165 } 166 ## Call RKPACK driver DMUDR 167 z <- .Fortran("dmudr0", 168 as.integer(code), 169 as.double(s), # s 170 as.integer(nobs), as.integer(nobs), as.integer(nnull), 171 as.double(q), # q 172 as.integer(nobs), as.integer(nobs), as.integer(nq), 173 as.double(y), # y 174 as.double(0), as.integer(0), 175 as.double(prec), as.integer(maxiter), 176 theta=double(nq), nlambda=double(1), 177 score=double(1), varht=as.double(varht), 178 c=double(nobs), d=double(nnull), 179 integer(nnull+nq), 180 double(nobs*nobs*(nq+2)), 181 info=integer(1),PACKAGE="gss")[c("theta","info")] 182 ## Check info for error 183 if (info<-z$info) { 184 if (info>0) 185 stop("gss error in mspreg: matrix s is rank deficient") 186 if (info==-2) 187 stop("gss error in mspreg: matrix q is indefinite") 188 if (info==-1) 189 stop("gss error in mspreg: input data have wrong dimensions") 190 if (info==-3) 191 stop("gss error in mspreg: unknown method for smoothing parameter selection.") 192 if (info==-4) 193 stop("gss error in mspreg: iteration fails to converge, try to increase maxiter") 194 if (info==-5) 195 stop("gss error in mspreg: iteration fails to find a reasonable descent direction") 196 } 197 qwk <- 10^z$theta[1]*q[,,1] 198 for (i in 2:nq) qwk <- qwk + 10^z$theta[i]*q[,,i] 199 ## Call RKPACK driver DSIDR 200 zz <- .Fortran("dsidr0", 201 as.integer(code), 202 swk=as.double(s), as.integer(nobs), 203 as.integer(nobs), as.integer(nnull), 204 as.double(y), 205 qwk=as.double(qwk), as.integer(nobs), 206 as.double(0), as.integer(0), double(2), 207 nlambda=double(1), score=double(1), varht=as.double(varht), 208 c=double(nobs), d=double(nnull), 209 qraux=double(nnull), jpvt=integer(nnull), 210 double(3*nobs), 211 info=integer(1),PACKAGE="gss") 212 ## Return the fit 213 c(list(method=method,theta=z$theta), 214 zz[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")]) 215} 216 217## Obtain c & d for new y's 218getcrdr <- function(obj,r) 219{ 220 ## Check inputs 221 if (is.vector(r)) r <- as.matrix(r) 222 if (!(any(class(obj)=="ssanova0")&is.matrix(r))) { 223 stop("gss error in getcrdr: inputs are of wrong types") 224 } 225 nobs <- length(obj$c) 226 nnull <- length(obj$d) 227 nr <- dim(r)[2] 228 if (!((dim(r)[1]==nobs)&(nr>0))) { 229 stop("gss error in getcrdr: inputs have wrong dimensions") 230 } 231 ## Call RKPACK ulitity DCRDR 232 z <- .Fortran("dcrdr", 233 as.double(obj$swk), as.integer(nobs), 234 as.integer(nobs), as.integer(nnull), 235 as.double(obj$qraux), as.integer(obj$jpvt), 236 as.double(obj$qwk), as.integer(nobs), 237 as.double(obj$nlambda), 238 as.double(r), as.integer(nobs), as.integer(nr), 239 cr=double(nobs*nr), as.integer(nobs), 240 dr=double(nnull*nr), as.integer(nnull), 241 double(2*nobs), integer(1),PACKAGE="gss")[c("cr","dr")] 242 ## Return cr and dr 243 z$cr <- matrix(z$cr,nobs,nr) 244 z$dr <- matrix(z$dr,nnull,nr) 245 z 246} 247 248## Obtain var-cov matrix for unpenalized terms 249getsms <- function(obj) 250{ 251 ## Check input 252 if (!any(class(obj)=="ssanova0")) { 253 stop("gss error in getsms: inputs are of wrong types") 254 } 255 nobs <- length(obj$c) 256 nnull <- length(obj$d) 257 ## Call RKPACK ulitity DSMS 258 z <- .Fortran("dsms", 259 as.double(obj$swk), as.integer(nobs), 260 as.integer(nobs), as.integer(nnull), 261 as.integer(obj$jpvt), 262 as.double(obj$qwk), as.integer(nobs), 263 as.double(obj$nlambda), 264 sms=double(nnull*nnull), as.integer(nnull), 265 double(2*nobs), integer(1),PACKAGE="gss")["sms"] 266 ## Return the nnull-by-nnull matrix 267 matrix(z$sms,nnull,nnull) 268} 269