Package 'spaMM'

Description

This is used in tests/test-adjacency-long.R

Usage

data("adjlg")

Format

Includes an adjacency matrix adjlgMat. and a data frame adjlg with 5474 observations on the following 8 variables.

ID

a factor with levels 1 to 1000

months

a numeric vector

GENDER

a character vector

AGE

a numeric vector

X1

a numeric vector

X2

a numeric vector

month

a numeric vector

BUY

a numeric vector

Source

The simulation code shown below is derived from an example produced by Jeroen van den Ochtend. Following a change incorporated in spaMM version 3.8.0, that implied stricter checks of the input matrix, it appeared that the precision matrix generated by this example had inappropriate (repeated) dimnames. This example was then updated to reproduce past fitting results with a correctly formatted matrix. Note that changing the names of an adjacency matrix (as below) is generally unwise as it generally changes the statistical model because these names are matched whenever possible to levels of the grouping factor in the data.

The code was also modified to compensate for changes in R's default random number generator.

Examples

data(adjlg)
## See further usage in tests/test-adjacency-long.R
## Not run: 
# as produced by:
  library(data.table) ## Included data produced using version 1.10.4.3
  library(igraph) ## Included data produced using version 1.2.1
  
  rsample <- function(N=100, ## size of implied adjacency matrix
                      month_max=10,seed) {
    if (is.integer(seed)) set.seed(seed)
    dt <- data.table(ID=factor(1:N))
    dt$months <- sample(1:month_max,N,replace=T) ## # of liens for each level of ID
    dt$GENDER <- sample(c("MALE","FEMALE"),N,replace=TRUE)
    dt$AGE <- sample(18:99,N,replace=T)
    dt$X1 <- sample(1000:9900,N,replace=T)
    dt$X2 <-  runif(N)
    
    dt <- dt[, c(.SD, month=data.table(seq(from=1, to=months, by = 1))), by = ID] 
    dt[,BUY := 0]
    dt[month.V1==months,BUY := sample(c(0,1),1),by=ID]
    setnames(dt,"month.V1","month")
    
    #### create adjacency matrix
    Network <- data.table(OUT=sample(dt$ID,N*month_max*4/10))
    Network$IN <- sample(dt$ID,N*month_max*4/10)
    Network <- Network[IN != OUT]
    Network <- unique(Network)
    g <- graph.data.frame(Network,directed=F)
    g <- add_vertices(g,sum(!unique(dt$ID) %in% V(g)),
             name=unique(dt[!dt$ID %in% V(g),list(ID)])) # => improper names
    Network <- as_adjacency_matrix(g,sparse = TRUE,type="both")
    colnames(Network) <- rownames(Network) <- seq(nrow(Network)) # post-v3.8.0 names 
    return(list(data=dt,adjMatrix=Network))
  }

  RNGkind("Mersenne-Twister", "Inversion", "Rounding"  )
  set.seed(123)
  adjlg_sam <- rsample(N=1000,seed=NULL) 
  RNGkind("Mersenne-Twister", "Inversion", "Rejection"  )
  #
  adjlg <- as.data.frame(adjlg_sam$data)
  adjlgMat <- adjlg_sam$adjMatrix

## End(Not run)

Description

get_any_IC computes model selection/information criteria such as AIC. See Details for more information about these criteria. The other extractors AIC and extractAIC are methods for HLfit objects of generic functions defined in other packages: AIC is equivalent to get_any_IC (for a single fitted-model object), and extractAIC returns the marginal AIC and the number of degrees of freedom for the fixed effects.

Usage

get_any_IC(object, nsim=0L, ..., verbose=interactive(),
           also_cAIC=TRUE, short.names=NULL)
## S3 method for class 'HLfit'
AIC(object, ..., nsim=0L, k, verbose=interactive(),
                    also_cAIC=TRUE, short.names=NULL)
## S3 method for class 'HLfit'
extractAIC(fit, scale, k, ..., verbose=FALSE)

Arguments

Details

The AIC is a measure (by Kullback-Leibler directed distance, up to an additive constant) of quality of prediction of new data by a fitted model. Comparing information criteria may be viewed as a fast alternative to a comparison of the predictive accuracy of different models by cross-validation. Further procedures for model choice may also be useful (e.g. Williams, 1970; Lewis et al. 2010).

The conditional AIC (Vaida and Blanchard 2005) applies the AIC concept to new realizations of a mixed model, conditional on the realized values of the random effects. Lee et al. (2006) and Ha et al (2007) defined a corrected AIC [i.e., AIC(D*) in their eq. 7] which is here interpreted as the conditional AIC.

Such Kullback-Leibler relative distances cannot generally be evaluated exactly and various estimates have been discussed. get_any_IC computes, optionally prints, and returns invisibly one or more of the following quantities:
* Akaike's classical AIC (marginal AIC, mAIC, i.e., minus twice the marginal log-likelihood plus twice the number of fitted parameters);
* a plug-in estimate (cAIC) and/or a bootstrap estimate (b_cAIC) of the conditional AIC;
* a focussed AIC for dispersion parameters (dispersion AIC, dAIC).

For the conditional AIC, Vaida and Blanchard's plug-in estimator involves the conditional likelihood, and degrees of freedom for (i) estimated residual error parameters and (ii) the overall linear predictor characterized by the Effective degrees of freedom already discussed by previous authors including Lee and Nelder (1996), which gave a plug-in estimator ($p_D$) for it in HGLMs. By default, the plug-in estimate of both the conditional AIC and of $n-p_D$ (GoFdf, where $n$ is the length of the response vector) are returned by get_any_IC. But these are biased estimates of conditional AIC and effective df, and an alternative procedure is available for GLM response families if a non-default positive nsim value is used. In that case, the conditional AIC is estimated by a bootstrap version of Saefken et al. (2014)'s equation 2.5; this involves refitting the model to each bootstrap samples, so it may take time, and a full cross-validation procedure might as well be considered for model selection.

The dispersion AIC has been defined from restricted likelihood by Ha et al (2007; eq.10). The present implementation will use restricted likelihood only if made available by an REML fit, otherwise marginal likelihood is used.

Value

get_any_IC, a numeric vector whose possible elements are described in the Details, and whose names are controlled by the short.names argument. Note that the bootstrap computation actually makes sense and works also for fixed-effect models (although it is not clear how useful it is in that case). The return value will still refer to its results as conditional AIC.

For AIC, If just one fit object is provided, the same return value as for get_any_IC. If multiple objects are provided, a data.frame built from such vectors, with rows corresponding to the objects.

For extractAIC, a numeric vector of length 2, with first and second elements giving

Likelihood is broadly defined up to a constant, which opens the way for inconsistency between different likelihood and AIC computations. In spaMM, likelihood is nothing else than the probability or probability density of the data as function of model parameters. No constant is ever added, in contrast to stats::extractAIC output, so there are discrepancies with the latter function (see Examples).

References

Ha, I. D., Lee, Y. and MacKenzie, G. (2007) Model selection for multi-component frailty models. Statistics in Medicine 26: 4790-4807.

Lee Y. and Nelder. J. A. 1996. Hierarchical generalized linear models (with discussion). J. R. Statist. Soc. B, 58: 619-678.

Lewis, F., Butler, A. and Gilbert, L. (2011), A unified approach to model selection using the likelihood ratio test. Methods in Ecology and Evolution, 2: 155-162. doi:10.1111/j.2041-210X.2010.00063.x

Saefken B., Kneib T., van Waveren C.-S., Greven S. (2014) A unifying approach to the estimation of the conditional Akaike information in generalized linear mixed models. Electron. J. Statist. 8, 201-225.

Vaida, F., and Blanchard, S. (2005) Conditional Akaike information for mixed-effects models. Biometrika 92, 351-370.

Williams D.A. (1970) Discrimination between regression models to determine the pattern of enzyme synthesis in synchronous cell cultures. Biometrics 26: 23-32.

Examples

data("wafers")
m1 <- fitme(y ~ X1+X2+X3+X1*X3+X2*X3+I(X2^2)+(1|batch), data=wafers, 
            family=Gamma(log))

get_any_IC(m1) 
# => The plug-in estimate is stored in the 'm1' object 
#    as a result of the previous computation, and is now returned even by: 
get_any_IC(m1, also_cAIC=FALSE)

if (spaMM.getOption("example_maxtime")>4) {
 get_any_IC(m1, nsim=100L, seed=123) # provides bootstrap estimate of cAIC.
 # (parallelisation options could be used, e.g. nb_cores=detectCores(logical=FALSE)-1L)
}

extractAIC(m1)

## Not run: 
# Checking (in)consistency with glm example from help("stats::extractAIC"):
utils::example(glm) # => provides 'glm.D93' fit object
logLik(glm.D93) # logL= -23.38066 (df=5)
dataf <- data.frame(counts=counts,outcome=outcome, treatment=treatment)
extractAIC(fitme(counts ~ outcome + treatment, family = poisson(), data=dataf))
# => 56.76132 = -2 logL + 2* df
extractAIC(glm.D93) # 56.76132 too
#
# But for LM:
lm.D93 <- lm(counts ~ outcome + treatment, data=dataf)
logLik(lm.D93) # logL=-22.78576 (df=6)
extractAIC(fitme(counts ~ outcome + treatment, data=dataf)) # 57.5715 = -2 logL + 2* df
extractAIC(lm.D93) # 30.03062

### Inconsistency also apparent in drop1 output for :
# Toy data from McCullagh & Nelder (1989, pp. 300-2), as in 'glm' doc:
clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100),
    lot1 = c(118,58,42,35,27,25,21,19,18),
    lot2 = c(69,35,26,21,18,16,13,12,12))
#    
drop1( fitme(lot1 ~ log(u), data = clotting), test = "F") # agains reports marginal AIC
# => this  may differ strongly from those returned by drop1( < glm() fit > ),
# but the latter are not even consistent with those from drop1( < lm() fit > )
# for linear models. Compare
drop1( lm(lot1 ~ log(u), data = clotting), test = "F") # consistent with drop1.HLfit()
drop1( glm(lot1 ~ log(u), data = clotting), test = "F") # inconsistent

## Discrepancies in drop1 output with Gamma() family:

gglm <- glm(lot1 ~ 1, data = clotting, family=Gamma())
logLik(gglm) # -40.34633 (df=2)

spgglm <-  fitme(lot1 ~ 1, data = clotting, family=Gamma())
logLik(spgglm)                      # -40.33777 (slight difference: 
#   see help("method") for difference in estimation method between glm() and fitme()).
# Yet this does not explain the following:

drop1(  fitme(lot1 ~ log(u), data = clotting, family=Gamma()), test = "F") 
# => second AIC is 84.676 as expected from above logLik(spgglm).
drop1( glm(lot1 ~ log(u), data = clotting, family=Gamma()), test = "F") 
# => second AIC is 1465.27, quite different from -2*logLik(gglm) + 2*df


## End(Not run)

Description

Autocorrelated gaussian random effects can be specified in terms of their covariance matrix, or in terms of the precision matrix (i.e. inverse covariance matrix). In a pre-processing step, spaMM may assess whether such precision matrices are sparse but the correlation matrix is dense, and if so, it may use “sparse-precision” algorithms efficient for this case. If the precision matrix does not appear sufficiently sparser than the correlation matrix, correlation matrices are used, and they can themselves be sparse or dense, with distinct algebraic methods used in each case.

For example, when the model includes a corrMatrix term specified by a covariance matrix, the precision matrix may be computed to assess its sparseness. The Example below illustrates a case where detecting sparsity of the precision matrix allows a faster fit. However, such a comparison of correlation and precision matrices takes time and is not performed for all types of random-effect structures. Instead, some fast heuristics may be used (see Details). The default selection of methods may not always be optimal, and may be overcome by using the control.HLfit argument of the fitting function (or by spaMM.options(), see Details). In particular one can use either control.HLfit=list(sparse_precision= <TRUE|FALSE>) or control.HLfit=list(algebra= <"spprec"|"spcorr"|"decorr">) with the obvious expected effects.

Such control may be useful when you already know that the precision matrix is sparse (as spaMM may even kindly remind you of, see Example below). In that case, it is also efficient to specify the precision matrix directly (see Example in Gryphon), as spaMM then assumes that sparse-precision methods are better without checking the correlation matrix.

Such control may also be useful when the correlation matrix is nearly singular so that computation of its inverse fails. This may occur if the model is poorly specified, but also occurs sometimes for valid correlation models because inversion of large matrices though Cholesky methods is not numerically accurate enough. In the latter case, you may be directed to this documentation by an error message, and specifying sparse_precision= FALSE may be useful.

Details

