Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feature/car #322

Open
wants to merge 69 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
69 commits
Select commit Hold shift click to select a range
33ea96b
adding stan_CARbym; print/summary methods work
imadmali Jul 10, 2017
6feee4f
get important post-estimation methods/functions working
imadmali Jul 11, 2017
6ad17ba
adding stan_bym (inla version)
imadmali Jul 13, 2017
c2e2212
adding stan_besag (icar)
imadmali Jul 13, 2017
5b3fda8
fixing conditional in stanreg for car models
imadmali Jul 13, 2017
55fb6c6
expanding prior options for CAR models
imadmali Jul 21, 2017
6247da9
changes to spatial.fit file
imadmali Jul 25, 2017
f6f853d
Merge branch 'master' into feature/car
imadmali Aug 4, 2017
e3ea475
adding linkinv chunks to spatial.stan
imadmali Aug 4, 2017
4c0f42d
using chunk functions for meanPPD
imadmali Aug 4, 2017
2070bb4
demeaning, QR decom, and posterior methods work
imadmali Aug 6, 2017
fcc8d71
get prior summary working for CAR models
imadmali Aug 7, 2017
34f03e6
adding documentation for CAR models
imadmali Aug 7, 2017
1869c9c
adding stan_besag tests; CAR models return 'terms'
imadmali Aug 7, 2017
44f7e4a
adding vb methods for CAR models; optimizing currently fails
imadmali Aug 8, 2017
8746893
adding tests for stan_bym
imadmali Aug 8, 2017
c329a6a
use sparse matrix ops for 'scaling_factor'
imadmali Aug 8, 2017
0ebd45c
adding light documentation for stan_spatial.fit
imadmali Aug 8, 2017
9859ddb
omit optimizing from supported algorithms in CAR models
imadmali Aug 8, 2017
7c36fef
use chunks stan code; use prior_aux
imadmali Aug 9, 2017
fbdeeeb
adding more likelihoods to CAR models
imadmali Aug 9, 2017
11c132e
refactor stan chunks for spatial.stan
imadmali Aug 10, 2017
aeb898b
draft of CAR vignette; fixing binomial lpmf bug in spatial.stan
imadmali Aug 11, 2017
81cf7f6
updating documentation/vignette for CAR models
imadmali Aug 13, 2017
b65253c
adding 2nd order spatial random walks
imadmali Aug 13, 2017
a039e2a
Q should be squared if 'order == 2'
imadmali Aug 13, 2017
40b248b
update stan_besag tests with inla comparisons
imadmali Aug 13, 2017
b846069
changing stan_bym to stan_bym2 (INLA convention)
imadmali Aug 13, 2017
9aafe0b
adding BYM model; changing spatial parameter names
imadmali Aug 14, 2017
3582ca6
change parameterization in spatial.stan from prec to sd
imadmali Aug 14, 2017
9cc2a37
convert spatial weight matrix to sparse if necessary
imadmali Aug 14, 2017
87ec6a9
prior must be on theta_raw when model_type = 2
imadmali Aug 15, 2017
4ea4674
remove INLA dependency for bym2 model
imadmali Aug 30, 2017
41555ff
updating CAR spatial tests
imadmali Aug 30, 2017
31e8d1c
update CAR vignette
imadmali Aug 30, 2017
c3a55d0
changing names for spatial priors
imadmali Aug 31, 2017
2a9001a
updates to besag unit tests
imadmali Jan 12, 2018
557179b
fixing conflicts in feature/car
imadmali Jan 12, 2018
ee55e48
getting spatial models to compile
imadmali Jan 21, 2018
2da6c77
get posterior predict to work with CAR models (hacky)
imadmali Jan 22, 2018
7ab7a9a
updating spatial docs
imadmali Jan 30, 2018
a641852
return error when using loo on CAR models (tentative)
imadmali Sep 5, 2018
5a1173e
drop loo tests
imadmali Sep 5, 2018
a12bb9a
resolving merge conflicts in feature/car
imadmali Sep 6, 2018
4824fd3
modified namespace
imadmali Sep 7, 2018
02aa765
Merge branch 'master' into feature/car
imadmali Sep 7, 2018
3f52c0c
spatial.Rmd: fix encoding
bgoodri Sep 7, 2018
2876ac3
update from master
imadmali Sep 7, 2018
a129012
Merge branch 'master' into feature/car
imadmali Sep 7, 2018
b498fc0
making CAR family func consistent with other stan models
imadmali Sep 9, 2018
19a809d
fixing awkward family/link pairs; neg_bin(link=identity) still broken
imadmali Sep 10, 2018
a564548
fix count family/link issues in spatial.stan
imadmali Sep 12, 2018
8fe3ffd
cleanup loo
imadmali Oct 21, 2018
1556945
revise CAR tests
imadmali Oct 22, 2018
4ab4d1f
use sparseMatrix when dealing with CAR adjacency matrix
imadmali Oct 22, 2018
9b5bedc
simplify spatial tests
imadmali Oct 22, 2018
e65b2c7
fixing unit tests
imadmali Oct 23, 2018
04fc47a
Merge branch 'master' into feature/car
imadmali Oct 23, 2018
abd0a72
adding is_car instead of is(object, 'car')
imadmali Oct 23, 2018
557958c
revising spatial model vignette
imadmali Oct 28, 2018
05eed19
delete src/include
bgoodri Oct 28, 2018
6efcfe4
rm src/Modules*
bgoodri Oct 28, 2018
9eddcc2
dropping INLA requirement
imadmali Oct 28, 2018
e0dfbcc
fix family/link bounds in spatial.stan
imadmali Oct 28, 2018
ba88f97
Merge branch 'feature/car' of https://github.com/stan-dev/rstanarm in…
bgoodri Oct 28, 2018
5b08936
refactoring in-line declarations and likelihood contributions
imadmali Oct 28, 2018
8eb3c7f
NAMESPACE: update
bgoodri Oct 28, 2018
b3d968b
Merge branch 'feature/car' of https://github.com/stan-dev/rstanarm in…
bgoodri Oct 28, 2018
c809d14
Merge branch 'master' into feature/car
bgoodri Jul 31, 2020
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
adding 2nd order spatial random walks
  • Loading branch information
