Prediction from a model fit
The following functions can be used to compute point predictions and/or various measures of uncertainty associated to such predictions. predict
can be used for prediction of the response variable by its expected value obtained as (the inverse link transformation of) the linear predictor (η) and more generally for terms of the form X_nβ+Z_nLv, for new design matrices X_n and Z_n. Various components of prediction variances and predictions intervals can also be computed using predict
. The get_
... functions are convenient extractors for such components. get_predCov_var_fix
extracts a block of a prediction covariance matrix. It was conceived for the specific purpose of computing the spatial prediction covariances between two “new” sets of geographic locations, without computing the full covariance matrix for both the new locations and the original (fitted) locations. When one of the two sets of new locations is fixed while the other varies, some expensive computations can be performed once for all sets of new locations, and be provided as the fix_X_ZAC.object
argument. The preprocess_fix_corr
extractor is designed to compute this argument.
## S3 method for class 'HLfit' predict(object, newdata = newX, newX = NULL, re.form = NULL, variances=list(), binding = FALSE, intervals = NULL, level = 0.95, blockSize = 2000L, type = "response", verbose=c(showpbar=eval(spaMM.getOption("barstyle"))), control=list(), ...) get_predCov_var_fix(object, newdata = NULL, fix_X_ZAC.object, fixdata, re.form = NULL, variances=list(disp=TRUE,residVar=FALSE,cov=FALSE), control=list(), ...) preprocess_fix_corr(object, fixdata, re.form = NULL, variances=list(residVar=FALSE, cov=FALSE), control=list()) get_fixefVar(...) get_predVar(..., variances=list(), which="predVar") get_residVar(...) get_respVar(...) get_intervals(..., intervals="predVar")
object |
The return object of fitting functions |
newdata |
Either NULL, a matrix or data frame, or a numeric vector. If If |
newX |
equivalent to newdata, available for back-compatibility |
re.form |
formula for random effects to include. By default, it is NULL, in which case all random effects are included. If it is NA, no random effect is included. If it is a formula, only the random effects it contains are retained. The other variance components are removed from both point prediction and |
variances |
A list whose elements control the computation of different estimated variances.
In particular, |
intervals |
NULL or character string or vector of strings. Provides prediction intervals with nominal level |
which |
any of |
level |
Coverage of the intervals. |
binding |
If |
fixdata |
A data frame describing reference data whose covariances with variable |
fix_X_ZAC.object |
The return value of calling |
blockSize |
Mainly for development purposes. For original or new data with many rows, it may be more efficient to split these data in small blocks, and this gives the maximum number or rows of the blocks. However, this will be ignored if a prediction covariance matrix is requested. |
type |
character string; The returned point prediction is on the response scale if |
control |
A list; a warning will direct you to relevant usage when needed. |
verbose |
A vector of booleans; it single currently used element is |
... |
further arguments passed to or from other methods. For the |
See the predVar
help page for information about the different concepts of prediction variances handled by spaMM (uncertainty of point prediction vs. of response) and about options controlling their computation.
If newdata
is NULL, predict
returns the fitted responses, including random effects, from the object.
Otherwise it computes new predictions including random effects as far as possible.
For spatial random effects it constructs a correlation matrix C between new locations and locations in the original fit. Then it infers the random effects in the new locations as C (L')^{-1} v (see spaMM
for notation). For non-spatial random effects, it checks whether any group (i.e., level of a random effect) in the new data was represented in the original data, and it adds the inferred random effect for this group to the prediction for individuals in this group.
In the point prediction of the linear predictor, the unconditional expected value of u is assigned to the realizations of u for unobserved levels of non-spatial random effects (it is zero in GLMMs but not for non-gaussian random effects), and the inferred value of u is assigned in all other cases. Corresponding values of v are then deduced. This computation yields the classical “BLUP” or empirical Bayes predictor in LMMs, but otherwise it may yield less well characterized predictors, where “unconditional” v may not be its expected value when the rand.family
link is not identity.
Intervals computations use the relevant variance estimates plugged in a Gaussian approximation, except for the simple linear model where it uses Student's t distribution.
See Details in Tpoisson
for questions specific to truncated distributions.
For predict
, a matrix or data frame (according to the binding
argument), with optional attributes frame
, intervals
, predVar
, fixefVar
, residVar
, and/or respVar
, the last four holding one or more variance vector or covariance matrices. The further attribute fittedName
contains the binding name, if any.
The get_
... extractor functions call predict
and extract from its result the attribute implied by the name of the extractor. By default, get_intervals
will return prediction intervals using predVar
.
get_predVar
with non-default which
argument has the same effect as the get_
... function whose name is implied by which
.
predVar
for information specific to prediction variances sensu lato;
get_cPredVar
for a bootstrap-corrected version of get_predVar
;
residVar
for an alternative extractor for residual variances, more general than get_residVar
.
data("blackcap") fitobject <- fitme(migStatus ~ 1 + Matern(1|longitude+latitude),data=blackcap, fixed=list(nu=4,rho=0.4,phi=0.05)) predict(fitobject) #### multiple controls of prediction variances ## (1) fit with an additional random effect grouped <- cbind(blackcap,grp=c(rep(1,7),rep(2,7))) fitobject <- fitme(migStatus ~ 1 + (1|grp) +Matern(1|longitude+latitude), data=grouped, fixed=list(nu=4,rho=0.4,phi=0.05)) ## (2) re.form usage to remove a random effect from point prediction and variances: predict(fitobject,re.form= ~ 1 + Matern(1|longitude+latitude)) ## (3) comparison of covariance matrices for two types of new data moregroups <- grouped[1:5,] rownames(moregroups) <- paste0("newloc",1:5) moregroups$grp <- rep(3,5) ## all new data belong to an unobserved third group cov1 <- get_predVar(fitobject,newdata=moregroups, variances=list(linPred=TRUE,cov=TRUE)) moregroups$grp <- 3:7 ## all new data belong to distinct unobserved groups cov2 <- get_predVar(fitobject,newdata=moregroups, variances=list(linPred=TRUE,cov=TRUE)) cov1-cov2 ## the expected off-diagonal covariance due to the common group in the first fit. ## Not run: #### Other extractors: fix_X_ZAC.object <- preprocess_fix_corr(fitobject,fixdata=blackcap) # ... for use in multiple calls to get_predCov_var_fix(): get_predCov_var_fix(fitobject,newdata=blackcap[14,],fix_X_ZAC.object=fix_X_ZAC.object) #### Prediction with distinct given phi's in different locations, # as specified by a resid.model: varphi <- cbind(blackcap,logphi=runif(14)) vphifit <- fitme(migStatus ~ 1 + Matern(1|longitude+latitude), resid.model = list(formula=~0+offset(logphi)), data=varphi, fixed=list(nu=4,rho=0.4)) # For respVar computation (i.e., response variance, often called prediction variance), # one then also needs to provide the variables used in 'resid.model', here 'logphi': get_respVar(vphifit,newdata=data.frame(latitude=1,longitude=1,logphi=1)) # For default 'predVar' computation (i.e., uncertainty in point prediction), # this is not needed: get_predVar(vphifit,newdata=data.frame(latitude=1,longitude=1)) #### point predictions and variances with new X and Z if(requireNamespace("rsae", quietly = TRUE)) { data("landsat", package = "rsae") fitobject <- fitme(HACorn ~ PixelsCorn + PixelsSoybeans + (1|CountyName), data=landsat[-33,]) newXandZ <- unique(data.frame(PixelsCorn=landsat$MeanPixelsCorn, PixelsSoybeans=landsat$MeanPixelsSoybeans, CountyName=landsat$CountyName)) predict(fitobject,newdata=newXandZ,variances = list(predVar=TRUE)) get_predVar(fitobject,newdata=newXandZ,variances = list(predVar=TRUE)) } ## End(Not run)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.