Currently the sparse-precision methods are selected by default in two cases (with possible exceptions indicated by specific messages): (1) for models including IMRF random effects; and (2) when the corrMatrix (or covStruct) syntax is used to provide a fixed precision matrix. Further, for models including autoregressive terms other than IMRF (i.e., adjacency, AR1), sparse-precision methods may or may not be selected on a simple heuristic based on the likely structure of the correlation matrix.

Algebraic methods can be controlled globally over all further fits by using
spaMM.options(sparse_precision= <TRUE|FALSE>)
and, among the correlation-based methods,
spaMM.options(QRmethod= <"sparse"|"dense">)
to select "spcorr" vs. "decorr" methods. Fit-specific controls (by control.HLfit) override these global ones.

See Also

pedigree

Examples

if (spaMM.getOption("example_maxtime")>6) {
 data("Gryphon")
 
 gry_df <- fitme(BWT ~ 1 + corrMatrix(1|ID), corrMatrix = Gryphon_A, 
                 data = Gryphon_df, method = "REML")
 how(gry_df)
               
 # => Note the message about 'Choosing matrix methods...'. 
 # Using control.HLfit=list(algebra="spprec") would indeed 
 # save the time used to select this method.
 
 # Conversely, using a correlation-based method would be a waste of time:
 
 gry_dn <- fitme(BWT ~ 1 + corrMatrix(1|ID), corrMatrix = Gryphon_A, 
               data = Gryphon_df, method = "REML",
               control.HLfit=list(sparse_precision=FALSE))
 how(gry_dn) # forced dense-correlation methods, which is slower here.
}

Description

For 948 “accessions” from European Arabidopsis thaliana populations, this data set merges the genotypic information at four single nucleotide polymorphisms (SNP) putatively involved in adaptation to climate (Fournier-Level et al, 2011, Table 1), with 13 climatic variables from Hancock et al. (2011).

Usage

data("arabidopsis")

Format

The data frame includes 948 observations on the following variables:

pos1046738, pos5510910, pos6235221, pos8132698

Genotypes at four SNP loci

LAT

latitude

LONG

longitude

seasonal, tempWarmest, tempColdest, preciWettest, preciDriest, preciCV, PAR_SPRING,

growingL, conseqCold, conseqFrFree, RelHumidSp, dayLSp, aridity

Thirteen climatic variables. See Hancock et al. (2011) for details about these variables.

Details

The response is binary so method="PQL/L" seems warranted (see Rousset and Ferdy, 2014).

Source

The data were retrieved from http://bergelson.uchicago.edu/regmap-data/climate-genome-scan on 22 February 2013 (they may no longer be available from there).

References

Fournier-Level A, Korte A., Cooper M. D., Nordborg M., Schmitt J., Wilczek AM (2011). A map of local adaptation in Arabidopsis thaliana. Science 334: 86-89.

Hancock, A. M., Brachi, B., Faure, N., Horton, M. W., Jarymowycz, L. B., Sperone, F. G., Toomajian, C., Roux, F., and Bergelson, J. 2011. Adaptation to climate across the Arabidopsis thaliana genome, Science 334: 83-86.

Rousset F., Ferdy, J.-B. (2014) Testing environmental and genetic effects in the presence of spatial autocorrelation. Ecography, 37: 781-790. doi:10.1111/ecog.00566

Examples

data("arabidopsis")
if (spaMM.getOption("example_maxtime")>2.5) {
  fitme(cbind(pos1046738,1-pos1046738)~seasonal+Matern(1|LAT+LONG),
        fixed=list(rho=0.119278,nu=0.236990,lambda=8.599),
        family=binomial(),method="PQL/L",data=arabidopsis)
}
## The above 'fixed' values are deduced from the following fit:
if (spaMM.getOption("example_maxtime")>46) {
  SNPfit <- fitme(cbind(pos1046738,1-pos1046738)~seasonal+Matern(1|LAT+LONG),
              verbose=c(TRACE=TRUE), 
              family=binomial(),method="PQL/L",data=arabidopsis)
  summary(SNPfit) # p_v=-125.0392
}

Description

