The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Stan based functions to estimate CAR-MM models. These models allow to estimate Generalised Linear Models with CAR (conditional autoregressive) spatial random effects for spatially and temporally misaligned data, provided a suitable Multiple Membership matrix.
The main references are:
# CRAN installation
install.packages("CARME")
# Or the development version from GitHub:
# install.packages("devtools")
::install_github("markgrama/CARME") devtools
We start by setting the number of membership we wish to use. In this case we take the adjacency matrix from the South East London (SEL) dataset used in Gramatica et al. (2023). It consists of \(m = 153\) memberships and \(n = 152\) areas.
set.seed(455)
#---- Load data
data(W_sel)
## Number of areas
<- nrow(W_sel)
n ## Number of memberships
<- 153
m #---- Simulate covariates
<- cbind(rnorm(nrow(W_sel)), rnorm(nrow(W_sel)))
X ## Min-max normalisation
<- apply(X, 2, function(x) (x - min(x))/diff(range(x))) X_cent
Then, we simulate a \((153 \times 15)2\) Multiple Membership matrix, together with the observed outcomes generated using a poisson distribution and a CAR prior.
#---- Simulate MM matrix
<- c(.5, .35, .15) # Weight of each neighbours orders
w_ord <- length(w_ord) - 1 # Order of neighbours to include
ord <- sim_MM_matrix(
H_sel_sim W = W_sel, m = m, ord = ord, w_ord = w_ord, id_vec = rep(1, nrow(W_sel))
)
#---- Simulate outcomes
## Linear term parameters
<- -.5 # Intercept
gamma <- c(1, .5) # Covariates coefficients
beta ## CAR random effects
<- sim_car(W = W_sel, alpha = .9, tau = 5)
phi_car # Areal log relative risks
<- X_cent %*% beta + phi_car
l_RR ## Membership log relative risks
<- as.numeric(apply(H_sel_sim, 1, function(x) x %*% l_RR))
l_RR_mm ## Expected rates
<- rpois(m, lambda = 20)
exp_rates ## Outcomes
<- rpois(m, lambda = exp_rates*exp(l_RR_mm))
y #---- Create dataset for stan function
<- list(
d_sel # Number of areas
n = nrow(W_sel),
# Covariates
k = ncol(X_cent),
X_cov = X_cent,
# Adjacency
W_n = sum(W_sel) / 2,
# Number of neighbour pairs
W = W_sel,
# Memberships
m = nrow(H_sel_sim),
H = H_sel_sim,
# Outcomes
y = y,
log_offset = log(exp_rates),
# Prior parameters
## Intercept (mean and sd of normal prior)
mu_gamma = 0, sigma_gamma = 1,
## Covariates (mean and sd of normal prior)
mu_beta = 0, sigma_beta = 1,
## Marginal precision gamma prior
tau_shape = 2,
tau_rate = 0.2
)
Finally, we can run the model:
#---- HMC parameters
<- 1E4
niter <- 4
nchains #---- Stan sampling
<- car_mm(
fit d_list = d_sel,
# arguments passed to sampling
iter = niter, chains = nchains, refresh = 500,
control = list(adapt_delta = .99, max_treedepth = 15)
)
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.