feature/spdep

[imadmali]

I'm just starting an issue to track the feature/spdep branch which incorporates (some) spatial models from spdep in rstanarm.

Currently it looks something like this:

       ┌─ stan_lagsarlm.R ───┐
user ──┤                     ├─ stan_sp.fit.R ── spatial.stan
       └─ stan_errorsarlm.R ─┘

spdep doesn't use .fit.R files for doing the heavy lifting for their modeling functions so currently I deal with all of the models in stan_sp.fit.R (to stay somewhat consistent with the other rstanarm models). As for spatial weights, the user can include a a listw object (as required in spdep) or a symmetric matrix (usually row normalized) in the call.

There are two other models that we might also want to include,

• spdep::sacsarlm() (Spatial simultaneous autoregressive SAC model estimation)
• spdep::lagmess() (matrix exponential spatial lag model)

Stuff to implement for stan_lagsarlm.R and stan_errorsarlm.R:

• [x] check spatial weights (row normalized, convert from listw, etc.)
• [x] center data
• [x] QR decomposition
• [x] prior for spatial autocorrelation (beta)
• [x] prior for intercept
• [x] prior summary
• [x] optimize/vb
• [x] print/summary stuff (in stanreg.R)
• [x] predict/posterior_predict
• [ ] loo
• [x] examples
• [x] documentation
• [ ] vignette(s)
• [ ] tests

The response variable is basically modeled as multivariate normal with an NxN covariance matrix,

Sigma = solve( (I - rho * W) %*% t(I - rho * W) )

(where I is a diagonal matrix, W is a spatial weight matrix, and rho is spatial correlation term) and linear predictor,

mu = X * beta (times other stuff, depending on the model)

[bgoodri]

Great!

On Wed, Mar 1, 2017 at 10:19 AM, Imad Ali notifications@github.com wrote:

I'm just starting an issue to track the feature/spdep branch which incorporates (some) spatial models from spdep https://cran.r-project.org/web/packages/spdep/index.html in rstanarm.

Currently it looks something like this:

   ┌─ stan_lagsarlm.R ───┐

user ──┤ ├─ stan_sp.fit.R ── spatial.stan └─ stan_errorsarlm.R ─┘

spdep doesn't use .fit.R files for doing the heavy lifting for their modeling functions so currently I deal with all of the models in stan_sp.fit.R (to stay somewhat consistent with the other rstanarm models). As for spatial weights, the user can include a a listw object (as required in spdep) or a symmetric matrix (usually row normalized) in the call.

There are two other models that we might also want to include,

• spdep::sacsarlm() (Spatial simultaneous autoregressive SAC model estimation)
• spdep::lagmess() (matrix exponential spatial lag model)

Stuff to implement for stan_lagsarlm.R and stan_errorsarlm.R:

• check spatial weights (row normalized, convert from listw, etc.)
• center data
• QR decomposition
• prior for spatial autocorrelation (beta)
• prior for intercept
• optimize/vb
• print/summary stuff (in stanreg.R)
• predict/posterior_predict
• loo
• examples
• documentation
• vignette(s)
• tests

The response variable is basically modeled as multivariate normal with an NxN covariance matrix,

Sigma = solve( (I - rho W) %% t(I - rho * W) )

(where I is a diagonal matrix, W is a spatial weight matrix, and rho is spatial correlation term) and linear predictor,

mu = X * beta (times other stuff, depending on the model)

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/166, or mute the thread https://github.com/notifications/unsubscribe-auth/ADOrqi2GVVpwvkA3grPszZnikw5jvRNQks5rhYxmgaJpZM4MPyU_ .

[imadmali]

The only difficult thing that remains for this is to get loo working with these spatial models. More generally it's a problem of how to get loo to work on a N-dimensional vector of dependent random variables. Any suggestions @bgoodri, @jgabry, @avehtari?

[avehtari]

@imadmali if N is the number of observations and "The response variable is basically modeled as multivariate normal with an NxN covariance matrix" then it's not difficult as discussed in Bayesian Leave-One-Out Cross-Validation Approximations for Gaussian Latent Variable Models http://jmlr.org/papers/v17/14-540.html . Conditional loo can be computed by computing conditional log_lik's using 1a) equation 31 and normal_lpdf, or 1b) equations 32-33 and normal_lpdf, or 1c) equation 34 and then 2) overall loo is computed using the conditional log_lik's as described in section 3.6.1.

[imadmali]

Thanks @avehtari, and sorry for the delay in getting around to this. For equation 31 in the paper, how do I ensure that the standard deviation parameter is positive? My understanding from what I've read is that the variables are defined as follows:

sigma = standard deviation of errors (scalar) Sigma = covariance matrix of N observations from multivariate normal mu = mean vector of N observations from multivariate normal

In the case of the spatial autocorrelation model I have:

mu = solve(I - rho * W) %*% beta
Sigma = sigma * solve(I - rho * W) %*% t(solve(I - rho * W))

where rho is the spatial autocorrelation term, W is a weight matrix, and beta is the vector of coefficients on the predictors.

In some cases I get negative values of v_i (based on how it's calculated in the paper).

Also, since I don't have latent variables f_i, is it correct to calculate the point-wise log likelihood with the outcome: normal_lpdf(y | mu_i, v_i)?

[avehtari]

Negative values really close to zero or negative values far away from zero? What is the condition number of Sigma? It's possible to get small negative values if if Sigma is ill-conditioned. Maybe you should add a nugget term?

Yes, you have to compute the likelihood with outcome. Remember also that in the paper v_i is variance, but Stan normal_lpdf takes standard deviation, ie, sqrt(v_i).

Also you should not compute solve(I - rho * W) twice (compute it once and assign it to a variable)

[imadmali]

Negative values really close to zero or negative values far away from zero? What is the condition number of Sigma? It's possible to get small negative values if if Sigma is ill-conditioned. Maybe you should add a nugget term?

They're close to zero and seem to consistently be around -0.05. The condition number of Sigma for one of the samples of sigma is kappa(Sigma) = 32.58. I haven't dealt with "nugget terms" much. From what I understand it's a scalar value that I would have to add to the diagonal of the covariance matrix Sigma.

Yes, you have to compute the likelihood with outcome. Remember also that in the paper v_i is variance, but Stan normal_lpdf takes standard deviation, ie, sqrt(v_i).

Thanks, I overlooked that.

Also you should not compute solve(I - rho * W) twice (compute it once and assign it to a variable)

Yeah I literally transcribed the math here but in R/Stan I've made sure to do these sorts of computations once as you've suggested.

[imadmali]

So it seems that the issue with negative variances was because I forgot to square the standard deviation sigma when computing the covariance matrix Sigma. i.e. I was calculating,

weight_stuff = solve(I - rho * W)
Sigma = sigma * weight_stuff %*% t(weight_stuff)

instead of

weight_stuff = solve(I - rho * W)
Sigma = sigma^2 * weight_stuff %*% t(weight_stuff)