Package 'spatemR'

Generalized Estimating Equations with Spatial Autoregressive Components

Description

'GSCIMC' estimates generalized estimating equations (GEE) incorporating spatial autoregressive (SAR) components. It extends GEE models to account for spatial dependence in the response variable.

Usage

GSCIMC(
  formula,
  family = gaussian(),
  weights = NULL,
  data,
  W,
  start = NULL,
  toler = 1e-04,
  maxit = 200,
  trace = FALSE,
  eps = 1e-06
)

Arguments

formula	A formula specifying the model structure (response ~ predictors).
family	A description of the error distribution and link function. Default is 'gaussian()'.
weights	Optional vector of prior weights. Must be positive.
data	A data frame containing the variables in the model.
W	A spatial weights matrix defining the spatial dependence structure.
start	Optional starting values for parameter estimation.
toler	Convergence tolerance for iterative optimization. Default is '1e-05'.
maxit	Maximum number of iterations for model fitting. Default is '50'.
trace	Logical; if 'TRUE', prints iteration details. Default is 'FALSE'.
eps	Minimun value for variance diference of zero

Details

The function estimates a spatially autoregressive GEE model by iteratively updating the spatial dependence parameter ('rho') and regression coefficients ('beta'). The estimation follows a quasi-likelihood approach using iterative weighted least squares (IWLS).

The function supports common GLM families ('gaussian', 'binomial', 'poisson', 'Gamma', 'inverse.gaussian') and their quasi-likelihood equivalents.

Value

A list of class '"GSCIMC"' containing:

coefficients	Estimated regression coefficients.
rho	Estimated spatial autoregressive parameter.
fitted.values	Predicted values from the model.
linear.predictors	Linear predictor values ('X * beta').
prior.weights	Weights used in estimation.
y	Observed response values.
formula	Model formula.
call	Function call used to fit the model.
data	Data used in the model.
converged	Logical indicating whether the algorithm converged.
logLik	Quasi-log-likelihood of the fitted model.
deviance	Residual deviance.
df.residual	Residual degrees of freedom.
phi	Dispersion parameter estimate.
R	Robust Variance Estimation.
CIC	Corrected Information Criterion.
RJC	Robust Jackknife Correction.

Source

https://doi.org/10.48550/arXiv.2412.00945

References

Cruz, N. A., Toloza-Delgado, J. D., & Melo, O. O. (2024). Generalized spatial autoregressive model. arXiv preprint arXiv:2412.00945.

See Also

glm, gee, spdep

Examples

library(spdep)
library(sp)
data(meuse)
sp::coordinates(meuse) <- ~x+y
W <- spdep::nb2mat(knn2nb(knearneigh(meuse, k=5)), style="W")
fit <- GSCIMC(cadmium ~ dist + elev, family=poisson(), data=meuse, W=W)
summary_SAR(fit)

---

Hurdle Model using GSCIMC

Description

This function fits a hurdle model using GSCIMC, consisting of: (1) A logit model for zero vs. non-zero responses. (2) A truncated Poisson model for positive counts.

Usage

Hurdle_GSCIMC(
  formula,
  data,
  W,
  weights = NULL,
  toler = 1e-05,
  maxit = 200,
  trace = FALSE
)

Arguments

formula	A formula specifying the model.
data	The dataset.
W	The spatial weight matrix.
weights	Optional weights.
toler	Convergence tolerance.
maxit	Maximum number of iterations.
trace	Logical. If TRUE, prints progress.

Value

A list containing the logit and Poisson-truncated models.

Examples

set.seed(123)
n <- 100
x <- rnorm(n)
y <- rpois(n, lambda = exp(0.5 * x))
y[rbinom(n, 1, 1/(1+exp(-0.5*x)))] <- 0  # Introduce zeros
W <- matrix(rbinom(n^2,1,0.2), n, n)  # Example spatial weight matrix
diag(W) <- 0
rtot <- rowSums(W)
W <- W/ifelse(rtot==0, 0.1, rtot)
model <- Hurdle_GSCIMC(y ~ x, data = data.frame(y, x), W = W)
summary_SAR(model$logit_model)
summary_SAR(model$poisson_truncated_model)