These times-series correlation models can be declared as correlation models for random effect. The AR(p) model is here parametrized by the partial correlation coefficients of the levels of the random effect, {U_t}, corr(U_s,U_t|U_(s+1),...,U_(t-1)), with valid values in the hypercube ]-1,,1[^p (Barndorff-Nielsen and Schou, 1973). In the autoregressive-moving average ARMA(p,q) model, the AR part is parametrized in the same way. AR parameters are named "p1", "p2"..., and MA parameters are named "q1", "q2"... .

Implementation of the AR(p) model uses the sparsity of the inverse covariance matrix. In the ARMA(p,q) model, neither the covariance nor its inverse are sparse, so fits are expected to be more time- and memory-consuming.

Usage

# corrFamily constructors:
ARp(p=1L, fixed=NULL, corr=TRUE, tpar=1/(1+seq(p)))
ARMA(p=1L, q=1L, fixed=NULL, tpar=c(1/(1+seq_len(p)),1/(1+seq_len(q))))

Arguments

Value

The ARp and ARMA functions return a corrFamily descriptor, hence a list including element Cf, a function returning, for given ARMA or AR parameters, the correlation matrix for ARMA, or its inverse for ARp.

The fitted correlation matrix can be extracted from a fit object, as for any autocorrelated random effect, by Corr(<fit object>)[[<random-effect index>]].

References

Barndorff-Nielsen 0. and Schou G., 1973 On the parametrization of autoregressive models by partial autocorrelations. J. Multivariate Analysis 3: 408-419. doi:10.1016/0047-259X(73)90030-4

Examples

if (spaMM.getOption("example_maxtime")>2) {
ts <- data.frame(lh=lh,time=seq(48)) ## using 'lh' data from 'stats' package

## Default 'tpar' => AR1 model 
#
(ARpfit <-  fitme(lh ~ 1 + ARp(1|time), data=ts, method="REML"))
#
## which is equivalent to
#
(AR1fit <- fitme(lh ~ 1 +AR1(1|time), data=ts, method="REML"))

## AR(3) model 
#
(AR3fit <- fitme(lh ~ 1 + ARp(1|time, p=3), data=ts, method="REML"))

## Same but with fixed 2-lag partial autocorrelation 
#
(AR3fix <- fitme(lh ~ 1 + ARp(1|time, p=3, fixed=c(p2=0)), data=ts, method="REML"))
#      
# The fit should be statistically equivalent to
#
(AR3_fix <- fitme(lh ~ 1 + ARp(1|time, p=3), data=ts, method="REML",
                 fixed=list(corrPars=list("1"=c(p2=0)))))
#                 
# with subtle differences in the structure of the fit objects:
#
get_ranPars(AR3fix)$corrPars     # p2 was not a parameter of the model
get_ranPars(AR3_fix)$corrPars    # p2 was a fixed parameter of the model
#
# get_fittefPars() expectedly ignores 'p2' whichever way it was fixed.    


## Same as 'AR3fix' but with an additional MA(1) component
#
(ARMAfit <- fitme(lh ~ 1 + ARMA(1|time, p=3, q=1, fixed=c(p2=0)), 
                  data=ts, method="REML"))
}

Description

The lmerTest::contest function, drop1 and anova methods implement a number of tests for linear mixed models, e.g. using effective degrees of freeedom based on (a generalization of) Satterthwaite's method. These tests can be performed using spaMM fits through the conversion of the fit object, by the as_LMLT function, to an ad-hoc format acceptable as input to lmerTest's internal procedures. The separately documented drop1.HLfit and (optionally) anova.HLfit methods, when called on a single LMM fit object, perform the conversion by as_LMLT and call drop1 or anova methods defined by lmerTest.

Only the tests using lmerTest's default method ddf="Satterthwaite" are formally supported, as the converted object do not have the required format for the other methods. Only LMMs are handled by lmerTest, and residual-dispersion models are not yet handled by the conversion. However, the conversion extends lmerTest's functionality by handling all random-effect parameters handled by numInfo, therefore including (e.g.) spatial-correlation parameters not handled by lme4.

Usage

as_LMLT(fitobject, nuisance=NULL, verbose=TRUE, transf=TRUE, check_deriv=NULL, ...)

Arguments

Value

The value is returned invisibly. It is an S4 object of class "LMLT" with slots matching those required in objects of S4 class "lmerModLmerTest" when used by package lmerTest with ddf="Satterthwaite" (many additional slots of a formal "lmerModLmerTest" object are missing). The additional nuisance slot contains the nuisance list.

References

Alexandra Kuznetsova, Per B. Brockhoff and Rune H. B. Christensen (2017) lmerTest Package: Tests in Linear Mixed Effects Models. Journal of Statistical Software, 82(13), 1–26. doi:10.18637/jss.v082.i13

Examples

## Reproducing an example from the doc of lmerTest::contest.lmerModLmerTest,
#    using a spaMM fit as input.
## Not run: 
  data("sleepstudy", package="lme4")
  
  ## The fit:
  spfit <- fitme(Reaction ~ Days + I(Days^2) + (1|Subject) + (0+Days|Subject),
             sleepstudy, method="REML")
  
  ## Conversion:           
  spfit_lmlt <-  as_LMLT(spfit) 

  ## Functions from package lmerTest can then be called on this object:
  lmerTest::contest(spfit_lmlt, L=diag(3)[2:3, ])   # Test of 'Days + I(Days^2)'.
  #
  anova(spfit_lmlt, type="1")         # : using lmerTest:::anova.lmerModLmerTest()
  drop1(spfit_lmlt)                   # : using lmerTest:::drop1.lmerModLmerTest()


## End(Not run)

Description

Diverse autoregressive (AR) models are implemented in spaMM. This documentation describe the adjacency model (a conditional AR, i.e., CAR), and the AR1 model for time series. Other documentation deals with more or less distantly related models: ARp for more general AR(p) and ARMA(p,q) models for time series, and IMRF and MaternIMRFa for mesh-based approximations of geostatistical models.

An AR1 random effect is specified as AR1(1|<grouping factor>). It describes correlations between realizations of the random effect for (typically) successive time-steps by a correlation $\phi$, denoted ARphi in function calls. Nested AR1 effects can be specified by a nested grouping factor, as in AR1(1|<time index> %in% <nesting factor>).

A CAR random effect is specified as adjacency(1|<grouping factor>). The correlations among levels of the random effect form a matrix (I$-\rho$ adjMatrix$)^{-1}$, in terms of an adjMatrix matrix which must be provided, and of the scalar $\rho$, denoted rho in function calls. The rows and columns of adjMatrix must have names matching those of levels of the random effect or else are assumed to match a sequence, from 1 to the number of columns, of values of the geographic location index specifying the spatial random effect. For example, if the model formula is
y ~ adjacency(1|geo.loc) and <data>$geo.loc is 2,4,3,1,... the first row/column of the matrix refers to geo.loc=1, i.e. to the fourth row of the data.

Details

Efficient algorithms for CAR models have been widely discussed in particular in the econometric literature (e.g., LeSage and Pace 2009), but these models are not necessarily recommended for irregular lattices (see Wall, 2004 and Martellosio, 2012 for some insights on the implications of autoregressive models).

In CAR models, the covariance matrix of random effects u can be described as $\lambda$(I$-\rho$ W$)^{-1}$ where W is the (symmetric) adjacency matrix. HLCor uses the spectral decomposition of the adjacency matrix, written as W=VDV' where D is a diagonal matrix of eigenvalues $d_i$. The covariance of V'u is $\lambda$(I$-\rho$ D$)^{-1}$, which is a diagonal matrix with elements $\lambda_i$=$\lambda$/(1$-\rho d_i$). Hence $1/\lambda_i$ is in the linear predictor form $\alpha$+$\beta d_i$ This can be used to fit $\lambda$ and $\rho$ efficiently. A call to corrHLfit with the additional argument init.HLfit=list(rho=0) should be equivalent in speed and result to the HLCor call.

This is fast for small datasets (as in the example below) but more generic maximization algorithms may be preferable for large ones. It is suggested to use fitme generally unless one has a large number of small data sets to analyze. A call to fitme or corrHLfit without that initial value does not use the spectral decomposition. It performs numerical maximization of the likelihood (or restricted likelihood) as function of the correlation parameter $\rho$. The choice of fitting function may slightly impact the results. The ML fits by corrHLfit and HLCor should be practically equivalent. The REML fits should slightly differ from each other, due to the fact that the REML approximation for GLMMs does not maximize a single likelihood function.

If HLCor is used, the results are reported as the coefficients $\alpha$ ((Intercept)) and $\beta$ (adjd) of the predictor for $1/\lambda_i$, in addition to the resulting values of $\rho$ and of the common $\lambda$ factor.

Different fits may also differ in using or not algorithms that exploit the sparsity of the precision matrix of the autoregressive random effect. By default, spaMM tends to select sparse-precision algorithms for large datasets and large (i.e. many-level) random effects (details are complex). However, for AR1 models, the dimension of the implied precision matrix is determined by the extreme values of grouping factor (typically interpreted as a time index), as all intermediate values must be considered. Then, the correlation-based algorithms may be more efficient if only a few levels are present in the data, as only a small correlation matrix is required in that case.

References

LeSage, J., Pace, R.K. (2009) Introduction to Spatial Econometrics. Chapman & Hall/CRC.

Martellosio, F. (2012) The correlation structure of spatial autoregressions, Econometric Theory 28, 1373-1391.

Wall M.M. (2004) A close look at the spatial structure implied by the CAR and SAR models: Journal of Statistical Planning and Inference 121: 311-324.

Examples

##### AR1 random effect:
ts <- data.frame(lh=lh,time=seq(48)) ## using 'lh' data from stats package
fitme(lh ~ 1 +AR1(1|time), data=ts, method="REML")
# With fixed parameters:
# HLCor(lh ~ 1 +AR1(1|time), data=ts, ranPars=list(ARphi=0.5,lambda=0.25,phi=0.001))

##### CAR random effect:
data("scotlip")
# CAR by Laplace with 'outer' estimation of rho
if (spaMM.getOption("example_maxtime")>0.8) {          
  fitme(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
          adjMatrix=Nmatrix, family=poisson(), data=scotlip) 
}

# CAR by Laplace with 'inner' estimation of rho
HLCor(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
          adjMatrix=Nmatrix, family=poisson(), data=scotlip, method="ML")

Description

Returns a family object for beta-response models. The model described by such a family is characterized by a linear predictor, a link function, and the beta density for the residual variation.

The precision parameter prec of this family is a positive value such that the variance of the response given its mean $\mu$ is $\mu(1-\mu)/(1+$prec). prec is thus the precision parameter $\phi$ of Ferrari & Cribari-Neto (2004) and of the betareg package (Cribari-Neto & Zeileis 2010).

A fixed-effect residual-dispersion model can be fitted, using the resid.model argument, which is used to specify the form of the logarithm of the precision parameter (see Examples). Thus the variance of the response become $\mu(1-\mu)/(1+$exp(<specified linear expression>)).

Usage

beta_resp(prec = stop("beta_resp's 'prec' must be specified"), link = "logit")

Arguments

Details

Prior weights are meaningful for this family and handled as a factor of the precision parameter (as for GLM families) hence here not as a divisor of the variance (in contrast to GLM families): the variance of the response become $\mu(1-\mu)/(1+$prec*<prior weights>). However, this feature is experimental and may be removed in the future. The fitting function's resid.model argument may be preferred to obtain the same effect, by specifying an offset(log(<prior weights>)) in its formula (given the log link used in that model). As usual in spaMM, the offset(.) argument should be a vector and any variable necessary for evaluating it should be in the data.

Value

A list, formally of class c("LLF", "family"). See LL-family for details about the structure and usage of such objects.

References

Cribari-Neto, F., & Zeileis, A. (2010). Beta Regression in R. Journal of Statistical Software, 34(2), 1-24. doi:10.18637/jss.v034.i02

Ferrari SLP, Cribari-Neto F (2004). “Beta Regression for Modelling Rates and Proportions.” Journal of Applied Statistics, 31(7), 799-815.

See Also

Further examples in LL-family.

Examples

set.seed(123)
  beta_dat <- data.frame(y=runif(100),grp=sample(2,100,replace = TRUE), x_het=runif(100))
  
  fitme(y ~1+(1|grp), family=beta_resp(), data= beta_dat)
  ## same logL, halved 'prec' when prior weights=2 are used: 
  # fitme(y ~1+(1|grp), family=beta_resp(), data= beta_dat, prior.weights=rep(2,100))
  
  ## With model for residual dispersion:
  # fitme(y ~1+(1|grp), family=beta_resp(), data= beta_dat, resid.model= ~ x_het)

Description

Returns a family object for beta-binomial models. The model described by such a family is characterized by a linear predictor, a link function, and the beta-binomial distribution for residual variation.

The precision parameter prec of this family is a positive value such that the variance of the beta-distributed latent variable given its mean $\mu$ is $\mu(1-\mu)/(1+$prec). prec is thus the same precision parameter as for the beta response family (see beta_resp. The variance of the beta-binomial sample of size $n$ is response is $\mu(1-\mu)n(n+$prec$)/(1+$prec).

A fixed-effect residual-dispersion model can be fitted, using the resid.model argument, which is used to specify the form of the logarithm of the precision parameter (see Examples). Thus the variance of the latent beta-distributed variable becomes $\mu(1-\mu)/(1+$exp(<specified linear expression>)).

Usage

betabin(prec = stop("betabin's 'prec' must be specified"), link = "logit")

Arguments

Details

Prior weights are meaningful for this family and handled as a factor of the precision parameter of the latent beta-distributed variable: the variance of the latent variable become $\mu(1-\mu)/(1+$prec*<prior weights>). However, this feature is experimental and may be removed in the future. The fitting function's resid.model argument may be preferred to obtain the same effect, by specifying an offset(log(<prior weights>)) in its formula (given the log link used in that model). As usual in spaMM, the offset(.) argument should be a vector and any variable necessary for evaluating it should be in the data.

Value

A list, formally of class c("LLF", "family"). See LL-family for details about the structure and usage of such objects.

Examples

if (requireNamespace("agridat", quietly = TRUE)) {
  data("crowder.seeds", package = "agridat")
  fitme(cbind(germ,n-germ) ~ gen+extract+(1|plate), data=crowder.seeds, family=betabin())
} else {
  data(clinics)
  fitme(cbind(npos,nneg)~1+(1|clinic), family=betabin(), data=clinics)
}

Description

This data set is extracted from a study of genetic polymorphisms potentially associated to migration behaviour in the blackcap (Sylvia atricapilla). Across different populations in Europe and Africa, the average migration behaviour was found to correlate with average allele size (dependent on the number of repeats of a small DNA motif) at the locus ADCYAP1, encoding a neuropeptide. This data set is quite small and ill-suited for separating random-effect variance from residual variance. The likelihood surface for the Matérn model actually has local maxima.

Usage

data("blackcap")

Format

The data frame includes 14 observations on the following variables:

latitude

latitude, indeed.

longitude

longitude, indeed.

migStatus

migration status as determined by Mueller et al, from 0 (resident populations) to 2.5 (long-distance migratory populations)

means

Mean allele sizes in each population

pos

Numerical index for the populations

Details

Migration status was coded as : pure resident populations as '0', resident populations with some migratory restlessness as '0.5', partial migratory populations as '1', completely migratory populations migrating short-distances as '1.5', intermediate-distance migratory populations as '2' and distinct long-distance migratory populations as '2.5'.

Source

Data from Mueller et al. (2011), including supplementary material now available from doi:10.1098/rspb.2010.2567.

References

Mueller, J. C., Pulido, F., and Kempenaers, B. 2011. Identification of a gene associated with avian migratory behaviour, Proc. Roy. Soc. (Lond.) B 278, 2848-2856.

Examples

## see 'fitme', 'corrHLfit' and 'fixedLRT' for examples involving these data

Description

The Cauchy family of correlation functions is useful to describe spatial processes with power-law decrease of correlation at long distance. It is valid for Euclidean distances in spaces of any dimension, and for great-circle distances on spheres of any dimension. It has a scale parameter (rho, as in the Matérn correlation function), a shape (or “smoothness”, Gneiting 2013) parameter, and a long-memory dependence (or, more abstractly, “shape”; Gneiting 2013) parameter (Gneiting and Schlater 2004). The present implementation also accepts a Nugget parameter. The family can be invoked in two ways. First, the CauchyCorr function evaluates correlations, using distances as input. Second, a term of the form Cauchy(1|<...>) in a formula specifies a random effect with Cauchy correlation function, using coordinates found in a data frame as input. In the latter case, the correlations between realizations of the random effect for any two observations in the data will be the value of the Cauchy function at the scaled distance between coordinates specified in <...>, using “+” as separator (e.g., Cauchy(1|longitude+latitude)). A syntax of the form Cauchy(1|longitude+latitude %in% grp) can be used to specify a Cauchy random effect with independent realizations for each level of the grouping variable grp.

Usage

## Default S3 method:
CauchyCorr(d, rho=1, shape, longdep, Nugget=NULL)
# Cauchy(1|...)

Arguments

Details

The correlation at distance $d>0$ is

$(1-\textrm{Nugget}) (1+(\rho d)^\textrm{shape})^(-\textrm{longdep/shape})$

Value

Scalar/vector/matrix depending on input.

References

Gneiting, T. and Schlater M. (2004) Stochastic models that separate fractal dimension and the Hurst effect. SIAM Rev. 46: 269–282.

Gneiting T. (2013) Strictly and non-strictly positive definite functions on spheres. Bernoulli 19: 1327-1349.

Examples

data("blackcap")
fitme(migStatus ~ means+ Cauchy(1|longitude+latitude),
      data=blackcap,
      #  fixed=list(longdep=0.5,shape=0.5,rho=0.05)
     )
## The Cauchy family can be used in Euclidean spaces of any dimension:
set.seed(123)
randpts <- matrix(rnorm(20),nrow=5)
distMatrix <- as.matrix(proxy::dist(randpts))
CauchyCorr(distMatrix,rho=0.1,shape=1,longdep=10)

# See ?MaternCorr for examples of syntaxes for group-specific random effects,
#  also handled by Cauchy().

Description

A small data set used by Booth & Hobert (1998).

Usage

data("clinics")

Format

A data frame with 16 observations on the following 4 variables.

npos

a numeric vector

nneg

a numeric vector

treatment

a numeric vector

clinic

a numeric vector

References

Booth, J.G., Hobert, J.P. (1998) Standard errors of prediction in generalized linear mixed models. J. Am. Stat. Assoc. 93: 262-272.

Examples

data(clinics)
## Not run: 
# The dataset was built as follows
npos <- c(11,16,14,2,6,1,1,4,10,22,7,1,0,0,1,6)
ntot <- c(36,20,19,16,17,11,5,6,37,32,19,17,12,10,9,7)
treatment <- c(rep(1,8),rep(0,8))
clinic <-c(seq(8),seq(8))
clinics <- data.frame(npos=npos,nneg=ntot-npos,treatment=treatment,clinic=clinic)

## End(Not run)

Description

The COM-Poisson family is a generalization of the Poisson family which can describe over-dispersed as well as under-dispersed count data. It is indexed by a parameter nu that quantifies such dispersion. For nu>1, the distribution is under-dispersed relative to the Poisson distribution with same mean. It includes the Poisson, geometric and Bernoulli as special (or limit) cases (see Details). The COM-Poisson family is here implemented as a family object, so that it can be fitted by glm, and further used to model conditional responses in mixed models fitted by this package's functions (see Examples). nu is distinct from the dispersion parameter $\nu=1/\phi$ considered elsewhere in this package and in the GLM literature, as $\nu$ affects in a more specific way the log-likelihood.

Several links are now allowed for this family, corresponding to different versions of the COMPoisson described in the literature (e.g., Sellers & Shmueli 2010; Huang 2017).

Usage

COMPoisson(nu =  stop("COMPoisson's 'nu' must be specified"), 
           link = "loglambda")

Arguments

Details

The $i$th term of the distribution can be written $q_i/Z$ where $q_i=\lambda^i / (i!)^\nu$ and $Z=\sum_{(i=0)}^\infty q_i$, for $\lambda=\lambda(\mu)$ implied by its inverse relationship, the expectation formula $\mu=\mu(\lambda)=\sum_{(i=0)}^\infty i q_i(\lambda)/Z$. The case nu=0 is the geometric distribution with parameter $\lambda$; nu=1 is the Poisson distribution with mean $\lambda$; and the limit as nu -> $\infty$ is the Bernoulli distribution with expectation $\lambda/(1+\lambda)$.

From this definition, this is an exponential family model with canonical parameters $log(\lambda)$ and $\nu$. When the linear predictor $\eta$ specifies $log(\lambda(\mu))$, the canonical link is used (e.g., Sellers & Shmueli 2010). It is here nicknamed "loglambda" and does not have a known expression in terms of elementary functions. To obtain $\mu$ as the link inverse of the linear predictor $\eta$, one then first computes $\lambda=e^\eta$ and then $\mu(\lambda)$ by the expectation formula. For other links (Huang 2017), one directly computes $\mu$ by the link inverse (e.g., $\mu=e^\eta$ for link "log"), and then one may solve for $\lambda= \lambda(\mu)$ to obtain other features of the distribution.

The relationships between $\lambda$ and $\mu$ or other moments of the distribution involve infinite summations. These sums can be easily approximated by a finite number of terms for large nu but not when nu approaches zero. For this reason, the code may fail to fit distributions with nu approaching 0 (strong residual over-dispersion). The case nu=0 (the geometric distribution) is fitted by an ad hoc algorithm devoid of such problems. Otherwise, spaMM truncates the sum, and uses numerical integrals to approximate missing terms (which slows down the fitting operation). In addition, it applies an ad hoc continuity correction to ensure continuity of the result in nu=1 (Poisson case). These corrections affect numerical results for the case of residual overdispersion but are negligible for the case of residual underdispersion. Alternatively, spaMM uses Gaunt et al.'s (2017) approximations when the condition defined by spaMM.getOption("CMP_asympto_cond") is satisfied. All approximations reduces the accuracy of computations, in a way that can impede the extended Levenberg-Marquardt algorithm sometimes needed by spaMM.

The name COMP_nu should be used to set initial values or bounds on nu in control arguments of the fitting functions (e.g., fitme(.,init=list(COMP_nu=1))). Fixed values should be set by the family argument (COMPoisson(nu=.)).

Value

A family object.

References

Gaunt, Robert E. and Iyengar, Satish and Olde Daalhuis, Adri B. and Simsek, Burcin. (2017) An asymptotic expansion for the normalizing constant of the Conway–Maxwell–Poisson distribution. Ann Inst Stat Math doi:10.1007/s10463-017-0629-6.

Huang, Alan (2017) Mean-parametrized Conway-Maxwell-Poisson regression models for dispersed counts. Stat. Modelling doi:10.1177/1471082X17697749

G. Shmueli, T. P. Minka, J. B. Kadane, S. Borle and P. Boatwright (2005) A useful distribution for fitting discrete data: revival of the Conway-Maxwell-Poisson distribution. Appl. Statist. 54: 127-142.

Sellers KF, Shmueli G (2010) A Flexible Regression Model for Count Data. Ann. Appl. Stat. 4: 943–961

Examples

if (spaMM.getOption("example_maxtime")>0.9) {
  # Fitting COMPoisson model with estimated nu parameter:
  #
  data("freight") ## example from Sellers & Shmueli, Ann. Appl. Stat. 4: 943–961 (2010)
  fitme(broken ~ transfers, data=freight, family = COMPoisson())
  fitme(broken ~ transfers, data=freight, family = COMPoisson(link="log"))

  # glm(), HLCor() and HLfit() handle spaMM::COMPoisson() with fixed overdispersion:
  #
  glm(broken ~ transfers, data=freight, family = COMPoisson(nu=10))
  HLfit(broken ~ transfers+(1|id), data=freight, family = COMPoisson(nu=10),method="ML")
  
  # Equivalence of poisson() and COMPoisson(nu=1):
  #
  COMPglm <- glm(broken ~ transfers, data=freight, family = poisson())
  coef(COMPglm)
  logLik(COMPglm)
  COMPglm <- glm(broken ~ transfers, data=freight, family = COMPoisson(nu=1))
  coef(COMPglm)
  logLik(COMPglm)
  HLfit(broken ~ transfers, data=freight, family = COMPoisson(nu=1))

}

Description

An example of a composite random effect is corrMatrix(sex|pair). It combines features of a random-coefficient model (sex|pair) and of an autocorrelated random effect corrMatrix(.|.). The random-coefficient model is characterized by a $2*2$ covariance matrix C for the random effects $u_{1,pair}$ and $u_{2,pair}$ both affecting each of the two sexes for each pair, and the corrMatrix random effect assumes that elements of each of the two vectors $u_i=(u_{i,pair})$ for pair=1,...,$P$ are correlated according to a given $P*P$ correlation matrix A. Then the composite random effect is defined as the one with 2$P*$2$P$ covariance matrix kronecker(C,A).

The definition of composite random effects through the kronecker product may be motivated and understood in light of a quantitative-genetics application (see help("Gryphon") for an example). In this context the two response variables are typically two individual traits. Each trait is affected by two sets of genes, the effect of each set being represented by a gaussian random effect (u_1 or u_2). The effect of genetic relatedness on the correlation of random effects u_i,ID among individuals ID within each set $i$ of genes is described by the corrMatrix A. The effects on the two traits for each individual are interpreted as different linear combinations of these two random effects (the coefficients of these linear combinations determining the C matrix). Under these assumptions the correlation matrix of the responses (in order (trait, individual)=(1,1)...(1,ID)... (2,1)...(2,ID)...) is indeed kronecker(C,A).

Composite random effects are not restricted to corrMatrix terms, and can also be fitted for multivariate-response models. For example, Matern(mv(1,2)|longitude+latitude) terms can be fitted, in which case the correlation model is still defined through the Kronecker product, where A will be a (fitted) Matérn correlation matrix, and C will be the correlation matrix of the fitted random-coefficient model for the mv virtual factor for multivariate response.

Implementation and tests of composite random effects is work in progress, with the following ones having been tested: corrMatrix, Matern, Cauchy, adjacency, IMRF, AR1, and to a lesser extent MaternIMRFa. Fits of other composite terms may fail. Even if they succeed, not all post-fit procedures may be operational: in particular, prediction (and then, simulation) with newdata may not yet work. Further, as for random-coefficient terms in univariate-response models, some components of the computed prediction variance depend on a poorly characterized approximation, yielding different results for different fitting algorithms (see Details in predVar).

The summary of the model provides fitted parameters for A if this matrix derives from a parametric correlation model (e.g., Matern), and a description of the C matrix where it is viewed as a covariance matrix, parametrized by its variances and its correlation coefficient(s). In a standard random-coefficient model these variances would be those of the correlated random effects (see summary.HLfit). In the composite random-effect model this is not necessarily so as the variance of the correlated random effects also depend on the variances implied by the A matrix, which are not necessarily 1 if A is a covariance matrix rather than simply a correlation matrix.

A <prefix>(<LHS>|<RHS>) term is *not* a composite random effect when the LHS in boolean, a factor from boolean, or “0+numeric”. See the Matérn examples, and the corrMatrix “<LHS> is 0+numeric” example, below.

Examples

if (spaMM.getOption("example_maxtime")>1.8) {

########################
#### corrMatrix examples
########################

## Toy data preparation

data("blackcap")
toy <- blackcap
toy$ID <- gl(7,2)
grp <- rep(1:2,7)
toy$migStatus <- toy$migStatus +(grp==2)
toy$loc <- rownames(toy) # to use as levels matching the corrMatrix dimnames

toy$grp <- factor(grp)
toy$bool <- toy$grp==1L
toy$boolfac <- factor(toy$bool)
toy$num <- seq(from=1, to=2, length.out=14)

## Build a toy corrMatrix as perturbation of identity matrix:
n_rhs <- 14L
eps <- 0.1
set.seed(123)
rcov <- ((1-eps)*diag(n_rhs)+eps*rWishart(1,n_rhs,diag(n_rhs)/n_rhs)[,,1])
# eigen(rcov)$values
colnames(rcov) <- rownames(rcov) <- toy$loc # DON'T FORGET NAMES


##### Illustrating the different LHS types

### <LHS> is logical (TRUE/FALSE) => No induced random-coefficient C matrix; 
#   corrMatrix affects only responses for which <LHS> is TRUE:
#
(fit1 <- fitme(migStatus ~ bool + corrMatrix(bool|loc), data=toy, corrMatrix=rcov))
#
# Matrix::image(get_ZALMatrix(fit1))


### <RHS> is a factor built from a logical => same a 'logical' case above:
#
(fit2 <- fitme(migStatus ~ boolfac + corrMatrix(boolfac|loc), data=toy, corrMatrix=rcov))
#
# Matrix::image(get_ZALMatrix(fit2))


### <RHS> is a factor not built from a logical: 
# (grp|.) and (0+grp|.) lead to equivalent fits of the same composite model, 
#   but contrasts are not used in the second case and the C matrices differ,
#   as for standard random-coefficient models.
#
(fit1 <- fitme(migStatus ~ grp +  corrMatrix(grp|loc), data=toy, corrMatrix=rcov))
(fit2 <- fitme(migStatus ~ grp +  corrMatrix(0+grp|loc), data=toy, corrMatrix=rcov))
# 
# => same fits, but different internal structures:
Matrix::image(fit1$ZAlist[[1]]) # (contrasts used) 
Matrix::image(fit2$ZAlist[[1]]) # (contrasts not used)
# Also compare ranef(fit1) versus ranef(fit2) 
#
#
## One can fix the C matrix, as for standard random-coefficient terms 
#
(fit1 <- fitme(migStatus ~ grp +  corrMatrix(0+grp|loc),data=toy, corrMatrix=rcov, 
               fixed=list(ranCoefs=list("1"=c(1,0.5,1)))))
#       
# same result without contrasts hence different 'ranCoefs':             
#
(fit2 <- fitme(migStatus ~ grp +  corrMatrix(grp|loc), data=toy, corrMatrix=rcov, 
               fixed=list(ranCoefs=list("1"=c(1,-0.5,1)))))


### <LHS> is numeric (but not '0+numeric'):
# composite model with C being 2*2 for Intercept and numeric variable
#
(fitme(migStatus ~ num +  corrMatrix(num|loc), data=toy, corrMatrix=rcov))

### <LHS> is 0+numeric: no random-coefficient C matrix 
#  as the Intercept is removed, but the correlated random effects 
#  arising from the corrMatrix are multiplied by sqrt(<numeric variable>)
#
(fitme(migStatus ~ num +  corrMatrix(0+num|loc), data=toy, corrMatrix=rcov))


### <LHS> for multivariate response (see help("Gryphon") for more typical example)
## More toy data preparation for multivariate response
ch <- chol(rcov)
set.seed(123)
v1 <- tcrossprod(ch,t(rnorm(14,sd=1)))
v2 <- tcrossprod(ch,t(rnorm(14,sd=1)))
toy$status <- 2*v1+v2
toy$status2 <- 2*v1-v2

## Fit:
fitmv(submodels=list(mod1=list(status ~ 1+ corrMatrix(0+mv(1,2)|loc)),
                     mod2=list(status2 ~ 1+ corrMatrix(0+mv(1,2)|loc))), 
      data=toy, corrMatrix=rcov)
      
##################################################
#### Matern examples: sex-dependent random effects
##################################################

if (spaMM.getOption("example_maxtime")>2) {
data(Leuca)
subLeuca <- Leuca[c(1:10,79:88),] # subset for faster examples

# The random effects in the following examples are composite because 'sex' is not 
# boolean nor factor from boolean. If 'Matern(factor(female)|x+y)' were used, the effect 
# would be the same 'Matern(female|x)', fitting 

  fitme(fec_div ~ 1 + Matern(sex|x+y),data=subLeuca)  # => builds a random effect
# correlated across sexes, from 2 independent realizations u_1 and u_2 of 20 values 
# (for the 20 locations in the data). In a (sex|x) term the 20 values would be 
# independent from each other within each u_i. In the Matern(sex|x+y) such 20 values 
# are autocorrelated within each u_i.

# For pedagogic purposes, one can also fit
  fitme(fec_div ~ 1 + Matern(sex|x + y %in% sex),data=subLeuca)
# which again builds a random effect from 2 independent realizations 
# u_1 and u_2, but each u_i now contains two realizations u_i1 and u_i2 of 10 values,
# autocorrelated within each u_ij following the Matern model, 
# but independent between u_i1 and u_i2. As a result, the overall random effect in each sex, 
# v_m or v_f, is a weighted sum of two sex-specific Matern random effect, 
# so that v_m and v_f are independent from each other. 
}

}

Description

This function interfaces two procedures: a profile confidence interval procedure implemented for fixed-effects coefficients only; and a parametric bootstrap procedure that can be used to provide confidence interval for any parameter, whether a canonical parameter of the model or any function of one or several such parameters. The bootstrap is performed if the parm argument is a function or a quoted expression or if the boot_args argument is a list. The profile confidence interval is computed if neither of these conditions is true. In that case parm must be the name(s) of some fixed-effect coefficient, and the (p_v approximation of the) profile likelihood ratio for the given parameter is used to define the interval, where the profiling is over all other fitted parameters, including other fixed-effects coefficients, as well as variances of random effects and spatial correlations if these were fitted.

Of related interest, see numInfo which evaluates numerically the information matrix for given sets of canonical model parameters, from which asymptotic confidence intervals can be deduced.

Usage

## S3 method for class 'HLfit'
confint(object, parm, level=0.95, verbose=TRUE, 
                          boot_args=NULL, format="default", ...)

Arguments

Value

The format of the value varies, but in all cases distinguished below, one or more tables are included, as a table attribute, in the format of the stats::confint result, to facilitate consistent extraction of results. By default confint returns invisibly the full values described below, but if format="stats", only the table attribute is returned.

If a profile CI has been computed for a single parameter, a list is returned including the confidence interval as shown by verbose=TRUE, and the fits lowerfit and upperfit giving the profile fits at the confidence bounds. This list bears the table attribute.

If a profile CI has been computed for several parameters, a structured list, named according to the parameter names, of such single-parameter results is returned, and a single table attribute for all parameters is attached to the upper level of the list.

If a bootstrap was performed, for a single parameter the result of the boot.ci call is returned, to which a table attribute is added. This table is now a list of tables for the different bootstrap CI types (default being normal, percent, and basic), each such table in the format of the stats::confint results. For several parameters, a named list of boot.ci results is returned, its names being the parameter names, and the table attribute is attached to the upper level of the list.

The boot.ci return value for each parameter includes the call to boot.ci. This call is typically shown including a long t vector, which makes a bulky display. spaMM hacks the printing to abbreviate long ts.

See Also

numInfo for information matrix.

Examples

data("wafers")
wfit <- HLfit(y ~X1+(1|batch), family=Gamma(log), data=wafers, method="ML")
confint(wfit,"X1")  # profile CI
if (spaMM.getOption("example_maxtime")>30) {
   
   # bootstrap CI induced by 'boot_args':
   confint(wfit,names(fixef(wfit)), boot_args=list(nsim=99, seed=123)) 
   
   # bootstrap CI induced by 'parm' being a function:
   confint(wfit,parm=function(v) fixef(v), 
           boot_args=list(nb_cores=10, nsim=199, seed=123))
   
   # Same effect if 'parm' is a quoted expression in terms of 'hlfit':
   confint(wfit,parm=quote(fixef(hlfit)), 
           boot_args=list(nb_cores=10, nsim=199, seed=123))
           
   # CI for the variance of the random effect:          
   ci <- confint(wfit,parm=function(fit){get_ranPars(fit)$lambda[1]}, 
        boot_args=list(nb_cores=10, nsim=199, seed=123))
   # The distribution of bootstrap replicates:
   plot(ecdf(ci$call$t))
   # We may be far from ideal condition for accuracy of bootstrap intervals;
   # for variances, a log transformation may sometimes help, but not here.
   
   # Passing arguments to child processes, as documented in help("spaMM_boot"):
   set.seed(123)
   rvar <- runif(nrow(wafers))
   wrfit <- fitme(y ~X1+(1|batch), family=Gamma(log), data=wafers, fixed=list(phi=rvar))
   confint(wrfit, parm = "(Intercept)", boot_args = list(nsim = 100, nb_cores = 2,
           fit_env = list(rvar=rvar)))
   
}

Description

A list of parameters controlling the HLfit fitting algorithm (potentially called by all fitting functions in spaMM), which should mostly be ignored in routine use. Possible controls are:

algebra, sparse_precision: see algebra;

conv.threshold and spaMM_tol: spaMM_tol is a list of tolerance values, with elements Xtol_rel and Xtol_abs that define thresholds for relative and absolute changes in parameter values in iterative algorithms (used in tests of the form “d(param)< Xtol_rel * param + Xtol_abs”, so that Xtol_abs is operative only for small parameter values). conv.threshold is the older way to control Xtol_rel. Default values are given by spaMM.getOption("spaMM_tol");

break_conv_logL: a boolean specifying whether the iterative algorithm should terminate when log-likelihood appears to have converged (roughly, when its relative variation over on iteration is lower than 1e-8). Default is FALSE (convergence is then assessed on the parameter estimates rather than on log-likelihood).

iter.mean.dispFix: the number of iterations of the iterative algorithm for coefficients of the linear predictor, if no dispersion parameters are estimated by the iterative algorithm. Defaults to 200 except for Gamma(log)-family models;

iter.mean.dispVar: the number of iterations of the iterative algorithm for coefficients of the linear predictor, if some dispersion parameter(s) is estimated by the iterative algorithm. Defaults to 50 except for Gamma(log)-family models;

max.iter: the number of iterations of the iterative algorithm for joint estimation of dispersion parameters and of coefficients of the linear predictor. Defaults to 200. This is typically much more than necessary, unless there is little information to separately estimate $\lambda$ and $\phi$ parameters;

resid.family: was a previously documented control (before version 2.6.40), and will still operate as previously documented, but should not be used in new code.

Usage

# <fitting function>(., control.HLfit=list(...))

Description

spaMM fits can produce various convergence warnings or messages.

Messages referring to convergence issues in initialization can generally be ignored but may help to diagnose other apparent problems, if any.

Warnings referring to .calc_dispGammaGLM (for residual-dispersion fits) can generally be ignored when they show a small criterion (say <1e-5) but may otherwise suggest that the final fit did not optimize its objective.

Messages pointing to slow convergence and drawing users to this doscumentation do not necessarily mean the fit is incorrect. Rather, they suggest that another fitting strategy could be tried. Keep in mind that several parameters (notably the dispersion parameters: the variance of random effects and the residual variance parameter, if any) can be estimated either by the iterative algorithms, or by generic optimization methods. In my experience, slow convergence happens in certain cases where a large random-effect variance is considered by the algorithm used.

How to know which algorithm has been selected for each parameter? fitme(., verbose=c(TRACE=TRUE)) shows successive values of the variables estimated by optimization (See Examples; if no value appears, then all are estimated by iterative methods). The first lines of the summary of a fit object should tell which variances are estimated by the “outer” method.

If the iterative algorithm is being used, then it is worth trying to use the generic optimization methods. In particular, if you used HLfit, try using fitme; if you already use fitme, try to enforce generic optimization of the random-effect variance(s) (see inits). Conversely, if generic optimization is being used, the maximum lambda value could be controlled (say, upper=list(lambda=c(10,NA))), or the iterative algorithm can be called (see inits again).

For the largest datasets, it may be worth comparing the speed of the "spcorr" and "spprec" choices of the algebra control, in case spaMM has not selected the most appropriate by default. However, this will not be useful for geostatistical models with many spatial locations.

Examples

# See help("inits") for examples of control by initial values.

Description

This implements the Matérn correlation structure (see Matern) for use with lme or glmmPQL. Usage is as for others corSpatial objects such as corGaus or corExp, except that the Matérn family has an additional parameter. This function was defined for comparing results obtained with corrHLfit to those produced by lme and glmmmPQL. There are problems in fitting (G)LMMs in the latter way, so it is not a recommended practice.

Usage

corMatern(value = c(1, 0.5), form = ~1, nugget = FALSE, nuScaled = FALSE, 
          metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)

Arguments

Details

This function is a constructor for the corMatern class, representing a Matérn spatial correlation structure. See MaternCorr for details on the Matérn family.

Value

an object of class corMatern, also inheriting from class corSpatial, representing a Matérn spatial correlation structure.

Note

The R and C code for the methods for corMatern objects builds on code for corSpatial objects, by D.M. Bates, J.C. Pinheiro and S. DebRoy, in a circa-2012 version of nlme.

References

Mixed-Effects Models in S and S-PLUS, José C. Pinheiro and Douglas M. Bates, Statistics and Computing Series, Springer-Verlag, New York, NY, 2000.

See Also

glmmPQL, lme

Examples

## LMM
data("blackcap")
blackcapD <- cbind(blackcap,dummy=1) ## obscure, isn't it? 
## With method= 'ML' in lme, The correlated random effect is described 
##  as a correlated residual error and no extra residual variance is fitted:
nlme::lme(fixed = migStatus ~ means, data = blackcapD, random = ~ 1 | dummy, 
    correlation = corMatern(form = ~ longitude+latitude  | dummy), 
    method = "ML", control=nlme::lmeControl(sing.tol=1e-20))

## Binomial GLMM
if (spaMM.getOption("example_maxtime")>32) {
 data("Loaloa")
 LoaloaD <- cbind(Loaloa,dummy=1) 
 MASS::glmmPQL(fixed =cbind(npos,ntot-npos)~elev1+elev2+elev3+elev4+maxNDVI1+seNDVI, 
        data = LoaloaD, random = ~ 1 | dummy,family=binomial, 
        correlation = corMatern(form = ~ longitude+latitude | dummy))
}

Description

corr_family objects provide a convenient way to implement correlation models handled by spaMM, analogous to family objects. These objects are undocumented (but there are documentation pages for each of the models implemented).

Usage

# Matern(...)           # see help(Matern)
# Cauchy(...)           # see help(Cauchy)
# corrMatrix(...)       # see help(corrMatrix)
# AR1(...)              # see help(AR1)
# adjacency(...)        # see help(adjacency)
# IMRF(...)             # see help(IMRF)
## S3 method for class 'corr_family'
print(x,...)

Arguments

Description

One can declare and fit correlated random effects belonging to a user-defined correlation (or covariance) model (i.e., a parametric family of correlation matrices, although degenerate case with no parameter are also possible). This documentation is a first introduction to this feature. It is experimental in the sense that its design has been formalized only from a limited number of corrFamily examples, and that the documentation is not mature. Implementing prediction for random-effects defined in this way may be tricky. A distinct documentation corrFamily-design provides more information for the efficient design of new correlation models to be fitted in this way.

A simple example of random-effect model implemented in this way is the autoregressive model of order $p$ (AR(p) in the literature; specifically documented elsewhere, see ARp). It can be used as a formula term like other autocorrelated random-effects predefined in spaMM, to be fitted by fitme or fitmv:

fitme(lh ~ 1 + ARp(1|time, p=3),  # <= declaration of random effect
  < data and other possible arguments >)

User-defined correlation models should be registered for this simple syntax to work (see Details for an alternative syntax):

myARp <- ARp                   # 'myARP' is thus a user-defined model
register_cF("myARp")        # Register it so that the next call works
fitme(lh ~ 1 + myARp(1|time, p=3),  
  < data and other possible arguments >)

The ARp object here copied in myARp is a function (the corrFamily constructor) which returns a list (the corrFamily descriptor) which contains the necessary information to fit a random effect with an AR(p) correlation. The p argument in the myARp(1|time, p=3) term enforces evaluation of myARp(p=3), producing the descriptor for the AR(3) model. The structure of this descriptor is

List of 5
 $ Cf            :function (parvec)  
  ..- < with some attribute >
 $ tpar         : num [1:3] 0.5 0.333 0.25
 $ type         : chr "precision"
 $ initialize   :function (Zmatrix, ...) 
  ..- < with some attribute >
 $ fixed        : NULL
 $ calc_moreargs:function (corrfamily, ...)  
  ..- < with some attribute >
 $ levels_type        : chr "time_series"
 $ calc_corr_from_dist:function (ranFix, char_rd, distmat, ...)  
  ..- < with some attribute >
 < and possibly other elements > 

The meaning of these elements and some additional ones is explained below.

Only Cf and tpar are necessary elements of a corrFamily object. If one designs a new descriptor where some other elements are absent, spaMM will try to provide plausible defaults for these elements. Further, if the descriptor does not provide parameter names (as the names of tpar, or in some more cryptic way), default names "p1", "p2"... will be provided.

Usage

## corrFamily descriptor provided as a list of the form
#
# list(Cf=<.>, tpar=<.>, ...)

## corrFamily constructor: any function that returns 
#    a valid corrFamily descriptor
#
# function(tpar=<.>, fixed=<.>, ...) # typical but not mandatory arguments

## There is a distinct documentation page for 'register_cF'.

Arguments

Elements of the corrFamily descriptor:

fitting-function arguments:

lower, upper, init and fixed optimization controls can be used to control optimization of continuous parameters as for other random-effect parameters. They are specified as numeric vectors, themselves being element of the corrPars list (see Examples in corrFamily-design). Parameter names (consistent with those to be used for the tpar argument) may be required to disambiguate incomplete vectors (e.g., to specify only its second element). Apart from fixed ones, any of the values not specified through the fitting-function arguments will be sought in the return value of the calc_moreargs function, if provided in the descriptor. If the lower or upper information is missing there, it must be provided throught the fitting-function call. If the init information is missing, a default value will be deduced from the bounds. The init specification is thus always optional while the bounds specification is optional only if the descriptor provides default values.

Arguments of the corrFamily constructor

These may be ad libitum, as design rules are defined only for the returned descriptor. However, arguments tpar, fixed, and public of predefined constructors, such as ARp, are designed to match the above respective elements of the descriptor.

Details

* Constructor elements for prediction:

For prediction of autocorrelated random effects, one must first assess whether levels of the random effect not represented in the fit are possible in new data (corresponding to new spatial locations in geostatistical models, or new time steps in time series). In that case need_Cnn must be TRUE (Interpolated MRFs do not require this as all required random-effect levels are determined by the IMRF mesh argument rather than by the fitted data or new data).

Next, for autocorrelated random effects where need_Cnn is TRUE, a make_new_corr_lists function must be provided, except when a calc_corr_from_dist function is instead provided (which may be sufficient for models that construct the correlation from a spatial distance matrix). When need_Cnn is FALSE, a make_new_corr_lists function may still be needed.

The Examples section provides a simple example of such design, and the source code of the ARp or ARMA constructors provide further examples. They show that the make_new_corr_lists function may assign matrices or vectors as elements of several lists contained in a newLv_env environment. A matrix is assigned in the cov_newLv_oldv_list, specifying correlations between “new” levels of the random effect (implied by the new data) and “old” levels (those already included in the design matrix of the random effect for the fit). If need_Cnn is TRUE, a second matrix may be assigned in the cov_newLv_newLv_list, specifying correlation between “new” levels, and the diagonal of this matrix is assigned in the diag_cov_newLv_newLv_list. The overall structure of the code (the conditions where these assignments are performed, and the list indices), should be conserved.

When calling simulate(., newdata=<non-NULL>, type="marginal", a fourth matrix may be useful, assigned into a L_newLv_newLv_list, specifying the matrix root (as a tcrossprod factor) of the correlation matrix stored in cov_newLv_newLv_list. The relevant spaMM procedure will however try to compute it on the fly when it has not been provided by the make_new_corr_lists function.

* Fitting a corrFamily without a constructor:

It is possible to use an unregistered corrFamily, as follows:

AR3 <- ARp(p=3)          # Generate descriptor of AR model of order 3

fitme(lh ~ 1 + corrFamily(1|time),  # <= declaration of random effect
  covStruct=list( 
    corrFamily= AR3     # <= declaration of its correlation structure
  ), 
  < data and other possible arguments >)

Here the fit only uses a descriptor list, not a constructor function. This descriptor is here provided to the fitting function as an element of the covStruct argument (using the general syntax of this argument), and in the model formula the corresponding random effect is declared as a term of the form
corrFamily(1|<grouping factor>).

This syntax is more complex than the one using a registered constructor, but it might be useful for development purposes (one only has to code the descriptor, not the constructor function). However, it is not general; in particular, using registered constructors may be required to obtain correct results when fitting multivariate-response models by fitmv.

See Also

See ARp, diallel, and MaternIMRFa for basic examples of using a predefined corrFamily descriptor, and corrFamily-design for more geeky stuff including examples of implementing simple new correlation families.

Examples

## Not run: 
### Minimal (with many features missing) reimplementation 
#     of corrMatrix() terms as a corrFamily 


corrMatrix_cF <- function(corrMatrix) {
  
  force(corrMatrix) # Makes it available in the environment of the functions next defined.
  oldZlevels <- NULL
  
  initialize <- function(Zmatrix, ...) {
    oldZlevels <<- colnames(Zmatrix) # Pass info about levels of the random effect in the data.
  }
  
  Cf <- function(newlevels=oldZlevels ) {
    if (length(newlevels)) {
      corrMatrix[newlevels,newlevels]
    } else corrMatrix[oldZlevels,oldZlevels] # for Cf(tpar=numeric(0L))
  }
  
  calc_moreargs <- function(corrfamily, ...) {
    list(init=c(),lower=c(),upper=c())
  }
  
  make_new_corr_lists <- function(newLv_env, which_mats, newZAlist, new_rd, ...) {
    newlevels <- colnames(newZAlist[[new_rd]])
    newLv_env$cov_newLv_oldv_list[[new_rd]] <- corrMatrix[newlevels,oldZlevels, drop=FALSE]
    if (which_mats$nn[new_rd]) {
      newLv_env$cov_newLv_newLv_list[[new_rd]] <- corrMatrix[newlevels,newlevels, drop=FALSE]
    } else { 
      newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- rep(1,length(newlevels)) 
    }
  }
  
  list(Cf=Cf, tpar=numeric(0L), initialize=initialize, calc_moreargs=calc_moreargs, 
       make_new_corr_lists=make_new_corr_lists,
       tag="corrMatrix_cF") 
}

register_cF("corrMatrix_cF")

# usage:

data("blackcap") 
MLcorMat <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
                       nu=0.6285603,rho=0.0544659)
corrmat <- proxy::as.matrix(MLcorMat, diag=1)

fitme(migStatus ~ means+ corrMatrix_cF(1|name, corrMatrix=corrmat),data=blackcap,
      corrMatrix=MLcorMat,method="ML")
      
unregister_cF("corrMatrix_cF") # Tidy things before leaving.         


## End(Not run)

Description

Tentative formal rules for definition of corrFamily descriptors (work in progress). This is likely to repeat and extend information partially given in corrFamily and corrFamily-design documentations.

User-level rules (not relevant fo corrFamily descriptors internally modified during a fit):

tpar

Should always be present. For trivial parameterless cases (e.g. ranGCA), it should be numeric(0L), not NULL.

Cf

function; should always be present. For trivial uncorrelated random effects (e.g. ranGCA, where only the Af function carries the information for the model), it should return an identity matrix, not NULL, with row names to be matched to the column names of the Z matrix for the random effect.

calc_moreargs

optional function. If present, it should have formal arguments including at least corrfamily and ....

Af

function; optional. If present, it should have row names to be matched to the column names of the Z matrix for the random effect, and also needs column names if it is to be matched with the row names of a correlation matrix (or its inverse).

initialize

Optional function. If present, should have formal arguments including at least Zmatrix and ....

In predefined corrFamily constructors, variables created by initialize for use by Cf or Af should be declared (typically as NULL) in the body of the constructor, so that R CMD check does not complain.

public

An environment. initialize may write into it. It might also read into it, for example read the result of a long previous computation by initialize during a previous fit, though this opens the door to various errors.

Description

This documentation describe additional design features to be taken into account when defining a new corrFamily descriptor for a correlation model. Using such a descriptor will be more efficient than the equally general method, of maximizing an objective function of the correlation parameters that calls (say) fitme() on a model including a corrMatrix itself function of the correlation parameters. But this may still be inefficient if a few issues are ignored.

For elements of the corrFamily descriptor for basic cases:

Cf

The function value should (1) be of constant class for all parameter values. For families of mathematically sparse matrices, the CsparseMatrix class is recommended (and more specifically the dsCMatrix class since the matrix is symmetric); (2) have row names that match the levels of the grouping factor (the nested random effect Example shows the code needed when this nested effect is defined from two variables).

tpar

In order to favor the automatic selection of suitable algorithms, tpar should be chosen so that Cf(tpar) is least sparse (i.e., has the minimal number of elements equal to zero) in the correlation family, in terms of its sparsity and of the sparsity of its inverse. A tpar yielding an identity matrix is often a *bad* template as least sparse correlation matrices and their inverses are denser for most families except diagonal ones. For degerate corrFamily objects that describe a constant correlation model rather than a parametric family, use tpar=numeric(0).

type

Do not forget type="precision" it if the return value of Cf is an inverse correlation matrix rather than a correlation matrix, in which case one should specify .

calc_moreargs

should have formal arguments including at least corrfamily and .... The source code of ARp, ARMA or diallel shows the expected structure of its return value.

For advanced features of the corrFamily descriptor:

Af

Af has (minimally) three formal arguments (newdata, term, ...). spaMM will call Af with distinct values of the newdata argument for the fit, and for predictions for new data. For the curious: the term argument that will be provided by spaMM to Af is the formula term for the random effect – an object of class call, as obtained e.g. by
( ~ 1+ corrFamily(1 | longitude + latitude))[[2]][[3]] –, which will provide the names of the variables that need to be taken from the newdata to construct the matrix returned by Af.

Details

• spaMM will regularize invalid or nearly-singular correlation or covariance matrices internally if the correlation function has not done so already, but it it better to control this in the correlation function. The regularize convenience function is available for that purpose, but parametrizations that avoid the need for regularization are even better, since fitting models with nearly-singular correlation matrices is prone to various difficulties (The Toeplitz example below is good to illustrate potential problems but is otherwise poor as it produces non-positive definite matrices; the ARp constructor illustrates a parametrization that avoids that problem).

• Users should make sure that any regularized matrix still belongs to the intended parametric family of matrices, and they should keep in mind that the spaMM output will show the input parameters of the unregularized matrix, not the parameters of the regularized one (e.g., in the Toeplitz example below, the fitted matrix is a regularized Toepliz matrix with slightly different coefficients than the input parameters).

  And for efficiency,

• Let us repeat that the correlation function should return matrices of constant class, and in sparse format when the matrices are indeed mathematically sparse. For mathematically dense matrices (as in the Toeplitz example below), the dsyMatrix class may be suitable.

• Let us repeat that in order to favor the automatic selection of suitable algorithms, tpar should be chosen so that Cf(tpar) is least sparse in the correlation family. For matrices of CsparseMatrix, a check is implemented to catch wrong choices of tpar.

• For challenging problems (large data, many parameters...) it may pay to optimize a bit the correlation function. The Example of nested effects with heterogenous variance below illustrates a possible trick. In the same cases, It may also pay to try the alternative algebraic methods, by first comparing speed of the different methods (control.HLfit=list(algebra= <"spprec"|"spcorr"|"decorr">)) for given correlation parameter values, rather than to assume that spaMM will find the best method (even if it often does so).

• The corrFamily descriptor may optionally contain booleans possiblyDenseCorr and sparsePrec to help spaMM select the most appropriate matrix algebraic methods. sparsePrec should be set to TRUE if sparse-precision methods are expected to be efficient for fitting the random effect. possiblyDenseCorr should be set to FALSE if the correlation matrix is expected to be sparse, which means here that less than 15% of its elements are non-zero.

Examples

if (spaMM.getOption("example_maxtime")>2 &&
      requireNamespace("agridat", quietly = TRUE)) {

data("onofri.winterwheat", package="agridat")

##### Fitting a Toeplitz correlation model for temporal correlations #####

# A Toeplitz correlation matrix of dimension d*d has d-1 parameters 
# (by symmetry, and with 1s on the main diagonal). These d-1 parameters 
# can be fitted as follows:

Toepfn <- function(v) {
  toepmat <- Matrix::forceSymmetric(toeplitz(c(1,v))) # dsyMatrix
  # Many of the matrices in this family are not valid correlation matrices;
  #   the regularize() function is handy here:
  toepmat <- regularize(toepmat, maxcondnum=1e12)
  # And don't forget the rownames!
  rownames(toepmat) <- unique(onofri.winterwheat$year)
  toepmat
}

(Toepfit <- spaMM::fitme(
  yield ~ gen + corrFamily(1|year), data=onofri.winterwheat, method="REML",
  covStruct=list(corrFamily=list(Cf=Toepfn, tpar=rep(1e-4,6))), 
         # (Note the gentle warning if one instead uses tpar=rep(0,6) here)  
  lower=list(corrPars=list("1"=rep(-0.999,6))), 
  upper=list(corrPars=list("1"=rep(0.999,6))))) 

# The fitted matrix is (nearly) singular, and was regularized:

eigen(Corr(Toepfit)[[1]])$values

# which means that the returned likelihood may be inaccurate, 
# and also that the actual matrix elements differ from input parameters:

Corr(Toepfit)[[1]][1,-1]  

### The usual rules for specifying covStruct, 'lower', 'upper' and 'init' apply
# here when the corrFamily term is the second random-effect:

(Toep2 <- spaMM::fitme(
        yield ~ 1 + (1|gen) + corrFamily(1|year), data=onofri.winterwheat, method="REML",
        covStruct=list("1"=NULL, corrFamily=list(Cf=Toepfn, tpar=rep(1e-4,6))), 
        , init=list(corrPars=list("2"=rep(0.1,6))),
        lower=list(corrPars=list("2"=rep(-0.999,6))), 
        upper=list(corrPars=list("2"=rep(0.999,6)))))

##### Fitting one variance among years per each of 8 genotypes. #####

# First, note that this can be *more efficiently* fitted by another syntax:

### Fit as a constrained random-coefficient model: 
    
# Diagonal matrix of NA's, represented as vector for its lower triangle:
ranCoefs_for_diag <- function(nlevels) { 
  vec <- rep(0,nlevels*(nlevels+1L)/2L)
  vec[cumsum(c(1L,rev(seq(nlevels-1L)+1L)))] <- NA
  vec
} 
    
(by_rC <- spaMM::fitme(yield ~ 1 + (0+gen|year), data=onofri.winterwheat, method="REML",
                       fixed=list(ranCoefs=list("1"=ranCoefs_for_diag(8)))))
                         
### Fit as a corrFamily model:   

gy_levels <- paste0(gl(8,1,length =56,labels=levels(onofri.winterwheat$gen)),":",
                        gl(7,8,labels=unique(onofri.winterwheat$year)))
                        
# A log scale is often suitable for variances, hence is used below;

# a correct but crude implementation of the model is
diagf <- function(logvar) {
  corr_map <- kronecker(Matrix::.symDiagonal(n=7),diag(x=exp(logvar)))
  rownames(corr_map) <- gy_levels
  corr_map
}

# but we can minimize matrix operations as follows:

corr_map <- Matrix::.symDiagonal(n=8,x=seq(8))
rownames(corr_map) <- unique(onofri.winterwheat$gen)
      
diagf <- function(logvar) {
  corr_map@x <- exp(logvar)[corr_map@x]
  corr_map
}                 # (and this returns a dsCMatrix)
      
(by_cF <- spaMM::fitme(
        yield ~ 1 + corrFamily(1|gen %in% year), data=onofri.winterwheat, method="REML",
        covStruct=list(corrFamily = list(Cf=diagf, tpar=rep(1,8))), 
        fixed=list(lambda=1),            # Don't forget to fix this redundant parameter!
        # init=list(corrPars=list("1"=rep(log(O.1),8))), # 'init' optional 
        lower=list(corrPars=list("1"=rep(log(1e-6),8))), # 'lower' and 'upper' required
        upper=list(corrPars=list("1"=rep(log(1e6),8)))))

# => The 'gen' effect is nested in the 'year' effect and this must be specified in the 
# right-hand side of corrFamily(1|gen %in% year) so that the design matrix 'Z' for the 
# random effects to have the correct structure. And then, as for other correlation
# structures (say Matern) it should be necessary to specify only the correlation matrix 
# for a given year, as done above. Should this fail, it is also possible to specify the 
# correlation matrix over years, as done below. spaMM will automatically detect, from   
# its size matching the number of columns of Z, that it must be the matrix over years.  

corr_map <- Matrix::forceSymmetric(kronecker(Matrix::.symDiagonal(n=7),diag(x=seq(8))))
rownames(corr_map) <- gy_levels

diagf <- function(logvar) {
  corr_map@x <- exp(logvar)[corr_map@x]
  corr_map
}                 # (and this returns a dsCMatrix)

(by_cF <- spaMM::fitme(
  yield ~ 1 + corrFamily(1|gen %in% year), data=onofri.winterwheat, method="REML",
  covStruct=list(corrFamily = list(Cf=diagf, tpar=rep(1,8))), 
  fixed=list(lambda=1),            # Don't forget to fix this redundant parameter!
  # init=list(corrPars=list("1"=rep(log(O.1),8))), # 'init' optional 
  lower=list(corrPars=list("1"=rep(log(1e-6),8))), # 'lower' and 'upper' required
  upper=list(corrPars=list("1"=rep(log(1e6),8)))))

exp(get_ranPars(by_cF)$corrPars[[1]]) # fitted variances 
  
}

Description

This was the first function for fitting all spatial models in spaMM, and is still fully functional, but it is recommended to use fitme which has different defaults and generally selects more efficient fitting methods, and will handle all classes of models that spaMM can fit, including non-spatial ones. corrHLfit performs the joint estimation of correlation parameters, fixed effect and dispersion parameters.

Usage

corrHLfit(formula, data, init.corrHLfit = list(), init.HLfit = list(), 
          ranFix, fixed=list(), lower = list(),  upper = list(), 
          objective = NULL, resid.model = ~1, 
          control.dist = list(), control.corrHLfit = list(),
          processed = NULL, family = gaussian(), method="REML",
          nb_cores = NULL, weights.form = NULL, ...)

Arguments

for default values.

Details

For approximations of likelihood, see method. For the possible structures of random effects, see random-effects,

By default corrHLfit will estimate correlation parameters by maximizing the objective value returned by HLCor calls wherein the dispersion parameters are estimated jointly with fixed effects for given correlation parameters. If dispersion parameters are specified in init.corrHLfit, they will also be estimated by maximizing the objective value, and HLCor calls will not estimate them jointly with fixed effects. This means that in general the fixed effect estimates may vary depending on init.corrHLfit when any form of REML correction is applied.

Correctly using corrHLfit for likelihood ratio tests of fixed effects may then be tricky. It is safe to perform full ML fits of all parameters (using method="ML") for such tests (see Examples). The higher level function fixedLRT is a safe interface for likelihood ratio tests using some form of REML estimation in corrHLfit.

attr(<fitted object>,"optimInfo")$lower and ...$upper gives the lower and upper bounds for optimization of correlation parameters. These are the default values if the user did not provide explicit values. For the adjacency model, the default values are the inverse of the maximum and minimum eigenvalues of the adjMatrix. For the Matérn model, the default values are not so easily summarized: they are intended to cover the range of values for which there is statistical information to distinguish among them.

Value

The return value of an HLCor call, with additional attributes. The HLCor call is evaluated at the estimated correlation parameter values. These values are included in the return object as its $corrPars member. The attributes added by corrHLfit include the original call of the function (which can be retrived by getCall(<fitted object>), and information about the optimization call within corrHLfit.

See Also

See more examples on data set Loaloa, to compare fit times by corrHLfit and fitme. See fixedLRT for likelihood ratio tests.

Examples

# Example with an adjacency matrix (autoregressive model):
if (spaMM.getOption("example_maxtime")>0.7) {          
  corrHLfit(cases~I(prop.ag/10) +adjacency(1|gridcode)+offset(log(expec)),
          adjMatrix=Nmatrix,family=poisson(),data=scotlip,method="ML") 
}

#### Examples with Matern correlations
## A likelihood ratio test based on the ML fits of a full and of a null model.
if (spaMM.getOption("example_maxtime")>1.4) {
 data("blackcap")
 (fullfit <- corrHLfit(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
                    method="ML") )
 (nullfit <- corrHLfit(migStatus ~ 1 + Matern(1|longitude+latitude),data=blackcap,
                    method="ML",init.corrHLfit=list(phi=1e-6))) 
 ## p-value:
 1-pchisq(2*(logLik(fullfit)-logLik(nullfit)),df=1)
}

Description

corrMatrix is an argument of HLCor, of class dist or matrix, with can be used if the model formula contains a term of the form corrMatrix(1|<...>). It describes a correlation matrix, possibly as a dist object. A covariance matrix can actually be passed through this argument, but then it must be a full matrix, not a dist object. The way the rows and columns of the matrix are matched to the rows of the data depends on the nature of the grouping term <...>.

The covStruct argument can be used for the same purpose and is much more general, in particular allowing to specify several correlation matrices.

Details

The simplest case is illustrated in the first two examples below: the grouping term is identical to a single variable which is present in the data, whose levels match the rownames of the corrMatrix. As illustrated by the second example, the order of the data does not matter in that case, because the factor levels are used to match the data rows to the appropriate row and columns of the corrMatrix. The corrMatrix may even contain rows (and columns) in excess of the levels of the grouping term, in which case these rows are ignored.

These convenient properties no longer hold when the grouping term is not a single variable from the data (third example below), or when its levels do not correspond to row names of the matrix. In these cases, (1) no attempt is made to match the data rows to the row and column names of the corrMatrix. Such attempt could succeed only if the user had given names to the matrix matching those that the called function could create from the information in the data, in which case the user should find easier to specify a single variable that can be matched; (2) the order of data and corrMatrix matter; Internally, a single factor variable is constructed from all levels of the variables in the grouping term (i.e., from all levels of latitude and longitude, in the third example), with levels 1,2,3... that are matched to rows 1,2,3... of the corrMatrix. Thus the first row of the data is always associated to the first row of the matrix; (3) further, the dimension of the matrix must match the number of levels implied by the grouping term. For example, one might consider the case of 14 response values but of correlations between only 7 levels of a random effect, with two responses for each level. Then the matrix must be of dimension 7x7.

Examples

data("blackcap") 
## Here we manually reconstruct the correlation matrix 
##  of the ML fit produced by corrHLfit:
MLcorMat <- MaternCorr(proxy::dist(blackcap[,c("longitude","latitude")]),
                        nu=0.6285603,rho=0.0544659)
blackcap$name <- as.factor(rownames(blackcap))                
#

## (1) Single variable present in the data 
#
HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
      corrMatrix=MLcorMat,method="ML")

## (2) Same, permuted: still gives correct result
#
perm <- sample(14)
# Permuted matrix (with permuted names) as 'dist' object
pmat <- as.matrix(MLcorMat)[perm,perm] 
HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
      corrMatrix=as.dist(pmat),method="ML")
#
# Permuted matrix (with permuted names) as correlation matrix
pcorr <- proxy::as.matrix(MLcorMat, diag=1)[perm,perm] 
HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
      corrMatrix=pcorr,method="ML")
