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