imadmali committed Aug 13, 2017
commit b65253c4c9ba7683e713ddd616af8caafe872056
3 changes: 1 addition & 2 deletions R/stan_besag.R
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 101,7 @@ stan_besag <- function(formula,
data,
trials = NULL,
W,
order = c(1,2),
order = 1,
...,
prior = normal(), prior_intercept = normal(),
prior_tau = normal(), prior_aux = NULL,
Expand All @@ -114,7 114,6 @@ stan_besag <- function(formula,
stop(paste("Please install and load the INLA package before using", stan_function))
mc <- match.call(expand.dots = FALSE)
algorithm <- match.arg(algorithm)
order <- match.arg(order)
family <- validate_family(family)
mf <- model.frame(mc, data)
mt <- terms(formula, data = data)
Expand Down
3 changes: 1 addition & 2 deletions R/stan_bym.R
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 109,7 @@ stan_bym <- function(formula,
data,
trials = NULL,
W,
order = c(1,2),
order = 1,
...,
prior = normal(), prior_intercept = normal(),
prior_tau = normal(), prior_rho = beta(0.5,0.5), prior_aux = NULL,
Expand All @@ -122,7 122,6 @@ stan_bym <- function(formula,
stop(paste("Please install the INLA package before using", stan_function))
mc <- match.call(expand.dots = FALSE)
algorithm <- match.arg(algorithm)
order <- match.arg(order)
family <- validate_family(family)
mf <- model.frame(mc, data)
mt <- terms(formula, data = data)
Expand Down
64 changes: 42 additions & 22 deletions R/stan_spatial.fit.R
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 87,9 @@ stan_spatial.fit <- function(x, y, w,
else if (stan_function == "stan_bym")
model_type <- 2

if (!(order %in% c(1,2)))
stop("Argument 'order' must be 1 or 2.")

x_stuff <- center_x(x, sparse = FALSE)
for (i in names(x_stuff)) # xtemp, xbar, has_intercept
assign(i, x_stuff[[i]])
Expand Down Expand Up @@ -228,33 231,29 @@ stan_spatial.fit <- function(x, y, w,
global_prior_scale_for_intercept,
num_normals = if(prior_dist == 7) as.integer(prior_df) else integer(0),
prior_dist_rho = prior_dist_for_rho)

# create scaling_factor a la Dan Simpson
create_scaling_factor <- function(dat) {
edges <- dat$edges
#Build the adjacency matrix
adj.matrix <- Matrix::sparseMatrix(i=edges[,1],j=edges[,2],x=1,symmetric=TRUE)
#The ICAR precision matrix (note! This is singular)
Q <- Matrix::Diagonal(dat$N, Matrix::rowSums(adj.matrix)) - adj.matrix

#Add a small jitter to the diagonal for numerical stability (optional but recommended)
Q_pert <- Q Matrix::Diagonal(dat$N) * max(Matrix::diag(Q)) * sqrt(.Machine$double.eps)

# Compute the diagonal elements of the covariance matrix subject to the
# constraint that the entries of the ICAR sum to zero.
# See the function help for further details.
Q_inv <- inla.qinv(Q_pert, constr=list(A = matrix(1,1,dat$N),e=0))

# Compute the geometric mean of the variances, which are on the diagonal of Q.inv
scaling_factor <- exp(mean(log(Matrix::diag(Q_inv))))
return(scaling_factor)
}

if (stan_function == "stan_bym")
standata$scaling_factor <- create_scaling_factor(standata)
else
standata$scaling_factor <- 0


standata$order <- order
if (order == 2) {
Q <- diag(rowSums(W) - W)
Q <- Q %*% Q
sparse_stuff <- rstan::extract_sparse_parts(Q)
standata$Q_n <- as.array(length(sparse_stuff$w), dim = 1)
standata$w <- sparse_stuff$w
standata$v <- sparse_stuff$v
standata$u <- sparse_stuff$u
}
if (order == 1) {
standata$Q_n <- array(0, dim = c(0))
standata$w <- array(0, dim = c(0))
standata$v <- array(0, dim = c(0))
standata$u <- array(0, dim = c(0))
}

pars <- c(if (has_intercept) "alpha", "beta", if(model_type == 2) "rho", "tau", if(has_aux) "aux",
"mean_PPD", "psi")

Expand Down Expand Up @@ -339,6 338,27 @@ stan_spatial.fit <- function(x, y, w,
}
}

# create scaling_factor a la Dan Simpson
create_scaling_factor <- function(dat) {
edges <- dat$edges
#Build the adjacency matrix
adj.matrix <- Matrix::sparseMatrix(i=edges[,1],j=edges[,2],x=1,symmetric=TRUE)
#The ICAR precision matrix (note! This is singular)
Q <- Matrix::Diagonal(dat$N, Matrix::rowSums(adj.matrix)) - adj.matrix

#Add a small jitter to the diagonal for numerical stability (optional but recommended)
Q_pert <- Q Matrix::Diagonal(dat$N) * max(Matrix::diag(Q)) * sqrt(.Machine$double.eps)

# Compute the diagonal elements of the covariance matrix subject to the
# constraint that the entries of the ICAR sum to zero.
# See the function help for further details.
Q_inv <- inla.qinv(Q_pert, constr=list(A = matrix(1,1,dat$N),e=0))

# Compute the geometric mean of the variances, which are on the diagonal of Q.inv
scaling_factor <- exp(mean(log(Matrix::diag(Q_inv))))
return(scaling_factor)
}

# Summarize spatial prior

summarize_spatial_prior <- function(user_prior,
Expand Down
31 changes: 21 additions & 10 deletions exec/spatial.stan
Original file line number Diff line number Diff line change
Expand Up @@ -12,20 12,28 @@ data {
int<lower=0> K; // number of predictors (inc intercept)
matrix[N,K] X; // model matrix
int<lower=1,upper=5> family;
int link;
int<lower=0,upper=1> is_continuous;
int<lower=0,upper=1> has_aux;
int<lower=1,upper=2> model_type; // Besag = 1; BYM = 2
int<lower=1,upper=2> model_type; // 0 = 1D RW; Besag = 1; BYM = 2
int<lower=0,upper=1> has_intercept;
vector[K] xbar;
int<lower=0> trials[family == 4 ? N : 0]; // binomial trials (0 1d array if not applicable)
int y_int[is_continuous == 1 ? 0 : N]; // outcome
vector[is_continuous == 1 ? N : 0] y_real; // outcome
int link;
int<lower=0> trials[family == 4 ? N : 0];
int y_int[is_continuous == 1 ? 0 : N];
vector[is_continuous == 1 ? N : 0] y_real;
// pairwise difference version of CAR
int E_n; // number of adjacency pairs
int edges[E_n, 2]; // adjacency pairs
// matrix version of CAR
int<lower=1,upper=2> order;
int<lower=0> Q_n[order == 2];
vector[order == 2 ? Q_n[1] : 0] w;
int v[order == 2 ? Q_n[1] : 0];
int u[order == 2 ? N 1 : 0];
// prior stuff
int<lower=0,upper=1> prior_dist_rho;
real<lower=0> shape1_rho; // priors
real<lower=0> shape2_rho; // priors
real<lower=0> shape1_rho;
real<lower=0> shape2_rho;
real scaling_factor;
int<lower=0> prior_dist_for_intercept;
int<lower=0> prior_dist;
Expand All @@ -39,8 47,8 @@ data {
real prior_mean_tau;
real<lower=0> prior_scale_tau;
real<lower=0> prior_df_tau;
real<lower=0> global_prior_df; // for hs priors only
real<lower=0> global_prior_scale; // for hs priors only
real<lower=0> global_prior_df;
real<lower=0> global_prior_scale;
int<lower=2> num_normals[prior_dist == 7 ? K : 0];
int<lower=0,upper=3> prior_dist_for_aux;
real<lower=0> prior_mean_for_aux;
Expand Down Expand Up @@ -114,7 122,10 @@ model {
target = GammaReg(y_real, eta, aux, link, sum_log_y);
}
// prior on spatial parameter vector (GMRF)
target = -0.5 * dot_self(phi[edges[,1]] - phi[edges[,2]]);
if (order == 1)
target = -0.5 * dot_self(phi[edges[,1]] - phi[edges[,2]]);
else if (order == 2)
target = -0.5 * dot_product(phi, csr_matrix_times_vector(N, N, w, v, u, phi));
// priors on coefficients
#include "priors.stan"
// model specific priors
Expand Down