#

## (3) Other grouping terms (note the messages):
#
HLCor(migStatus ~ means+ corrMatrix(1|longitude+latitude),data=blackcap,
      corrMatrix=MLcorMat,method="ML")

Description

covStruct is a formal argument of HLCor, also handled by fitme and corrHLfit, that allows one to specify the correlation structure for different types of random effects, It is an alternative to other ad hoc formal arguments such as corrMatrix or adjMatrix. It replaces the deprecated function Predictor(...) which has served as an interface for specifying the design matrices for random effects in early versions of spaMM.

The main use of covStruct is to specify the correlation matrix of levels of a given random effect term, or its inverse (a precision matrix). Assuming that the design matrix of each random effect term follows the structure ZAL described in random-effects, it is thus an indirect way of specifying the “square root” L of the correlation matrix. The optional A factor can also be given by the optional "AMatrices" attribute of covStruct.

covStruct is a list of matrices with names specifying the type of matrix considered:
covStruct=list(corrMatrix=<some matrix>) or covStruct=list(adjMatrix=<some matrix>), where the “corrMatrix” or “adjMatrix” labels are used to specify the type of information provided (accordingly, the names can be repeated: covStruct=list(corrMatrix=<.>,corrMatrix=<.>)). NULL list members may be necessary, e.g.
covStruct=list(corrMatrix=<.>,"2"=NULL,corrMatrix=<.>))
when correlations matrices are required only for the first and third random effect.

