copCAR (version 2.0-4)

copCAR: Fit copCAR model to discrete areal data.

Description

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

Usage

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())

Arguments

formula

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".

family

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.

data

an optional data frame, list or environment (or object coercible by as.data.frame 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.

offset

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.

A

the symmetric binary adjacency matrix for the underlying graph.

method

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

confint

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).

model

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

x

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

y

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

verbose

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

control

a list of parameters for controlling the fitting process.

bootit

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.

m

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

rho.max

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.

epsilon

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.

parallel

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

type

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

nodes

the number of slave nodes to create.

Value

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

boot.sample

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

call

the matched call.

coefficients

a named vector of parameter estimates.

confint

the value of confint supplied in the function call.

control

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

convergence

the integer code returned by optim.

cov.hat

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

data

the data argument.

family

the family object used.

fitted.values

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

formula

the formula supplied.

linear.predictors

the linear fit on link scale.

message

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

method

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

model

if requested (the default), the model frame.

npar

the number of model parameters.

offset

the offset vector used.

residuals

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

terms

the terms object used.

value

the value of the objective function at its minimum.

x

if requested, the model matrix.

xlevels

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

y

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

Details

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<U+00E9> (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.

References

Besag, J., York, J., and Molli<U+00E9>, 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.

Examples

# 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))
summary(fit.ce)

# 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")
summary(fit.cml)
# }