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