revising spatial model vignette

stan-dev · imadmali · Jul 10, 2017 · Jul 11, 2017 · Jul 13, 2017 · Jul 13, 2017
commit 557958c44c799d17f15f84ad440906b10744adb5
diff --git a/R/log_lik.R b/R/log_lik.R
@@ -271,7  271,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL,
         trials <- 1
         if (is.factor(y))
           y <- fac2bin(y)
-        if (!is(object, "car"))
         if (!is_car(object))
           stopifnot(all(y %in% c(0, 1)))
         else
           trials <- object$trials
@@ -359,7  359,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL,
     data <- cbind(data, as.matrix(z)[1:NROW(x),, drop = FALSE])
     draws$beta <- cbind(draws$beta, b)
   }
-  if (is(object, "car")) {
   if (is_car(object)) {
     psi_indx <- grep("^psi\\[[[:digit:]] \\]", colnames(stanmat))
     psi <- stanmat[, psi_indx, drop = FALSE]
     data$psi <- t(psi)

diff --git a/R/stan_besag.R b/R/stan_besag.R
@@ -151,7  151,7 @@ stan_besag <- function(formula,
                stan_function = stan_function)
 
   if (family$family == "binomial") {
-    fit$family <- binomial(link = "logit")
     # fit$family <- binomial(link = "logit")
     fit$trials <- trials
   }
   out <- stanreg(fit)

diff --git a/R/stanreg.R b/R/stanreg.R
@@ -31,7  31,7 @@ stanreg <- function(object) {
   nobs <- NROW(y)
   ynames <- if (is.matrix(y)) rownames(y) else names(y)
 
-  # is_car <- object$stan_function %in% c("stan_besag", "stan_bym", "stan_bym2")
   is_car <- object$stan_function %in% c("stan_besag", "stan_bym", "stan_bym2")
   is_betareg <- is.beta(family$family)
   if (is_betareg) { 
     family_phi <- object$family_phi  # pull out phi family/link
@@ -92,7  92,7 @@ stanreg <- function(object) {
 
   # linear predictor, fitted values
   eta <- linear_predictor(coefs, x, object$offset)
-  if (is_car(object)) {
   if (is_car) {
     psi_indx <- grep("psi", colnames(as.matrix(object$stanfit)))
     psi <- as.matrix(object$stanfit)[,psi_indx]
     psi <- unname(colMeans(psi))
@@ -156,7  156,7 @@ stanreg <- function(object) {
     out$eta_z <- eta_z
     out$phi <- phi
   }
-  if (is_car(object)) {
   if (is_car) {
     out$psi <- psi
     out$trials <- object$trials
   }

diff --git a/vignettes/spatial.Rmd b/vignettes/spatial.Rmd
@@ -41,27  41,31 @@ The linear predictor takes the following form,
 $$
 \boldsymbol{\eta} = \alpha   \mathbf{X}\boldsymbol{\beta}   \boldsymbol{\psi}
 $$
-where $\alpha$ is the intercept, $\mathbf{X}$ is an $N$-by-$K$ matrix of predictors, $\boldsymbol{\beta}$ is a $K$-dimensional vector of regression coefficients, and $\boldsymbol{\psi}$ is a $N$-dimensional vector representing the spatial effect (and $N$ denotes the number of spatial units). The construction of $\boldsymbol{\psi}$ depends on the model, which is discussed in the relevant sections below.
 where $\alpha$ is the intercept, $\mathbf{X}$ is an $N$-by-$K$ matrix of predictors ($N$ being the number of observations and $K$ being the number of predictors), $\boldsymbol{\beta}$ is a $K$-dimensional vector of regression coefficients, and $\boldsymbol{\psi}$ is a $N$-dimensional vector representing the spatial effect. The construction of $\boldsymbol{\psi}$ depends on the model, which is discussed in the relevant sections below.
 
 Depending on the choice of likelihood there may or may not be an additional auxiliary parameter $\gamma$ in the model (e.g. in a Gaussian likelihood this would be the variation of the data). With all these components, for some probability density/mass function $f$, we can state the general form of the likelihood as,
 
-Depending on the choice of likelihood there may or may not be an additional auxiliary parameter $\gamma$ in the model (e.g. in a Gaussian likelihood this would be the variation of the data). Thus, for some probability density/mass function $f$ we can state the general form of the likelihood as,
 $$
 \mathcal{L}(\alpha, \boldsymbol{\beta}, \gamma | \mathbf{y}) = \prod_{i=1}^N f(y_i | \alpha, \boldsymbol{\beta}, \gamma )
 $$
 
 ## GMRF Hierarchical Component
 
 CAR models require that you define the spatial component as a Gaussian Markov Random Field (GMRF). The random vector $\boldsymbol{\phi}$ is a GMRF with respect to the graph $\mathcal{G} = (\mathcal{V} = \{1,\ldots,n\},\mathcal{E})$ with mean vector $\boldsymbol{\mu}$ and precision matrix $\mathbf{W}$ if its probability density function takes the precision form of the multivariate normal distribution,
 
 $$
 f(\boldsymbol{\phi} | \boldsymbol{\mu}) = (2\pi)^{-n/2} |\mathbf{W}|^{1/2}
 \exp\bigg( -\frac{1}{2}(\boldsymbol{\phi-\mu})^{\top}\mathbf{W}(\boldsymbol{\phi-\mu}) \bigg)
 $$
 
 and $W_{i,j} \neq 0$ when $\{i,j\}\in\mathcal{E}$ for all $i \neq j$. Here, $\mathcal{V}$ refers to the verticies on the graph (i.e. the spatial units) and $\mathcal{E}$ refers to the edges on the graph (i.e. the spatial units that are neighbors). In other words, $\mathbf{W}$ is an $N$-by-$N$ matrix  (with zeros on the diagonal) which describes spatial adjacency. 
 
 Unfortunately there is no guarantee that $\mathbf{W}$ is positive definite so $\mathbf{Q} = \mbox{diag}(\mathbf{W1}) - W$ (which is guaranteed to be positive semi-definite) is used as the precision matrix.
 
 ## Besag (ICAR) Spatial Prior
 
 The `stan_besag` modeling function fits the data to an ICAR model. This means that the spatial effect enters the linear predictor as
 
 $$
 \begin{align}
 \boldsymbol{\psi} = \boldsymbol{\phi} \\
@@ -74,7  78,8 @@ where $\tau$ is a scalar that controls the overall spatial variation and has an
 
 ## BYM
 
-The ICAR model is limited in that it only accounts for spatial variation among the spatial units. Thus, the random variation is picked up by the spatial variation which results in misleading parameter estimates and invalid inferences. The `stan_bym` model explains spatial variation as the sum of a structured (spatial) component $\boldsymbol{\phi}$ and an unstructured (random) component $\boldsymbol{\theta}$. Therefore, the spatial effect takes the following form,
 The ICAR model is limited in that it only accounts for spatial variation among the spatial units. Thus, the random variation is picked up by the spatial variation which may result in misleading parameter estimates and invalid inferences. The `stan_bym` model explains spatial variation as the sum of a structured (spatial) component $\boldsymbol{\phi}$ and an unstructured (random) component $\boldsymbol{\theta}$. Therefore, the spatial effect takes the following form,
 
 $$
 \begin{align}
 \boldsymbol{\psi} &= \rho\boldsymbol{\phi}   \tau\boldsymbol{\theta} \\
@@ -83,11  88,13 @@ f(\boldsymbol{\phi} | \boldsymbol{\mu}) &= (2\pi)^{-n/2} |\mathbf{Q}|^{1/2}
 f(\theta_i) &= (2\pi)^{-1/2}\exp{\bigg( -\frac{\theta_i^2}{2} \bigg)}
 \end{align}
 $$
 
 Note that the unstructured effect $\boldsymbol{\theta}$ is distributed standard normal.
 
 ## BYM2 (Variant of the BYM Spatial Prior)
 
 The `stan_bym2` modeling function fits the data to a variant of the BYM model where the spatial effect enters as a convolution of the structured (spatial) effect and the unstructured (random) effect,
 
 $$
 \begin{align}
 \boldsymbol{\psi} &= \tau(\boldsymbol{\theta}\sqrt{1-\rho}   \boldsymbol{\phi}\sqrt{\rho}) \\
@@ -97,9  104,9 @@ f(\theta_i) &= (2\pi)^{-1/2}\exp{\bigg( -\frac{\theta_i^2}{2} \bigg)}
 \end{align}
 $$
 
-As in the BYM model $\boldsymbol{\theta}$ is distributed standard normal. However, the parameter $\rho$ is on the unit interval and is interpreted as the proportion of spatial variation that is contributed to overall variation, and $\tau$ explains the overall (convolved) variation.
 As in the BYM model $\boldsymbol{\theta}$ is distributed standard normal. However, now the parameter $\rho$ is on the unit interval and is interpreted as the proportion of spatial variation that is contributed to overall variation, and $\tau$ explains the overall (convolved) variation.
 
-Priors on $\rho$ should be chosen wisely as $\rho=0$ reduces to a model that not account for spatial variation and $\rho=1$ reduces to the ICAR model, which does not account for random variation among the spatial units. 
 Priors on $\rho$ should be chosen carefully as $\rho=0$ reduces to a model that not account for spatial variation and $\rho=1$ reduces to the ICAR model, which does not account for random variation among the spatial units. 
 
 ## Posterior
 
@@ -193,19  200,28 @@ spplot(grid_sim, "y_pred", at = var_range_y_pred, main = expression(y[pred]),
 ```
 
 Alternatively we can look at the conventional one-dimensional posterior predictive check with the `pp_check` function.
 
 ```{r ppcheck-1d, fig.align='center', fig.height=8}
 pp_check(fit_besag)
 ```
 
-In order to compare the predictive performance between the models we need to use the [loo](http://mc-stan.org/loo) package. Currently, however, this is not supported for CAR spatial models.
 In order to compare the predictive performance between the models we need to use the [loo](http://mc-stan.org/loo) package. Looking at the results below we can confirm that the correct model outperforms the incorrect model in terms of predictive accuracy.
 
 ```{r loo}
 library(loo)
 loo_correct <- loo(fit_besag)
 loo_incorrect <- loo(fit_besag_bad)
 compare(loo_correct, loo_incorrect)
 ```
 
 ## Smoothing the Spatial Random Walk
 
 In some cases modeling the GMRF spatial component with the precision matrix $\mathbf{Q}$ leads to rough spatial varation. This often occurs when dealing with spatial units on a fine lattice. Using the square of the precision matrix $\mathbf{Q}\mathbf{Q}$ allows us to smooth out the spatial variation across the spatial units using a spatial random walk of order 2. Below we use the `order = 2` argument to fit the model with smoothing of order 2.
 
 ```{r smoothing, results="hide"}
 fit_besag_smooth <- stan_besag(y ~ 1   x   z, data = spatial_data, W, order = 2,
-                               prior_intercept = normal(0,1), prior = normal(0,1), prior_rho = normal(0,1),
                                prior_intercept = normal(0,1), prior = normal(0,1),
                                prior_rho = normal(0,1),
                                family = binomial(link="logit"), trials = spatial_data$trials,
                                chains = CHAINS, cores = CORES, seed = SEED, iter = ITER)
 ```