The covariance structure of a corrMatrix(1|<grouping factor>) formula term can be specified in two ways (see Examples): either by a correlation matrix factor (covStruct=list(corrMatrix=<some matrix>)), or by a precision matrix factor Q such that the covariance factor is $\lambda$Q$^{-1}$, using the type name "precision": covStruct=list(precision=<some matrix>). The function as_precision can be used to perform the conversion from correlation information to precision factor (using a crude solve() that may not always be efficient), but fitting functions may also perform such conversions automatically.

"AMatrices" is a list of matrices. The names of elements of the list does not matter, but the ith A matrix, and its row names, should match the ith Z matrix, and its column names. This implies that NULL list members may be necessary, as for the covStruct list.

Usage

as_precision(corrMatrix, condnum=1e12)

Arguments

Details

covStruct can also be specified as a list with an optional "types" attribute, e.g.
structure(list(<some matrix>,types="corrMatrix")).

Value

as_precision returns a list with additional class precision and with single element a symmetric matrix of class dsCMatrix.

See Also

Gryphon and pedigree for a type of applications where declaring a precision matrix is useful.

Examples

## Not run: 
data("blackcap") 
# a 'dist' object can be used to specify a corrMatrix:  
MLdistMat <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
                        nu=0.6285603,rho=0.0544659) # a 'dist' object!
