predict: Prediction from a model fit

Description

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 (\(\eta\)) and more generally for terms of the form X_n\(\beta\)+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.

Usage

# S3 method for HLfit
predict(object, newdata = newX, newX = NULL, re.form = NULL,
        variances=list(), binding = FALSE, intervals = NULL,
        level = 0.95, blockSize = 2000L, type = "response", verbose=NULL,  
        control=list(), na.action=na.omit, cluster_args=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")

Value

See Details in Tpoisson for questions specific to truncated distributions.

predict returns a matrix or data frame (according to the binding argument), with optional attributes including 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 frame is a data frame of predictors, and includes information about any na.action effect on the new data, except for predictions from multivariate response fits, where the frame is a list of data frames, each equivalent to the frame attribute for the univariate case. In the multivariate case, the additional mv attribute gives the predictions as a list with elements being the predictions for each submodel.

For multivariate response fits, the additional mv attribute gives the predictions as a list with elements for each submodel, and the frame is a list of data frames, each equivalent to the frame attribute for the univariate case.

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.

Arguments

object

The return object of fitting functions HLfit,corrHLfit,HLCor... returning an object inheriting from HLfit class.

newdata

Either NULL, a matrix or data frame, or a numeric vector.

If NULL, the original data are reused. Otherwise, all variables required to evaluate model formulas must be included. Which variables are required may depend on other arguments: see “prediction with given phi's” example, also illustrating the syntax when formulas include an offset.

If newdata is a numeric vector, its names (if any) are ignored. This makes it easier to use predict as an objective function for an optimization procedure such as optim, which calls the objective function on unnamed vectors. However, one must make sure that the order of elements in the vector is the order of first occurrence of the variables in the model formula. This order can be checked in the error message returned when calling predict on a newX vector of clearly wrong size, e.g. predict(<object>,newdata=numeric(0)).

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 calculations. If you want to retain only the spatial effects in the point prediction, but all variances, either use re.form and add missing variances (on linear predictor scale) manually, or ignore this argument and see Details and Examples for different ways of controlling variances.

variances

A list whose elements control the computation of different estimated variances. predict can return four components of prediction variance: fixefVar, predVar, residVar and respVar, whose definitions is detailed in predVar. They are all returned as attributes of the point predictions.

In particular, variances=list(predVar=TRUE) is suitable for uncertainty in point prediction, distinguished from the response variance given by list(respVar=TRUE). See the predVar help page for further explanations and other options.

intervals

NULL or character string or vector of strings. Provides prediction intervals with nominal level level, deduced from the given prediction variance term, e.g. intervals="predVar". Currently only intervals from fixefVar and predVar (and for LMMs respVar including the residual variance) may have a probabilistic meaning. Intervals returned in other cases are (currently) meaningless.

which

any of "predVar","respVar","residVar", "fixefVar", "intervals", or "naive".

level

Coverage of the intervals.

binding

If binding is a character string, the predicted values are bound with the newdata and the result is returned as a data frame (except for multivariate-response fits: see Details). The predicted values column name is the given binding, or a name based on it if the newdata already include a variable with this name. If binding is FALSE, The predicted values are returned as a one-column matrix and the data frame used for prediction is returned as an attribute (unless it was NULL). If binding is NA, a vector is returned, without the previous attributes.

fixdata

A data frame describing reference data whose covariances with variable newdata may be requested.

fix_X_ZAC.object

The return value of calling preprocess_fix_corr (see trivial Example). This is a more efficient way of providing information about the fixdata for repeated calls to get_predCov_var_fix with variable newdata.

blockSize

For data with many rows, it may be more efficient to perform some operations on slices of the data, and this gives the maximum number or rows of each slice. Further, parallelisation of computations over the slices is possible, as controlled by the cluster_args argument. Slicing and parallelisation may operate only if covariance matrices are not requested.

type

character string; The returned point predictions are on the response scale if type="response" (the default; for binomial response, a frequency 0<.<1). It is on the linear predictor scale if type="link".
* The “prediction variance” (as opposed to the response variance, see predVar) that may be returned as a "predVar" attribute of the point predictions is always on the linear predictor scale, even when type="response". If you want to extract this predVar transformed to the response scale, use predict(.,variances=list(respVar=TRUE)) and take the difference between the respVar and residVar attributes of the result.
* Prediction intervals (as opposed to the response intervals) will be on the linear predictor or response scale depending on type (new to versions more recent than 3.12.0).

control

A list; a warning will direct you to relevant usage when needed.

cluster_args

Passed to makeCluster. Parallel computations are possible if the slicing mechanism (as controlled by argument blockSize) is effective.

verbose

a list or a vector; The NULL default is interpreted as
c(showpbar=eval(spaMM.getOption("barstyle")), na=TRUE, na_once=FALSE).
showpbar controls whether to show a progress bar (and its style) in certain prediction variance computations; na controls whether to output a message when values required for some predictions are missing from the data; and setting na_once=TRUE is useful to prevent repetitive messages when predict is called repeatedly by an higher-level function (e.g., pdep_effects sets na_once=TRUE internally).

na.action

Possible values of this argument are some of the functions dealing with NAs in data frames (see na.omit). if it is set to na.exclude, NAs will be included in the returned point predictions, for rows of the newdata which do not provide information for all required predictor variables. The effect of the default na.omit is to not include such NAs (this differs from the default of, e.g., predict.lm). Implementation is limited; in particular, na.exclude currently does not have the effect of including NAs in the optional attributes providing (co-)variance information, except the "mv" attribute for predictions of multivariate-response fits.

...

further arguments passed to or from other methods. For the get_... functions, they are passed to predict.

Details

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.

There are cases where prediction without a newdata argument may give results of different length than prediction with newdata=<original data>, as for predict. Notably, for multivariate-response fits, different subsets of lines of the data may be used for each submodel depending on the availability of all variables (including the response variable) for each submodel, and the resulting fitted values from each submodel will be used from prediction; while prediction with newdata does not check the availability of a response variable.

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.

Binding the predicted values to a single data frame of predictors is not implemented for multivariate-response fits. It would not generally be feasible in this case without various side effects, because there is a data frame of predictors for each submodels, each with possibly different predictors. Instead, binding of predictions for each submodel to their respective frames is attempted. This is an experimental feature which may not work for all options. If this fails, the returned object when binding=FALSE can nevertheless be used to achieve suitable forms of binding: in particular, the frames are returned as the "frame" attribute (a list of data frames), and a list, over submodels, of the predictions is returned as the "mv" attribute.

Examples

Run this code

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))) 
fitobject2 <- 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(fitobject2,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(fitobject2,newdata=moregroups,
                     variances=list(linPred=TRUE,cov=TRUE))
moregroups$grp <- 3:7 ## all new data belong to distinct unobserved groups
cov2 <- get_predVar(fitobject2,newdata=moregroups,
                     variances=list(linPred=TRUE,cov=TRUE))
cov1-cov2 ## the expected off-diagonal covariance due to the common group in the first fit.

if (FALSE) {
#### 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))
}

}