1## S3 method
2predict9 <- function (object,...) UseMethod("predict9")
3## Calculate prediction and Bayesian SE from ssanova objects
4predict9.gssanova <- function(object,newdata,ci=FALSE,level=.95,nu=NULL,...)
5{
6    est <- predict(object,newdata,se.fit=ci)
7    if (object$family=="binomial") {
8        if (!ci) return(plogis(est))
9        else {
10            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,-1,1))
11            z <- plogis(z)
12            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
13        }
14    }
15    if (object$family%in%c("poisson","Gamma","inverse.gaussian")) {
16        if (!ci) return(exp(est))
17        else {
18            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,-1,1))
19            z <- exp(z)
20            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
21        }
22    }
23    if (object$family=="nbinomial") {
24        if (!is.vector(model.response(object$mf))) {
25            if (is.null(nu))
26                stop("gss error: nu is missing for nbinomial family with 2 column response" )
27        }
28        else nu <- object$nu
29        if (!ci) return(nu/exp(est))
30        else {
31            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,1,-1))
32            z <- nu/exp(z)
33            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
34        }
35    }
36    if (object$family=="weibull") {
37        gg <- gamma(1+1/object$nu)
38        if (!ci) return(gg*exp(est))
39        else {
40            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,-1,1))
41            z <- gg*exp(z)
42            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
43        }
44    }
45    if (object$family=="lognorm") {
46        gg <- exp(1/2/object$nu^2)
47        if (!ci) return(gg*exp(est))
48        else {
49            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,-1,1))
50            z <- gg*exp(z)
51            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
52        }
53    }
54    if (object$family=="loglogis") {
55        gg <- pi/object$nu/sin(pi/object$nu)
56        if (!ci) return(gg*exp(est))
57        else {
58            z <- est$fit+outer(est$se.fit*qnorm(1-(1-level)/2),c(0,-1,1))
59            z <- gg*exp(z)
60            return(list(fit=z[,1],lcl=z[,2],ucl=z[,3]))
61        }
62    }
63    if (object$family=="polr") {
64        P <- plogis(outer(est,c(0,cumsum(object$nu)),"+"))
65        z <- P[,1]
66        J <- dim(P)[2]
67        for (i in 2:J) z <- cbind(z,P[,i]-P[,i-1])
68        z <- cbind(z,1-P[,J])
69        colnames(z) <- NULL
70        return(z)
71    }
72}
73