predict.plbpsm: Prediction from fitted PLBPSM model

In funstatpackages/GgAM: Generalized geoAdditive Models

Description

Takes a fitted plbpsm object produced by plbpsm() and produces predictions given a new set of values for the model covariates or the original values used for the model fit. Predictions can be accompanied by standard errors, based on the posterior distribution of the model coefficients. The routine can optionally return the matrix by which the model coefficients must be pre-multiplied in order to yield the values of the linear predictor at the supplied covariate values.

Usage

## S3 method for class 'plbpsm'
predict(object, newdata, type = "response",
  se.fit = FALSE, terms = NULL, exclude = NULL, block.size = NULL,
  newdata.guaranteed = FALSE, na.action = na.pass,
  unconditional = FALSE, newB = NULL, newind00 = NULL,
  backfitting = object$backfitting, ...)

Arguments

object	a fitted plbpsm object as produced by plbpsm()
newdata	A data frame or list containing the values of the model covariates at which predictions are required. If this is not provided then predictions corresponding to the original data are returned. If newdata is provided then it should contain all the variables needed for prediction: a warning is generated if not.
type	the type of prediction required. The default is on the scale of the linear predictors; the alternative "response" is on the scale of the response variable. Thus for a default binomial model the default predictions are of log-odds (probabilities on logit scale) and type = "response" gives the predicted probabilities. The "terms" option returns a matrix giving the fitted values of each term in the model formula on the linear predictor scale. When this has the value "link" the linear predictor (possibly with associated standard errors) is returned.
se.fit	when this is TRUE (not default) standard error estimates are returned for each prediction.
terms	if type=="terms" then only results for the terms given in this array will be returned.
exclude	if type=="terms" or type="iterms" then terms (smooth or parametric) named in this array will not be returned. Otherwise any smooth terms named in this array will be set to zero. If NULL then no terms are excluded. Note that this is the term names as it appears in the model summary, see example.
block.size	maximum number of predictions to process per call to underlying code: larger is quicker, but more memory intensive. Set to < 1 to use total number of predictions as this. If NULL then block size is 1000 if new data supplied, and the number of rows in the model frame otherwise.
newdata.guaranteed	Set to TRUE to turn off all checking of newdata: this can speed things up for large prediction tasks, but newdata must be complete, with no NA values for predictors required in the model.
na.action	what to do about NA values in newdata. With the default na.pass, any row of newdata containing NA values for required predictors, gives rise to NA predictions (even if the term concerned has no NA predictors). na.exclude or na.omit result in the dropping of newdata rows, if they contain any NA values for required predictors. If newdata is missing then NA handling is determined from object$na.action.
unconditional	if TRUE then the smoothing parameter uncertainty corrected covariance matrix is used to compute uncertainty bands, if available. Otherwise the bands treat the smoothing parameters as fixed.
newB	the given matrix of bivariate spline basis functions.
newind00	the given index of the data points in the triangulation.
backfitting	whether SBL estimation is obtained.
...	other arguments.

Details

See examples to see usages of different types.

Value

if se.fit is TRUE then a 2 item list is returned with items (both arrays) fit and se.fit containing predictions and associated standard error estimates, otherwise an array of predictions is returned. The dimensions of the returned arrays depends on whether type is "terms" or not: if it is then the array is 2 dimensional with each term in the linear predictor separate, otherwise the array is 1 dimensional and contains the linear predictor/predicted values (or corresponding s.e.s). The linear predictor returned termwise will not include the offset or the intercept.

Examples

library(MASS)
library(grpreg)
library(BPST)
data("eg1pop_dat")
eg1_V2 <- eg1pop_dat[['V2']]
eg1_T2 <- eg1pop_dat[['T2']]
eg1pop_rho03 <- eg1pop_dat[['rho03']]
sam <- eg1pop_rho03[sample(1:dim(eg1pop_rho03)[1],100),]

### Partial Linear Spatial Model ###
formula_d4 <- Y~z1+z2+z3+z4+z5+z6+z7+z8+b(x1,x2,d=4,r=1,V=eg1_V2,Tr=eg1_T2)
res <- plbpsm(formula=formula_d4,data=as.data.frame(sam),VS=TRUE)
# check deviance and the estimated error from predict function
newdata <- as.data.frame(sam[sample(1:dim(sam)[1],90),])
pred <- predict(res,newdata)

# exclude one covariate
pred0 <- predict(res,newdata,exclude="z2",type='terms')
# check

### Generalized Partially Linear Spatial Model ###
# Poisson family
data("eg_poi_pl")
sam <- as.data.frame(eg_poi[['sam_poi']])
sam[1,]
V <- eg_poi[['V1']]
Tr <- eg_poi[['T1']]
formula <- y~c1+c2+c3+c4+c5+c6+c7+c8+c9+c10+c11+c12+c13+c14+c15+b(loc1,loc2,V=V,Tr=Tr,d=2,r=1)
res_poi <- plbpsm(formula=formula,family=poisson(),data=as.data.frame(sam))

# with offset
res_poi2 <- plbpsm(formula=formula,family=poisson(),offset=rep(0.1,length(data$y)),
data <- as.data.frame(sam))
newdata <- as.data.frame(sam[sample(1:dim(sam)[1],90),])

# return the predicted value for each term #
predict(res_poi,newdata,type='terms',se.fit = TRUE)
predict(res_poi2,newdata)

data(eg1pop_bin_rho00)
formula <- Y~z1+z2+z3+b(x1,x2,V=eg1_V2,Tr=eg1_T2,d=2,r=1)
n <- 100
Npop <- nrow(eg1pop_bin_rho00)
# ind.pop <- (1:Npop)[!is.na(eg1pop_bin_rho00[,1])]
ind.pop <- (1:Npop)
sam.ind <- sort(sample(ind.pop,n,replace=FALSE))
sam <- eg1pop_bin_rho00[sam.ind,]
res_eg1_bin <- plbpsm(formula=formula,data=as.data.frame(sam),family=binomial())
eg1pop_bin_rho00=as.data.frame(eg1pop_bin_rho00)
newdata <- eg1pop_bin_rho00[1000:1100,]
predict(res_eg1_bin,as.data.frame(newdata),type='terms',se.fit = TRUE)

### Generalized Geoadditive Models with Model Identification ###
data(eg1pop_poi2)
formula <- Y~u(z1)+u(z2)+u(z3)+b(x1,x2,V=eg1_V2,Tr=eg1_T2,d=2,r=1)
n <- 100
Npop <- nrow(eg1pop_poi2)
# ind.pop <- (1:Npop)[!is.na(eg1pop_poi2[,1])]
ind.pop <- (1:Npop)
sam.ind <- sort(sample(ind.pop,n,replace=FALSE))
sam <- eg1pop_poi2[sam.ind,]
res_eg1_poi_add <- plbpsm(formula=formula,data=as.data.frame(sam),family='poisson')
newdata <- as.data.frame(sam[sample(1:dim(sam)[1],90),])

## backfitting = FALSE ##
pred_noBKF <- predict(res_eg1_poi_add,newdata,backfitting=FALSE)

## defualt backfitting for \code{res_eg1_poi_add}
pred_BKF <- predict(res_eg1_poi_add,newdata)

funstatpackages/GgAM documentation built on Nov. 4, 2019, 12:59 p.m.