copCAR {copCAR}R Documentation

Fit copCAR model to discrete areal data.


Fit the copCAR model to Bernoulli, negative binomial, or Poisson observations.


copCAR(formula, family, data, offset, A, method = c("CML", "DT", "CE"),
  confint = c("none", "bootstrap", "asymptotic"), model = TRUE, x = FALSE,
  y = TRUE, verbose = FALSE, control = list())



an object of class “formula” (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of the model specification are given under "Details".


the marginal distribution of the observations at the areal units and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See family for details of family functions.) Supported families are binomial, negbinomial, and poisson. When the negative binomial family is used, an initial value for \theta must be passed to the negbinomial family function.


an optional data frame, list or environment (or object coercible by to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which copCAR is called.


this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of observations. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See model.offset.


the symmetric binary adjacency matrix for the underlying graph.


the method for inference. copCAR supports the continous extension (“CE”), distributional transform (“DT”), and composite marginal likelihood (“CML”).


the method for computing confidence intervals. This defaults to “none”. The other options are “bootstrap” (for parametric bootstrap intervals using the quantile method) and “asymptotic” (for intervals computed using an estimate of the asymptotic covariance matrix).


a logical value indicating whether the model frame should be included as a component of the returned value.


a logical value indicating whether the model matrix used in the fitting process should be returned as a component of the returned value.


a logical value indicating whether the response vector used in the fitting process should be returned as a component of the returned value.


a logical value indicating whether to print various messages to the screen, including progress updates. Defaults to FALSE.


a list of parameters for controlling the fitting process.


the size of the (parametric) bootstrap sample. This applies when confint = "bootstrap", or when confint = "asymptotic" and method = "CML" or method = "DT". Defaults to 500.


the number of vectors of standard uniforms used to approximate the expected likelhood when method = "CE". Defaults to 1000.


the value \rho^{\max}, which is the maximum value of \rho used to approximate the CAR variances when method = "CE" or method = "DT". If missing, assumed to be 0.999.


the tolerance \epsilon > 0 used to approximate the CAR variances when method = "CE" or method = "DT". If missing, assumed to be 0.01.

itemmaxitthe maximum number of iterations to be used by optim when optimizing the objective function. Defaults to 1000.


a logical value indicating whether to parallelize the bootstrap. This defaults to TRUE if the parallel package can be loaded.


the cluster type, one of “SOCK” (default), “PVM”, “MPI”, or “NWS”.


the number of slave nodes to create.


This function performs frequentist inference for the copCAR model proposed by Hughes (2015). copCAR is a copula-based areal regression model that employs the proper conditional autoregression (CAR) introduced by Besag, York, and MolliƩ (1991). Specifically, copCAR uses the CAR copula, a Caussian copula based on the proper CAR.

The spatial dependence parameter \rho \in [0, 1), regression coefficients \beta = (\beta_1, \dots, \beta_p)' \in R^p, and, for negative binomial margins, dispersion parameter \theta>0 can be estimated using the continous extension (CE) (Madsen, 2009), distributional transform (DT) (Kazianka and Pilz, 2010), or composite marginal likelihood (CML) (Varin, 2008) approaches.

The CE approach transforms the discrete observations to continous outcomes by convolving them with independent standard uniforms (Denuit and Lambert, 2005). The true likelihood for the discrete outcomes is the expected likelihood for the transformed outcomes. An estimate (sample mean) of the expected likelihood is optimized to estimate the copCAR parameters. The number of standard uniform vectors, m, can be chosen by the user. The default value is 1,000. The CE approach is exact up to Monte Carlo standard error but is computationally intensive (the computational burden grows rapidly with increasing m). The CE approach tends to perform poorly when applied to Bernoulli outcomes, and so that option is not permitted.

The distributional transform stochastically "smoothes" the jumps of a discrete distribution function (Ferguson, 1967). The DT-based approximation (Kazianka and Pilz, 2010) for copCAR performs well for Poisson and negative binomial marginals but, like the CE approach, tends to perform poorly for Bernoulli outcomes.

The CML approach optimizes a composite marginal likelihood formed as the product of pairwise likelihoods of adjacent observations. This approach performs well for Bernoulli, negative binomial, and Poisson outcomes.

In the CE and DT approaches, the CAR variances are approximated. The quality of the approximation is determined by the values of control parameters \epsilon > 0 and \rho^{\max} \in [0, 1). The default values are 0.01 and 0.999, respectively.

When confint = "bootstrap", a parametric bootstrap is carried out, and confidence intervals are computed using the quantile method. Monte Carlo standard errors (Flegal et al., 2008) of the quantile estimators are also provided.

When confint = "asymptotic", confidence intervals are computed using an estimate of the asymptotic covariance matrix of the estimator. For the CE method, the inverse of the observed Fisher information matrix is used. For the CML and DT methods, the objective function is misspecified, and so the asymptotic covariance matrix is the inverse of the Godambe information matrix (Godambe, 1960), which has a sandwich form. The "bread" is the inverse of the Fisher information matrix, and the "meat" is the covariance matrix of the score function. The former is estimated using the inverse of the observed Fisher information matrix. The latter is estimated using a parametric bootstrap.


copCAR returns an object of class "copCAR", which is a list containing the following components:


(if confint = "bootstrap") the bootstrap sample.


the matched call.


a named vector of parameter estimates.


the value of confint supplied in the function call.


a list containing the names and values of the control parameters.


the integer code returned by optim.


(if confint = "asymptotic") the estimate of the asymptotic covariance matrix of the parameter estimator.


the data argument.


the family object used.


the fitted mean values, obtained by transforming the linear predictors by the inverse of the link function.


the formula supplied.


the linear fit on link scale.


A character string giving any additional information returned by the optimizer, or NULL.


the method (CE, CML, or DT) used for inference.


if requested (the default), the model frame.


the number of model parameters.


the offset vector used.


the response residuals, i.e., the outcomes minus the fitted values.


the terms object used.


the value of the objective function at its minimum.


if requested, the model matrix.


(where relevant) a record of the levels of the factors used in fitting.


if requested (the default), the response vector used.


Besag, J., York, J., and MolliĆ©, A. (1991) Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20.

Denuit, M. and Lambert, P. (2005) Constraints on concordance measures in bivariate discrete data. Journal of Multivariate Analysis, 93, 40–57.

Ferguson, T. (1967) Mathematical statistics: a decision theoretic approach, New York: Academic Press.

Flegal, J., Haran, M., and Jones, G. (2008) Markov Chain Monte Carlo: can we trust the third significant figure? Statistical Science, 23(2), 250–260.

Godambe, V. (1960) An optimum property of regular maximum likelihood estimation. The Annals of Mathmatical Statistics, 31(4), 1208–1211.

Hughes, J. (2015) copCAR: A flexible regression model for areal data. Journal of Computational and Graphical Statistics, 24(3), 733–755.

Kazianka, H. and Pilz, J. (2010) Copula-based geostatistical modeling of continuous and discrete data including covariates. Stochastic Environmental Research and Risk Assessment, 24(5), 661–673.

Madsen, L. (2009) Maximum likelihood estimation of regression parameters with spatially dependent discrete data. Journal of Agricultural, Biological, and Environmental Statistics, 14(4), 375–391.

Varin, C. (2008) On composite marginal likelihoods. Advances in Statistical Analysis, 92(1), 1–28.


## Not run: 
# Simulate data and fit copCAR to them.

# Use the 20 x 20 square lattice as the underlying graph.

m = 20
A = adjacency.matrix(m)

# Create a design matrix by assigning coordinates to each vertex
# such that the coordinates are restricted to the unit square.

x = rep(0:(m - 1) / (m - 1), times = m)
y = rep(0:(m - 1) / (m - 1), each = m)
X = cbind(x, y)

# Set the dependence parameter, regression coefficients, and dispersion parameter.

rho = 0.995      # strong dependence
beta = c(1, 1)   # the mean surface increases in the direction of (1, 1)
theta = 2        # dispersion parameter

# Simulate negative binomial data from the model.

z = rcopCAR(rho, beta, X, A, family = negbinomial(theta))

# Fit the copCAR model using the continous extension, and compute 95% (default)
# asymptotic confidence intervals. Give theta the initial value of 1. Use m equal to 100.

fit.ce = copCAR(z ~ X - 1, A = A, family = negbinomial(1), method = "CE", confint = "asymptotic",
                control = list(m = 100))

# Fit the copCAR model using the DT approximation, and compute 90% confidence
# intervals. Bootstrap the intervals, based on a bootstrap sample of size 100.
# Do the bootstrap in parallel, using ten nodes.

fit.dt = copCAR(z ~ X - 1, A = A, family = negbinomial(1), method = "DT", confint = "bootstrap",
                control = list(bootit = 100, nodes = 10))
summary(fit.dt, alpha = 0.9)

# Fit the copCAR model using the composite marginal likelihood approach.
# Do not compute confidence intervals.

fit.cml = copCAR(z ~ X - 1, A = A, family = negbinomial(1), method = "CML", confint = "none")

## End(Not run)

[Package copCAR version 2.0-4 Index]