The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.
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")
devtools::install_github("markgrama/CARME")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
n <- nrow(W_sel)
## Number of memberships
m <- 153
#---- Simulate covariates
X <- cbind(rnorm(nrow(W_sel)), rnorm(nrow(W_sel)))
## Min-max normalisation
X_cent <- apply(X, 2, function(x) (x - min(x))/diff(range(x)))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
w_ord <- c(.5, .35, .15) # Weight of each neighbours orders
ord <- length(w_ord) - 1 # Order of neighbours to include
H_sel_sim <- sim_MM_matrix(
W = W_sel, m = m, ord = ord, w_ord = w_ord, id_vec = rep(1, nrow(W_sel))
)
#---- Simulate outcomes
## Linear term parameters
gamma <- -.5 # Intercept
beta <- c(1, .5) # Covariates coefficients
## CAR random effects
phi_car <- sim_car(W = W_sel, alpha = .9, tau = 5)
# Areal log relative risks
l_RR <- X_cent %*% beta + phi_car
## Membership log relative risks
l_RR_mm <- as.numeric(apply(H_sel_sim, 1, function(x) x %*% l_RR))
## Expected rates
exp_rates <- rpois(m, lambda = 20)
## Outcomes
y <- rpois(m, lambda = exp_rates*exp(l_RR_mm))
#---- Create dataset for stan function
d_sel <- list(
# 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
niter <- 1E4
nchains <- 4
#---- Stan sampling
fit <- car_mm(
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.
Health stats visible at Monitor.