blackcap$name <- as.factor(rownames(blackcap))     
fitme(migStatus ~ means + corrMatrix(1|name), data=blackcap,
      corrMatrix=MLdistMat)

#### Same result by different input and algorithm:
fitme(migStatus ~ means + corrMatrix(1|name), data=blackcap,
      covStruct=list(precision=as_precision(MLdistMat)))

# Manual version of the same:
as_mat <- proxy::as.matrix(MLdistMat, diag=1) 
prec_mat <- solve(as_mat) ## precision factor matrix
fitme(migStatus ~ means + corrMatrix(1|name), data=blackcap,
      covStruct=list(precision=prec_mat))

# Since no correlation parameter is estimated, 
# HLcor(., method="ML")  is here equivalent to fitme()

## End(Not run)

Description

ranGCA and diallel are random-effect structures designed to represent the effet of symmetric interactions between pairs of individuals (order of individuals in the pair does not matter), while antisym represents anti-symmetric interactions (the effect of reciprocal ordered pairs on the outcome are opposite, as in the so-called Bradley-Terry models). These random-effect structures all account for multiple membership, i.e., the fact that the same individual may act as the first or the second individual among different pairs, or even within one pair if this makes sense).

More formally, the outcome of an interaction between a pair $i,j$ of agents is subject to a symmetric overall random effect $v_{ij}$ when the effect “on” individual $i$ (or viewed from the perspective of individual $i$) equals the effect on $j$: $v_{ij}=v_{ji}$. This may result from the additive effect of individual random effects $v_{i}$ and $v_{j}$: $v_{ij}=v_i+v_j$, but also from non-additive effects $v_{ij}=v_i+v_j+a_{ij}$ if the interaction term $a_{ij}$ is itself symmetric ($a_{ij}=a_{ji}$). ranGCA and diallel effects represent such symmetric effects, additive or non-additive respectively, in a model formula (see Details for the semantic origin of these names and how they can be changed). Conversely, antisymmetry is characterized by $v_{ij}=v_i-v_j=-v_{ji}$ and is represented by the antisym formula term.

If individual-level random effects of the form (1|ID1)+ (1|ID2) were included in the model formula instead of ranGCA(1|ID1+ID2) for symmetric additive interactions, this would result in different variances being fitted for each random effect (breaking the assumption of symmetry), and the value of the random effect would differ for an individual whether it appears as a level of the first random effect or of the second (which is also inconsistent with the idea that the random effect represents a property of the individual).

When ranGCA or antisym random effects are fitted, the individual effects are inferred. By contrast, when a diallel random effect is fitted, an autocorrelated random effect $v_{ij}$ is inferred for each unordered pair (no individual effect is inferred), with correlation $\rho$ between levels for pairs sharing one individual. This correlation parameter is fitted and is constrained by $\rho < 0.5$ (see Details). ranGCA is equivalent to the case $\rho= 0.5$. diallel fits can be slow for large data if the correlation matrix is large, as this matrix can have a fair proportion of nonzero elements. There may also be identifiability issues for variance parameters: in a LMM as shown in the examples, there will be three parameters for the random variation (phi, lambda and rho) but only two can be estimated if only one observation is made for each dyad.

Usage

## formula terms:

# ranGCA(1| <.> + <.>) 
# antisym(1| <.> + <.>) 
# diallel(1|  <.> + <.>, tpar, fixed = NULL, public = NULL)

