Instantly share code, notes, and snippets.

@dhicks /mwe.R
Created Jul 10, 2018

Embed
What would you like to do?
Spatial Durbin Model
library(rstan)
## options(mc.cores = parallel::detectCores())
## Data should be here: <https://drive.google.com/open?id=1BUV8mChrZ0UkyW2G4EHB5wc1taUpMZUM>
load('mwe.Rdata')
## dataf: The data; response is `w_use`
## W_3nn: row-normalized spatial weights matrix, based on 3 nearest neighbors, in triplet sparse form
## W_dnn: row-normalized spatial weights matrix, distance-based, in triplet sparse form
parsed = stanc(file = 'sdm.stan', verbose = TRUE)
compiled = stan_model(stanc_ret = parsed, verbose = FALSE)
tictoc::tic()
## ~300 sec on my laptop
samples = sampling(compiled, data = list(N = nrow(dataf),
p = ncol(dataf) - 1,
y = dataf$w_use,
X = dataf[-1],
# W = as.matrix(W_3nn)
Wn = length(W_3nn@i),
Wi = W_3nn@i+1,
Wj = W_3nn@j+1,
Ww = W_3nn@x
),
chains = 2, iter = 1000,
control = list(adapt_delta = 0.8))
tictoc::toc()
library(spdep)
freq = lagsarlm(w_use ~ hispanicP + densityE, data = dataf,
type = 'Durbin',
listw = mat2listw(W_3nn))
## Compare estimates
summary(samples, pars = c('rho', 'alpha', 'beta', 'gamma'))$summary
freq
// This gist contains Stan code for a Spatial Durbin Model
// The code is a bit slow, but appears to converge well enough and yields the same point estimates as spdep::lagsarlm(type = 'Durbin')
// As written, the spatial weights matrix W is passed as a triple sparse matrix;
// in R, use the `dgTMatrix` class from the `Matrix` package for this.
data {
int<lower=0> N; // number of observations / locations
int<lower=0> p; // number of covariates
matrix[N, p] X; // covariate matrix, w/o intercept
vector<lower=0>[N] y; // response
// matrix[N, N] W; // spatial weights matrix (non-sparse version)
int Wn; // num. non-zero entries in spatial weights
int Wi[Wn]; // i index of spatial weights
int Wj[Wn]; // j index of spatial weights
vector[Wn] Ww; // spatial weight
}
parameters {
real<lower = -1, upper = 1> rho; // spatial autoregression coef.
real alpha; // intercept
vector[p] beta; // covariate effects
vector[p] gamma; // lagged covariate effects
real<lower = 0> sigma; // sd of noise
}
transformed parameters {
vector[N] y_hat;
// Non-sparse version
// y_hat = rho*W*y + alpha*rep_vector(1, N) + X*beta + W*X*gamma;
// Sparse version
{
vector[N] y_hat_local;
vector[N] y_hat_lags;
y_hat_local = alpha*rep_vector(1, N) + X*beta;
y_hat_lags = rep_vector(0, N);
for (w_idx in 1:Wn) {
y_hat_lags[Wi[w_idx]] += rho*Ww[w_idx]*y[Wj[w_idx]];
for (col_idx in 1:p) {
y_hat_lags[Wi[w_idx]] += Ww[w_idx]*X[Wj[w_idx], col_idx]*gamma[col_idx];
}
}
y_hat = y_hat_local + y_hat_lags;
}
}
model {
y ~ normal(y_hat, sigma);
// priors
alpha ~ normal(0, 3);
beta ~ normal(0, 3);
gamma ~ normal(0, 3);
sigma ~ cauchy(0, 3);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment