spaMM (version 4.5.0)

fix_predVar: Prediction from models with nearly-singular covariance matrices

Description

This explains how to handle a warning occurring in computation of prediction variance, where the user is directed here.

For Matern or Cauchy correlation models with vanishing scale factor for distances, a warning may be produced when predict.HLfit (or get_predVar, etc.) is called with non-NULL newdata, because a nearly-singular correlation matrix of the random effect is met. To decide what to do in that case, users should compare the values of get_predVar(.) and get_predVar(., newdata=myfit$data) (see Example below). In the absence of numerical inaccuracies, The two values should be identical, and in the presence of such inaccuracies, the more reliable value is the first one. In really poor cases, the second syntax may yield negative prediction variances. If users deem the inaccuracies too large, they should use control=list(fix_predVar=TRUE) in the next call to predict.HLfit (or get_predVar, etc.) as shown in the Example. The drawback of this control is that the computation may be slower, and might even exceed memory capacity for large problems (some matrix operations being performed with exact rational arithmetic, which is memory-consuming for large matrices). it is also still experimental, in the sense that I fear that bugs (stop) may occur. If the user instead chooses control=list(fix_predVar=FALSE), the default standard floating-point arithmetic is used, but no warning is issued.

For fix_predVar left NULL (the default), standard floating-point arithmetic is also used. But in addition (with exceptions: see Details), the warning keeps being issued, and the (possibly costly) computation of the inverse of the correlation matrix is not stored in the fitted model object, hence is repeated for each new prediction variance computation. This is useful to remind users that something needs to be done, but for programming purposes where repeated warnings may be a nuisance, one can use control=list(fix_predVar=NA) which will issue a warning then perform as control=list(fix_predVar=FALSE), i.e. store an approximate inverse so the warning is not issued again. Finally, control=list(fix_predVar=NaN) will remove the inverse of the correlation matrix from the fitted model object, and start afresh as if the control was NULL.

Details

Nearly-singular correlation matrices of random effects occur in several contexts. For random-slope models, it commonly occurs that the fitted correlation between the random effects for Intercept and slope is 1 or -1, in which case the correlation matrix between these random effects is singular. This led to quite inaccurate computations of prediction variances in spaMM prior to version 3.1.0, but this problem has been fixed.

control=list(fix_predVar=NaN) may be more appropriate than control=list(fix_predVar=NULL) when predict.HLfit is called through code that one cannot control. For this reason, spaMM provides another mode of control of the default. It will convert control=list(fix_predVar=NULL) to other values when the call stack has call names matching the patterns given by
spaMM.getOption("fix_predVar") (as understood by grep). Thus if spaMM.getOption("fix_predVar")$"NA"=="MSL|bboptim", the default behaviour is that defined by control=list(fix_predVar=NA) when predict.HLfit is called through Infusion::MSL or blackbox::bboptim. FALSE or TRUE are handled in a similar way.

Examples

data("blackcap")
fitobject <- corrHLfit(migStatus ~ 1 + Matern(1|longitude+latitude),data=blackcap,
                       ranFix=list(nu=10,rho=0.001)) ## numerically singular C
get_predVar(fitobject,newdata=blackcap[6,]) 
## => warning => let us apply the recommended procedure:
get_predVar(fitobject) 
get_predVar(fitobject,newdata=fitobject$data) 
# Negative values again in the second case => easy decision:
get_predVar(fitobject,newdata=blackcap[1:6,], 
            control=list(fix_predVar=TRUE)) # now it's accurate
            # and the accuracy control is stored in the object:
get_predVar(fitobject,newdata=blackcap[1:6,]) 
# Clean and start afresh:
get_predVar(fitobject,newdata=blackcap[1:6,], 
            control=list(fix_predVar=NaN))