---

Print Method for Summary of GSCIMC Models

Description

This method prints a formatted summary of a 'summary.GSCIMC' object, including details of the model coefficients, rho, dispersion, and other statistics.

Usage

## S3 method for class 'summary.GSCIMC'
print(x, ...)

Arguments

x	An object of class 'summary.GSCIMC'.
...	Additional arguments (currently unused).

Value

Print a summary for the specified Generalized Spatial Autoregresive Model class.

---

Print Method for Summary of SARARgamlss Models

Description

This method prints a formatted summary of a 'summary.SARARgamlss' object, including details of the GAMLSS model, spatial parameters (rho and lambda), and Wald tests.

Usage

## S3 method for class 'summary.SARARgamlss'
print(x, ...)

Arguments

x	An object of class 'summary.SARARgamlss'.
...	Additional arguments (currently unused).

Value

Print a summary for the specified GAMLSS model.

---

Truncated Poisson Family for GLM

Description

This function defines a truncated Poisson family for use in Generalized Linear Models (GLMs), where zero values are not allowed. It modifies the Poisson likelihood by excluding zero-count observations.

Usage

ptfamily(link = "log")

Arguments

link	Character string or a link-glm object specifying the link function. Accepted values are "log", "identity", and "sqrt".

Value

An object of class "family" that can be used in glm().

Examples

set.seed(123)
y <- rpois(100, lambda = 3)
y <- y[y > 0]  # Truncate zeros
x <- rnorm(length(y))
model <- glm(y ~ x, family = ptfamily())
summary(model)

---

SARARgamlss: Spatial Autoregressive Generalized Additive Model for Location Scale (GAMLSS)

Description

This function estimates a Spatial Autoregressive Generalized Additive Model for Location Scale (SARARgamlss) using GAMLSS. The model includes both spatial dependencies and the possibility of non-parametric terms in the formulas for the mean and variance. The function supports SAR, SARAR, and SEM model types and performs the estimation through an iterative process that updates spatial dependence parameters. The variance of the spatial parameters $\hat{\rho}$ and $\hat{\lambda}$ is estimated using the inverse of the Hessian matrix from the optimization.

Usage

SARARgamlss(
  formula,
  sigma.formula = ~1,
  W1 = diag(0, nrow(data)),
  W2 = diag(0, nrow(data)),
  data,
  tol = 1e-04,
  maxiter = 20,
  type = c("SAR", "SARAR", "SEM"),
  weights = NULL
)

Arguments

formula	A formula specifying the mean structure of the model (response ~ explanatory variables).
sigma.formula	A formula specifying the variance structure of the model (default: ~1).
W1	A spatial weights matrix for the SAR term (default: identity matrix).
W2	A spatial weights matrix for the SARAR term (default: identity matrix).
data	A data.frame containing the variables used in the model.
tol	Convergence tolerance (default: 1E-4).
maxiter	Maximum number of iterations for optimization (default: 20).
type	The type of spatial model to fit: one of "SAR", "SARAR", or "SEM".
weights	Optional weights for the observations (default: NULL).

Value

A fitted GAMLSS model object with spatial autoregressive terms. The model object also includes the variance of the spatial parameters $\hat{\rho}$ and $\hat{\lambda}$

References

Toloza-Delgado, J. D., Melo, O. O., & Cruz, N. A. Joint spatial modeling of mean and non-homogeneous variance combining semiparametric SAR and GAMLSS models for hedonic prices. Spatial Statistics, 65, 100864 (2025) @source https://doi.org/10.1016/j.spasta.2024.100864

Examples