## where the <.> are two factor identifiers, ** whose levels
## must be identical when representing the same individual **

## corrFamily constructors:
ranGCA() # no argument
antisym() # no argument
diallel(tpar, fixed = NULL, public = NULL)

Arguments

Details

Although the symmetric random-effect structures may be used in many different contexts (including social network analysis, or “round robin” experiments in psychology; another possibly relevant literature keyword here is “multi membership”), their present names refer to the semantics established for diallel experiments (e.g., Lynch & Walsh, 1998, p. 611), because it is not easy to find a more general yet intuitive semantics. If the names ranGCA and diallel sound inappropriate for your context of application, you can declare and use an alternative name for them, taking advantage of the fact that they are random-effect structures defined through corrFamily constructors, which are functions named as the formula term. For example, symAdd(1|ID1+ID2) can be used in a model formula after the following two steps:

# Define the 'symAdd' corrFamily constructor (a function) by copy:
symAdd <- ranGCA     
# Associate the 'symAdd' function to 'symAdd' formula terms:
register_cF("symAdd")  

In diallel experiments, one analyzes the phenotypes of offspring from multiple crosses among which the mother in a cross can be the father in another, so this is an example of multiple-membership. The additive genetic effects of each parent's genotypes are described as “general combining abilities” (GCA). In case of non-additive effect, the half-sib covariance is not half the full-sib covariance and this is represented by the interaction $a_{ij}$ described as “specific combining abilities” (SCA). The sum of GCA and SCA defines a synthetic random effect “received” by the offspring, with distinct levels for each unordered parental pair, and with correlation $\rho$ between effects received by half-sibs (one shared parent). $\rho$ corresponds to var(GCA)/[2*var(GCA)+var(SCA)] and is necessarily $\le 0.5$.

See the X.GCA documentation for similar constructs for fixed effects.

Value

The functions return corrFamily descriptors whose general format is described in corrFamily. The ones produced by ranGCA and antisym are atypical in that only their Af element is non-trivial.

References

Lynch, M., Walsh, B. (1998) Genetics and analysis of quantitative traits. Sinauer, Sunderland, Mass.

Examples

#### Simulate dyadic data

set.seed(123)
nind <- 10       # Beware data grow as O(nind^2)
x <- runif(nind^2) 
id12 <- expand.grid(id1=seq(nind),id2=seq(nind))
id1 <- id12$id1
id2 <- id12$id2
u <-  rnorm(nind,mean = 0, sd=0.5)

## additive individual effects:
y <-  0.1 + 1*x + u[id1] +  u[id2] + rnorm(nind^2,sd=0.2)

## Same with non-additive individual effects:
dist.u <- abs(u[id1] -  u[id2])
z <- 0.1 + 1*x + dist.u + rnorm(nind^2,sd=0.2)

## anti-smmetric individual effects:
t <-  0.1 + 1*x + u[id1] - u[id2] + rnorm(nind^2,sd=0.2)

dyaddf <- data.frame(x=x, y=y, z=z, t=t, id1=id1,id2=id2)
# : note that this contains two rows per dyad, which avoids identifiability issues.

# Enforce that interactions are between distinct individuals (not essential for the fit):
dyaddf <- dyaddf[- seq.int(1L,nind^2,nind+1L),] 


# Fits:

(addfit <- fitme(y ~x +ranGCA(1|id1+id2), data=dyaddf))
#
# practically equivalent to:
#
(fitme(y ~x +diallel(1|id1+id2, fixed=0.49999), data=dyaddf))

(antifit <- fitme(t ~x +antisym(1|id1+id2), data=dyaddf))

(distfit <- fitme(z ~x +diallel(1|id1+id2), data=dyaddf))

Description

This experimental function displays information about parameter values for which some numerical problems have occurred. Some warnings suggest its use.

Numerical problems may occur if the information matrix (for the augmented linear model used in the iteratively reweighted least-squares algorithm) is nearly singular. spaMM may try to check whether such singularity occurs when this algorithm has not converged. But itself may be slow so it is not performed systematically for large matrices. spaMM.options(diagnose_conv=<integer>) may be used to control the maximum size of matrices for which the check is performed.

When “outer” generic optimization is performed, information is reported about the range of parameter values for which problems occurred, (see Value). The fit object divinfo element may also contain more informative tables of parameter values. This information is currently missing for “inner”-optimized parameters.

Usage

div_info(object, ...)

Arguments

Value

Used mainly for the side effects (printed output) but returns invisibly either a single parameter vector (if a single numerical problem occurred) or a matrix of parameter ranges, or NULL if there is no problem to report.

Examples

# Tragically ;-), no simple example of numerical problems 
# that can be diagnosed by div_info() is currently available.

Description

This extracts the number of degrees of freedom for a model, in the usual sense for likelihood-ratio tests: a count of number of fitted parameters, distinguishing different classes of parameters (see Value).

Usage

DoF(object)

Arguments

Details

The output distinguishes counts of random-effect vs residual-dispersion parameters, following the conceptual distinction between effects that induce correlations between different levels of the resonse vs. observation-level effects. However, a residual-dispersion component can be declared as a random effect, so that the counts for logically equivalent models may differ according to the way a model was declared. For example if residual dispersion for an LLM is declared as an observation-level random effect while phi is fixed, the p_lambda component will include 1 df for what would otherwise be accounted by the p_rdisp component. A more involved case where the same contrast happens is when a negative-binomial model (with a residual-dispersion shape parameter) is declared as a Poisson-gamma mixture model (with a varaince parameter for the Gamma-distributed individual-level random effect).

Value

A vector with possible elements p_fixef, p_lambda, p_corrPars and p_rdisp for, respectively, the number of fixed-effect coefficients of the main-response model, the number of random-effect variance parameters, the number of random-effect correlation parameters, and the number of residual dispersion parameters (the latter being itself, for a mixed-effect residual-dispersion model, the sum of such components).

See Also

df.residual.HLfit; get_any_IC for extracting effective degrees of freedom considered in the model-selection literature; as_LMLT for access to the effective degrees of freedom considered in Satterthwaite's test and its extentions.

Description

interface to apply some function fn in parallel on columns of a matrix. It is not logically restricted to mixed-effect applications, hence it can be used more widely. Depending on the nb_cores argument, parallel or serial computation is performed, calling the future.apply::future_apply function. A socket cluster is used by default for parallel computations, but a fork cluster can be requested on linux and alike operating systems by using argument cluster_args=list(type="FORK").

Usage

dofuture(newresp, fn, nb_cores=NULL, fit_env, control=list(), 
      cluster_args=NULL, debug.=FALSE, iseed=NULL, 
      showpbar="ignored", pretest_cores="ignored",
      ... )

Arguments

Value

The result of calling future.apply. If the progressr package is loaded, a side-effect of dofuture is to show a progress bar with character 'S' or 'P' or 'F' depending on parallelisation status (serial/socket/fork).

See Also

dopar for an alternative implementation of (essentially) the same functionalities, and wrap_parallel for its differences from dofuture.

Examples

## Not run: 
if (requireNamespace("future.apply", quietly = TRUE)) {

 # Useless function, but requiring some argument beyond the first
 foo <- function(y, somearg, ...) {
   if ( is.null(somearg) || TRUE ) length(y)
 }

 # Whether FORK can be used depends on OS and whether Rstudio is used:
   dofuture(matrix(1,ncol=4,nrow=3), foo, fit_env=list(), somearg=NULL, 
     nb_cores=2, cluster_args=list(type="FORK"))
}

## End(Not run)

Description

dopar and combinepar are interfaces primarily designed to apply some function fn in parallel on columns of a matrix, although other uses are possible. Depending on the nb_cores argument, parallel or serial computation is performed. A socket cluster is used by default for parallel computations, but a fork cluster can be requested on linux and alike operating systems by using argument cluster_args=list(type="FORK").

dopar has been designed to provide by default a progress bar in all evaluations contexts. A drawback is that different procedures are called depending e.g. on the type of cluster, with different possible controls. In particular, foreach is called in some cases but not others, so non-trivial values of its .combine control are not always enforced. The alternative interface combinepar will always use foreach, and will still try to provide by default a progress bar but may fail to do so in some cases (see Details).

Usage

dopar(newresp, fn, nb_cores = NULL, fit_env, 
      control = list(.final=function(v) if( ! is.list(v[[1]])) {do.call(cbind,v)} else v), 
      cluster_args = NULL, debug. = FALSE, iseed = NULL, 
      showpbar = eval(spaMM.getOption("barstyle")), 
      pretest_cores =NULL, ...)
combinepar(newresp, fn, nb_cores = NULL, cluster=NULL, fit_env, 
      control = list(.final=function(v) if( ! is.list(v[[1]])) {do.call(cbind,v)} else v), 
      cluster_args = NULL, debug. = FALSE, iseed = NULL, 
      showpbar = eval(spaMM.getOption("barstyle")), 
      pretest_cores =NULL, ...)

Arguments

Details

Control of random numbers through the "L'Ecuyer-CMRG" generator and the iseed argument is not sufficient for consistent results when the doSNOW parallel backend is used, so if you really need such control in a fn using random numbers, do not use doSNOW. Yet, it is fine to use doSNOW for bootstrap procedures in spaMM, because the fitting functions do not use random numbers: only sample simulation uses them, and it is not performed in parallel.

combinepar calls foreach::%dopar% which assumes that a cluster has been declared using a suitable backend such as doSNOW, doFuture or doParallel. If only the latter is available, no progress bar is displayed. A method to render a bar when doParallel is used can be found on the Web, but that bar is not a valid progress bar as it is displayed only after all the processes have been run.

Value

The result of calling foreach, pbapply or mclapply, as dependent on the control argument and the interface used. A side-effect of either interface is to show a progress bar whose character informs about the type of parallelisation performed: a "F" or default "=" character for fork clusters, a "P" for parallel computation via foreach and doSNOW, a "p" for parallel computation via foreach and doFuture or via pbapply, and "s" for serial computation foreach and doParallel or via pbapply.

See Also

dofuture is yet another interface with (essentially) the same functionalities as dopar. See the documentation of the wrap_parallel option for its differences from dopar.

Examples

## See source code of spaMM_boot()

## Not run: 
# Useless function, but requiring some argument beyond the first
foo <- function(y, somearg, ...) {
  if ( is.null(somearg) || TRUE ) length(y)
}

# Whether FORK can be used depends on OS and whether Rstudio is used:
dopar(matrix(1,ncol=4,nrow=3), foo, fit_env=list(), somearg=NULL, 
  nb_cores=2, cluster_args=list(type="FORK"))

## End(Not run)