Nothing
#' @title Predict method for \pkg{bcgam} fits
#'
#' @description Predicted values based on \pkg{bcgam} object.
#'
#' @param object Object of class inheriting from "bcgam".
#' @param newdata An optional data frame in which to look for variable with which to predict.
#' @param interval Type of interval calculation. It can be either \code{"credible"} or \code{"prediction"}.
#' If \code{interval=prediction}, then a summary of the predictive posterior distribution is shown.
#' @param level Tolerance/credible level. The default is \code{0.95}.
#' @param parameter The type of parameter to predict. If parameter=\code{"eta"}, then the systematic component
#' \eqn{\eta} is predicted. If parameter=\code{"mu"}, then the mean value \eqn{\mu} obtained by transforming
#' \eqn{\eta} using the inverse of the link function is predicted. The default is \eqn{\mu}. If
#' \code{interval="prediction"}, then this variable is ignored.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#'
#' @details \code{predict.bcgam} produces estimated values, standard errors,
#' and interval estimates; obtained by the
#' fitted model in the frame \code{newdata}. Setting \code{intervals}
#' specifies computation of credible or prediction (tolerance) intervals at the
#' specified \code{level}.
#'
#' All predictions are based on the posterior distribution in the \code{bcgam} object.
#'
#'
#' @return \code{predict.bcgam} produces a list of predictions, standard errors, and bounds based
#' on the posterior disitribution in the \code{bcgam} object.
#'
#' If \code{interval="credible"}, a list with the following components is returned:
#' \item{cred.mean}{the mean values.}
#' \item{cred.sd}{standard error of mean values.}
#' \item{cred.lower}{lower bound of mean values.}
#' \item{cred.upper}{upper bound of mean values.}
#' If \code{interval="prediction"}, a list with the following components is returned:
#' \item{pred.mean}{the predicted values.}
#' \item{pred.sd}{standard error of predicted values.}
#' \item{pred.lower}{lower bound of predicted values.}
#' \item{pred.upper}{upper bound of predicted values.}
#'
#' @author Cristian Oliva-Aviles and Mary C. Meyer
#'
#' @keywords internal
#'
#' @examples
#' \dontrun{
#' n<-50
#' x<-(1:n)^{1/3}
#' z<-as.factor(rbinom(n, 1, 0.6))
#' y<-x+7*as.numeric(z)+rnorm(n,sd=2)
#' bcgam.fit <- bcgam(y~sm.incr(x)+z, nloop=100)
#' predict(bcgam.fit, newdata=data.frame(x=0.5, z="1"), parameter="mu", interval="prediction")
#' }
predict.bcgam <- function(object, newdata, interval=c("credible"),
level=0.95, parameter=c("mu"), ...)
{
tt <- terms(object)
if (!inherits(object, "bcgam"))
warning("calling predict.bcgam(<fake-bcgam-object>) ...")
if (missing(newdata) || is.null(newdata)) {
mm <- X <- model.matrix(object)
mmDone <- TRUE
offset <- object$offsetn
}
else {
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.pass,
xlev = object$levels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
offset <- rep(0, nrow(X))
if (!is.null(off.num <- attr(tt, "offset")))
for (i in off.num) offset <- offset + eval(attr(tt,
"variables")[[i + 1]], newdata)
if (!is.null(object$call$offset))
offset <- offset + eval(object$call$offset, newdata)
mmDone <- FALSE
}
n=length(object$y)
L=length(object$xmat)/n
delta.newdata<-NULL
incr<-NULL
decr<-NULL
zmatb.newdata<-NULL
xmat.newdata<-cbind(X)
no.rows=dim(newdata)[1]
v=1:no.rows*0+1
for(i in 1:L){
Xi=X[,object$ind_nonparam[i]]
if(object$shapes[i]==1){
delta.newdata=rbind( delta.newdata,monincr(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==2){
delta.newdata=rbind( delta.newdata,mondecr(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==3){ delta.i=convex(Xi, object$knots[[i]], pred.new=TRUE)
delta.newdata=rbind( delta.newdata,delta.i$sigma-t(delta.i$x.mat%*%object$center.delta[[i]]) )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==4){ delta.i=concave(Xi, object$knots[[i]], pred.new=TRUE)
delta.newdata=rbind( delta.newdata,delta.i$sigma-t(delta.i$x.mat%*%object$center.delta[[i]]) )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==5){
delta.newdata=rbind( delta.newdata,incconvex(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==7){
delta.newdata=rbind( delta.newdata,incconcave(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==6){
delta.newdata=rbind( delta.newdata,decconvex(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==8){
delta.newdata=rbind( delta.newdata,decconcave(Xi,object$knots[[i]])$sigma-object$center.delta[[i]] )} # delta.newdata contains basis splines in rows
if(object$shapes[i]==1|object$shapes[i]==5|object$shapes[i]==7){incr=c(incr,1)}else{incr=c(incr,0)}
if(object$shapes[i]==2|object$shapes[i]==6|object$shapes[i]==8){decr=c(decr,1)}else{decr=c(decr,0)}
if(incr[i]==0 & decr[i]==0){
zmatb.newdata=cbind(zmatb.newdata,Xi)
}
}
ind.nonparaplusinter=c(1, object$ind_nonparam)
if(no.rows==1){
if(object$ind_intercept==1){
zmatb.newdata=c(v,xmat.newdata[,-(ind.nonparaplusinter)],zmatb.newdata) }
else{zmatb.newdata=c(xmat.newdata[,-(ind.nonparaplusinter)],zmatb.newdata)}
deltazmat=c(t(delta.newdata),zmatb.newdata)
pred.sims=object$coefs.sims%*%deltazmat
}
else{
if(object$ind_intercept==1){
zmatb.newdata=cbind(v,xmat.newdata[,-(ind.nonparaplusinter)],zmatb.newdata) }
else{zmatb.newdata=cbind(xmat.newdata[,-(ind.nonparaplusinter)],zmatb.newdata)}
deltazmat=cbind(t(delta.newdata),zmatb.newdata)
pred.sims=object$coefs.sims%*%t(deltazmat)
}
alpha.upper=1-(1-level)/2
alpha.lower=(1-level)/2
q.upper=function(x){quantile(x,alpha.upper)}
q.lower=function(x){quantile(x,alpha.lower)}
get.predgaussian=function(x){rnorm( rep(1,length(x)) ,x , object$sigma.sims)}
get.predbinomial=function(x){rbinom( rep(1,length(x)), rep(1,length(x)), x) }
get.predpoisson=function(x){rpois( rep(1,length(x)), x) }
#freqtoprop=function(x){prop.table(table(x))}
if(object$family=="gaussian"){
if(interval=="credible"|interval=="prediction"){
cred.mean=apply(pred.sims,2,mean)
if(interval=="credible"){
if(parameter=="eta"){
cred.sd=apply(pred.sims,2,sd)
cred.upper=apply(pred.sims,2,q.upper)
cred.lower=apply(pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated Mean")
}
else{
if(parameter=="mu"){
cred.sd=apply(pred.sims,2,sd)
cred.upper=apply(pred.sims,2,q.upper)
cred.lower=apply(pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated Mean")
}
else{stop("parameter not valid!")}
}
}
else{
pred.gaussian=apply(pred.sims, 2, get.predgaussian)
pred.mean=apply(pred.gaussian,2,mean)
pred.sd=apply(pred.gaussian,2,sd)
pred.upper=apply(pred.gaussian,2,q.upper)
pred.lower=apply(pred.gaussian,2,q.lower)
list(pred.mean=pred.mean, pred.sd=pred.sd,
pred.lower=pred.lower, pred.upper=pred.upper,
y.lab="Prediction Mean")
}
}
else{stop("incorrect type of interval!")}
} #end if gaussian
else{
if(object$family=="binomial"){
if(interval=="credible"|interval=="prediction"){
mu.pred.sims=exp(pred.sims)/(1+exp(pred.sims))
if(interval=="credible"){
if(parameter=="eta"){
cred.mean=apply(pred.sims,2,mean)
cred.sd=apply(pred.sims,2,sd)
cred.upper=apply(pred.sims,2,q.upper)
cred.lower=apply(pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated log(odds)")
}
else{
if(parameter=="mu"){
cred.mean=apply(mu.pred.sims,2,mean)
cred.sd=apply(mu.pred.sims,2,sd)
cred.upper=apply(mu.pred.sims,2,q.upper)
cred.lower=apply(mu.pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated probability")
}
else{stop("parameter not valid!")}
}
}
else{
pred.binomial<-apply(mu.pred.sims, 2, get.predbinomial)
pred.mean<-apply(pred.binomial,2,mean)
pred.sd<-apply(pred.binomial,2,sd)
list(pred.mean=pred.mean, pred.sd=pred.sd,
y.lab="Prediction probability")
}
}
else{stop("incorrect type of interval!")}
} #end if binomial
else{
if(object$family=="poisson"){
if(interval=="credible"|interval=="prediction"){
mu.pred.sims=exp(pred.sims)
if(interval=="credible"){
if(parameter=="eta"){
cred.mean=apply(pred.sims,2,mean)
cred.sd=apply(pred.sims,2,sd)
cred.upper=apply(pred.sims,2,q.upper)
cred.lower=apply(pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated log(counts)")
}
else{
if(parameter=="mu"){
cred.mean=apply(mu.pred.sims,2,mean)
cred.sd=apply(mu.pred.sims,2,sd)
cred.upper=apply(mu.pred.sims,2,q.upper)
cred.lower=apply(mu.pred.sims,2,q.lower)
list(cred.mean=cred.mean, cred.sd=cred.sd,
cred.lower=cred.lower, cred.upper=cred.upper,
y.lab="Estimated counts")
}
else{stop("parameter not valid!")}
}
}
else{
pred.poisson<-apply(mu.pred.sims, 2, get.predpoisson)
pred.mean<-apply(pred.poisson,2,mean)
pred.sd<-apply(pred.poisson,2,sd)
list(pred.mean=pred.mean, pred.sd=pred.sd,
y.lab="Prediction counts")
}
}
else{stop("incorrect type of interval!")}
}
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.