library(spdep)
library(gamlss)
data(oldcol)
# Create spatial weight matrices W1 and W2
W1 <- spdep::nb2mat(COL.nb, style = "W")
W2 <- W1  # In this case, assume the same spatial weights for both
# Fit a SARARgamlss model
result <- SARARgamlss(formula = CRIME ~ INC + cs(HOVAL), 
sigma.formula = ~ INC + pb(HOVAL), W1 = W1, W2 = W2,data = COL.OLD, 
tol = 1E-4,  maxiter = 20, type = "SARAR")
summary_SAR(result)
gamlss::term.plot(result$gamlss, what="mu")

---

Custom Summary Function for SARARgamlss and GSCIMC Models

Description

This function generates a summary for objects of class 'SARARgamlss' or 'GSCIMC'. It combines the summary outputs for both models, including GAMLSS model details, spatial parameters (rho and lambda), and Wald tests.

Usage

summary_SAR(object)

Arguments

object

An object of class 'SARARgamlss' or 'GSCIMC'.

Value

A list containing the summary for the specified model class.

Examples

library(spdep)
library(gamlss)
data(oldcol)
W1 <- spdep::nb2mat(COL.nb, style = "W")
W2 <- W1  # In this case, assume the same spatial weights for both
# Fit a SARARgamlss model
result_sarar <- SARARgamlss(formula = CRIME ~ INC + HOVAL, 
                            sigma.formula = ~ INC + pb(HOVAL), 
                            W1 = W1, W2 = W2, data = COL.OLD,
                            type="SAR")
summary_SAR(result_sarar)

# Example for GSCIMC model
result_GSCIMC <- GSCIMC(formula = CRIME ~ INC + HOVAL, data = COL.OLD, W = W1)
summary_SAR(result_GSCIMC)

---

Compute the Inverse Variance of the Spatial Autoregressive Parameter (rho)

Description

This function calculates the inverse of the variance of the spatial autoregressive parameter $\rho$ in a generalized spatial autoregressive (GSAR) or GEE-SAR model. The calculation is based on the quasi-likelihood derivatives with respect to $\rho$ for different exponential family distributions.

Usage

var_rho_inv(A, W, X, beta, family, y, offs = NULL, weights = NULL, phi = 1)

Arguments

A	Matrix. The spatial transformation matrix $\mathbf{A} = \mathbf{I} - \rho \mathbf{W}$, typically of class 'Matrix'.
W	Matrix. Row-standardized spatial weights matrix $\mathbf{W}$ of dimension $n \times n$.
X	Matrix. Design matrix of covariates, dimension $n \times p$.
beta	Numeric vector. Current estimates of regression coefficients $\pmb{\beta}$, length $p$.
family	GLM family object. The response distribution family (e.g., 'gaussian()', 'poisson()', 'binomial()', 'Gamma()', 'Negative Binomial()').
y	response variable
offs	Numeric vector. Optional offset vector, length $n$. Default is 0.
weights	Numeric vector. Observation weights $m_i$ (e.g., number of trials for binomial data), length $n$.
phi	Numeric. Dispersion parameter, used for 'gaussian', 'Gamma', or 'Negative Binomial' families. Default is 1.

Details

The function computes first and second derivatives of the mean $\mu_i$ with respect to $\rho$, and then applies the appropriate formula for the inverse variance based on the selected family. This generalizes the quasi-likelihood derivations for spatially correlated generalized linear models.

For binomial families with large $m_i$, it is recommended to truncate $\mu_i$ within $[1e-10, 1-1e-10]$ to avoid numerical instability.

Value

Numeric. The inverse of the variance of $\hat{\rho}$ ($\operatorname{Var}(\hat{\rho})^{-1}$).

Examples

## Not run: 
library(Matrix)
n <- 10
W <- Matrix(0,n,n)
diag(W[-1,]) <- 1
X <- matrix(rnorm(n*2), n, 2)
beta <- c(0.5, -0.2)
rho <- 0.3
A <- Diagonal(n) - rho*W
family <- binomial()
weights <- rep(1,n)
var_rho_inv(A, W, X, beta, family, weights)

## End(Not run)

---