1## Calculate prediction and Bayesian SE from ssllrm objects
2predict.ssllrm <- function (object,x,y=object$qd.pt,odds=NULL,se.odds=FALSE,...)
3{
4    if (class(object)!="ssllrm")
5        stop("gss error in predict.ssllrm: not a ssllrm object")
6    if ("random"%in%colnames(x)) {
7        zz <- x$random
8        x$random <- NULL
9    }
10    else zz <- NULL
11    if (!all(sort(object$xnames)==sort(colnames(x))))
12        stop("gss error in predict.ssllrm: mismatched x variable names")
13    if (!all(sort(object$ynames)==sort(colnames(y))))
14        stop("gss error in predict.ssllrm: mismatched y variable names")
15    mf <- object$mf
16    term <- object$term
17    qd.pt <- object$qd.pt
18    qd.wt <- object$qd.wt
19    nmesh <- dim(qd.pt)[1]
20    y.id <- NULL
21    for (i in 1:dim(y)[1]) {
22        if (!sum(duplicated(rbind(qd.pt,y[i,object$ynames,drop=FALSE]))))
23            stop("gss error in predict.ssllrm: y value is out of range")
24        wk <- FALSE
25        for (j in 1:nmesh) {
26            if (sum(duplicated(rbind(qd.pt[j,],y[i,object$ynames])))) y.id <- c(y.id,j)
27        }
28    }
29    if (!is.null(odds)) {
30        if (length(y.id)-length(odds))
31            stop("gss error in predict.ssllrm: odds is of wrong length")
32        if (!max(odds)|sum(odds))
33            stop("gss error in predict.ssllrm: odds is not a contrast")
34        if (sum(duplicated(y.id)))
35            stop("gss error in predict.ssllrm: duplicated y in contrast")
36        qd.pt <- qd.pt[y.id,,drop=FALSE]
37    }
38    ## Generate s, and r
39    nobs <- dim(x)[1]
40    nmesh <- dim(qd.pt)[1]
41    nbasis <- length(object$id.basis)
42    nnull <- length(object$d)
43    nZ <- length(object$b)
44    s <- NULL
45    r <- array(0,c(nmesh,nbasis,nobs))
46    nu <- nq <- 0
47    for (label in term$labels) {
48        vlist <- term[[label]]$vlist
49        x.list <- object$xnames[object$xnames%in%vlist]
50        y.list <- object$ynames[object$ynames%in%vlist]
51        xy.basis <- mf[object$id.basis,vlist]
52        qd.xy <- data.frame(matrix(0,nmesh,length(vlist)))
53        names(qd.xy) <- vlist
54        qd.xy[,y.list] <- qd.pt[,y.list]
55        if (length(x.list)) xx <- x[,x.list,drop=FALSE]
56        else xx <- NULL
57        nphi <- term[[label]]$nphi
58        nrk <- term[[label]]$nrk
59        if (nphi) {
60            phi <- term[[label]]$phi
61            for (i in 1:nphi) {
62                nu <- nu+1
63                if (is.null(xx)) {
64                    s.wk <- phi$fun(qd.xy[,,drop=TRUE],nu=i,env=phi$env)
65                    wk <- matrix(s.wk,nmesh,nobs)
66                }
67                else {
68                    wk <- NULL
69                    for (j in 1:nobs) {
70                        qd.xy[,x.list] <- xx[rep(j,nmesh),]
71                        wk <- cbind(wk,phi$fun(qd.xy,i,phi$env))
72                    }
73                }
74                s <- array(c(s,wk),c(nmesh,nobs,nu))
75            }
76        }
77        if (nrk) {
78            rk <- term[[label]]$rk
79            for (i in 1:nrk) {
80                nq <- nq+1
81                if (is.null(xx)) {
82                    r.wk <- rk$fun(qd.xy[,,drop=TRUE],xy.basis,nu=i,env=rk$env,out=TRUE)
83                    r <- r + as.vector(10^object$theta[nq]*r.wk)
84                }
85                else {
86                    wk <- NULL
87                    for (j in 1:nobs) {
88                        qd.xy[,x.list] <- xx[rep(j,nmesh),]
89                        wk <- array(c(wk,rk$fun(qd.xy,xy.basis,i,rk$env,TRUE)),
90                                    c(nmesh,nbasis,j))
91                    }
92                    r <- r + 10^object$theta[nq]*wk
93                }
94            }
95        }
96    }
97    ## random effects
98    if (nZ) {
99        nz <- object$Random$sigma$env$nz
100        if (is.null(zz)) z.wk <- matrix(0,nobs,nz)
101        else z.wk <- as.matrix(zz)
102        if (dim(z.wk)[2]!=nz)
103            stop("gss error in predict.ssllrm: x$random is of wrong dimension")
104        z <- nlvl <- NULL
105        for (ylab in object$ynames) {
106            y.wk <- mf[,ylab]
107            lvl.wk <- levels(y.wk)
108            nlvl.wk <- length(lvl.wk)
109            nlvl <- c(nlvl,nlvl.wk)
110            z.aux <- diag(1,nlvl.wk-1)
111            z.aux <- rbind(z.aux,rep(-1,nlvl.wk-1))
112            rownames(z.aux) <- lvl.wk
113            pt.wk <- qd.pt[,ylab]
114            for (i in 1:(nlvl.wk-1)) {
115                for (j in 1:nmesh) {
116                    z <- cbind(z,z.aux[pt.wk[j],i]*z.wk)
117                }
118            }
119        }
120        z <- aperm(array(z,c(nobs,nz,nmesh,nZ/nz)),c(3,2,4,1))
121        z <- array(z,c(nmesh,nZ,nobs))
122    }
123    ## return
124    if (is.null(odds)) {
125        pdf <- NULL
126        for (j in 1:nobs) {
127            wk <- matrix(r[,,j],nmesh,nbasis)%*%object$c
128            if (nnull) wk <- wk + matrix(s[,j,],nmesh,nnull)%*%object$d
129            if (nZ) wk <- wk + matrix(z[,,j],nmesh,nZ)%*%object$b
130            wk <- exp(wk)*qd.wt
131            pdf <- cbind(pdf,wk/sum(wk))
132        }
133        return(t(pdf[y.id,]))
134    }
135    else {
136        s.wk <- r.wk <- z.wk <- w.wk <- 0
137        for (i in 1:length(odds)) {
138            r.wk <- r.wk + odds[i]*r[i,,]
139            if (nnull) s.wk <- s.wk + odds[i]*s[i,,]
140            if (nZ) z.wk <- z.wk + odds[i]*z[i,,]
141            w.wk <- w.wk + odds[i]*log(qd.wt[y.id[i]])
142        }
143        s.wk <- matrix(s.wk,nobs,nnull)
144        r.wk <- t(matrix(r.wk,nbasis,nobs))
145        z.wk <- t(matrix(z.wk,nZ,nobs))
146        rs <- cbind(r.wk,z.wk,s.wk)
147        if (!se.odds) as.vector(rs%*%c(object$c,object$b,object$d))
148        else {
149            fit <- as.vector(rs%*%c(object$c,object$b,object$d)) + w.wk
150            se.fit <- .Fortran("hzdaux2",
151                               as.double(object$se.aux$v), as.integer(dim(rs)[2]),
152                               as.integer(object$se.aux$jpvt),
153                               as.double(t(rs)), as.integer(dim(rs)[1]),
154                               se=double(dim(rs)[1]), PACKAGE="gss")[["se"]]
155            return(list(fit=fit,se.fit=se.fit))
156        }
157    }
158}
159