Lab sessions 5 - Bayesian Statistics
May, 11 2020
Cockroaches
Recap
We are analysing data about complaints for the presence of cocroaches in 10 residential buildings. Using a Bayesian framework, we modelled the number of complaints using a Poisson model. Diagnostic analysis of the model highlights some residual variability. We decided to try to improve our estimates adding one predictor and an exposure term to the Poisson model, but the situation did not improve. The ‘Poisson’ assumption on the model variance turned out to be too restrictive to model our data. The Negative Binomial model gave better results, but examining the standardized residual plot we noticed that some of them were still very large.
A posterior predictive check considering that the data are clustered by building.
library(rstan)
library(dplyr)
library(lubridate)
library(ggplot2)
library(bayesplot)
theme_set(bayesplot::theme_default())
rstan_options(auto_write=TRUE)pest_data <- readRDS('pest_data.RDS')
str(pest_data)## 'data.frame': 120 obs. of 14 variables:
## $ mus : num 0.369 0.359 0.282 0.129 0.452 ...
## $ building_id : int 37 37 37 37 37 37 37 37 37 37 ...
## $ wk_ind : int 1 2 3 4 5 6 7 8 9 10 ...
## $ date : Date, format: "2017-01-15" "2017-02-14" ...
## $ traps : num 8 8 9 10 11 11 10 10 9 9 ...
## $ floors : num 8 8 8 8 8 8 8 8 8 8 ...
## $ sq_footage_p_floor : num 5149 5149 5149 5149 5149 ...
## $ live_in_super : num 0 0 0 0 0 0 0 0 0 0 ...
## $ monthly_average_rent: num 3847 3847 3847 3847 3847 ...
## $ average_tenant_age : num 53.9 53.9 53.9 53.9 53.9 ...
## $ age_of_building : num 47 47 47 47 47 47 47 47 47 47 ...
## $ total_sq_foot : num 41192 41192 41192 41192 41192 ...
## $ month : num 1 2 3 4 5 6 7 8 9 10 ...
## $ complaints : num 1 3 0 1 0 0 4 3 2 2 ...summary(pest_data)## mus building_id wk_ind date
## Min. :-0.1008 Min. : 5.0 Min. : 1.00 Min. :2017-01-15
## 1st Qu.: 0.5279 1st Qu.:26.0 1st Qu.: 3.75 1st Qu.:2017-04-07
## Median : 1.0422 Median :46.0 Median : 6.50 Median :2017-06-29
## Mean : 1.1116 Mean :49.6 Mean : 6.50 Mean :2017-06-29
## 3rd Qu.: 1.6837 3rd Qu.:70.0 3rd Qu.: 9.25 3rd Qu.:2017-09-19
## Max. : 2.8030 Max. :98.0 Max. :12.00 Max. :2017-12-11
## traps floors sq_footage_p_floor live_in_super
## Min. : 1.000 Min. : 4.0 Min. :4186 Min. :0.0
## 1st Qu.: 6.000 1st Qu.: 8.0 1st Qu.:4770 1st Qu.:0.0
## Median : 7.000 Median :10.0 Median :5097 Median :0.0
## Mean : 7.033 Mean : 9.9 Mean :4991 Mean :0.3
## 3rd Qu.: 8.000 3rd Qu.:13.0 3rd Qu.:5206 3rd Qu.:1.0
## Max. :11.000 Max. :15.0 Max. :5740 Max. :1.0
## monthly_average_rent average_tenant_age age_of_building total_sq_foot
## Min. :3029 Min. :41.14 Min. :39.0 Min. :19217
## 1st Qu.:3564 1st Qu.:45.14 1st Qu.:47.0 1st Qu.:41192
## Median :3813 Median :48.20 Median :49.0 Median :47096
## Mean :3687 Mean :49.92 Mean :49.4 Mean :49248
## 3rd Qu.:3864 3rd Qu.:53.88 3rd Qu.:51.0 3rd Qu.:59251
## Max. :4019 Max. :65.18 Max. :60.0 Max. :78093
## month complaints
## Min. : 1.00 Min. : 0.000
## 1st Qu.: 3.75 1st Qu.: 1.000
## Median : 6.50 Median : 2.000
## Mean : 6.50 Mean : 3.658
## 3rd Qu.: 9.25 3rd Qu.: 5.250
## Max. :12.00 Max. :18.000##number of buildings
N_buildings <- length(unique(pest_data$building_id))
N_buildings## [1] 10## arrange data into a list
stan_dat_simple <- list(
N = nrow(pest_data),
complaints = pest_data$complaints,
traps = pest_data$traps,
log_sq_foot = log(pest_data$total_sq_foot/1e4),
live_in_super = pest_data$live_in_super
)
str(stan_dat_simple)## List of 5
## $ N : int 120
## $ complaints : num [1:120] 1 3 0 1 0 0 4 3 2 2 ...
## $ traps : num [1:120] 8 8 9 10 11 11 10 10 9 9 ...
## $ log_sq_foot : num [1:120] 1.42 1.42 1.42 1.42 1.42 ...
## $ live_in_super: num [1:120] 0 0 0 0 0 0 0 0 0 0 ...comp_model_NB <- stan_model('multiple_NB_regression.stan')fitted_model_NB <- sampling(comp_model_NB, data = stan_dat_simple)##
## SAMPLING FOR MODEL 'multiple_NB_regression' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.781213 seconds (Warm-up)
## Chain 1: 0.658042 seconds (Sampling)
## Chain 1: 1.43926 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'multiple_NB_regression' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.828736 seconds (Warm-up)
## Chain 2: 0.711737 seconds (Sampling)
## Chain 2: 1.54047 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'multiple_NB_regression' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.804047 seconds (Warm-up)
## Chain 3: 0.646387 seconds (Sampling)
## Chain 3: 1.45043 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'multiple_NB_regression' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 4.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.798981 seconds (Warm-up)
## Chain 4: 0.678746 seconds (Sampling)
## Chain 4: 1.47773 seconds (Total)
## Chain 4:samps_NB <- rstan::extract(fitted_model_NB)y_rep <- samps_NB$y_rep
ppc_stat_grouped(
y = stan_dat_simple$complaints,
yrep = y_rep,
group = pest_data$building_id,
stat = 'mean',
binwidth = 0.2
)
Lab 5
Select some of the variables containing information at the building level (live in super, age of building, average tentant age, and monthly average rent) and arrange them into a matrix called building_data (nrow = number of buildings).
N_months <- length(unique(pest_data$date))
# Add some IDs for building and month
pest_data <- pest_data %>%
mutate(
building_fac = factor(building_id, levels = unique(building_id)),
building_idx = as.integer(building_fac),
ids = rep(1:N_months, N_buildings),
mo_idx = lubridate::month(date)
)
# Center and rescale the building specific data
building_data <- pest_data %>%
select(
building_idx,
live_in_super,
age_of_building,
average_tenant_age,
monthly_average_rent
) %>%
unique() %>%
arrange(building_idx) %>%
select(-building_idx) %>%
scale(scale=FALSE) %>%
as.data.frame() %>%
mutate( # scale by constants
age_of_building = age_of_building / 10,
average_tenant_age = average_tenant_age / 10,
monthly_average_rent = monthly_average_rent / 1000
) %>%
as.matrix()
str(building_data)## num [1:10, 1:4] -0.3 -0.3 -0.3 -0.3 0.7 -0.3 -0.3 0.7 -0.3 0.7 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:4] "live_in_super" "age_of_building" "average_tenant_age" "monthly_average_rent"stan_dat_hier <-
with(pest_data,
list(complaints = complaints,
traps = traps,
N = length(traps),
J = N_buildings,
log_sq_foot = log(pest_data$total_sq_foot/1e4),
building_data = building_data,
mo_idx = as.integer(as.factor(date)),
K = 4,
building_idx = building_idx
)
)
str(stan_dat_hier)## List of 9
## $ complaints : num [1:120] 1 3 0 1 0 0 4 3 2 2 ...
## $ traps : num [1:120] 8 8 9 10 11 11 10 10 9 9 ...
## $ N : int 120
## $ J : int 10
## $ log_sq_foot : num [1:120] 1.42 1.42 1.42 1.42 1.42 ...
## $ building_data: num [1:10, 1:4] -0.3 -0.3 -0.3 -0.3 0.7 -0.3 -0.3 0.7 -0.3 0.7 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:4] "live_in_super" "age_of_building" "average_tenant_age" "monthly_average_rent"
## $ mo_idx : int [1:120] 1 2 3 4 5 6 7 8 9 10 ...
## $ K : num 4
## $ building_idx : int [1:120] 1 1 1 1 1 1 1 1 1 1 ...Model the variation across buildings using a varying intercept Negative Binomial hierarchical model.
Using stan write, compile and fit the following varying intercept model
\[ \text{complaints}_{ib} \sim \text{Neg-Binomial}(\lambda_{ib}, \phi) \\ \lambda_{ib} = \exp{(\eta_{ib})} \\ \eta_{ib} = \mu_{b(i)} + \beta \, {\rm traps}_{i} + \text{log_sq_foot}_i \\ \mu_b \sim \text{Normal}(\alpha + \beta_{\rm super}\, {\rm super}_b + \ldots + \beta_{\rm mar}\, {\rm mar}_b , \sigma_{\mu}) \]
comp_model_NB_hier <- stan_model('hier_NB_regression.stan')fitted_model_NB_hier <- sampling(comp_model_NB_hier, data = stan_dat_hier,
chains = 4, #cores = 4,
iter = 4000)##
## SAMPLING FOR MODEL 'hier_NB_regression' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 5.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.52 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.13562 seconds (Warm-up)
## Chain 1: 3.45736 seconds (Sampling)
## Chain 1: 6.59298 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'hier_NB_regression' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 5.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.56 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 3.15948 seconds (Warm-up)
## Chain 2: 2.69421 seconds (Sampling)
## Chain 2: 5.85369 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'hier_NB_regression' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.88282 seconds (Warm-up)
## Chain 3: 1.91586 seconds (Sampling)
## Chain 3: 4.79867 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'hier_NB_regression' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 3.02218 seconds (Warm-up)
## Chain 4: 2.11161 seconds (Sampling)
## Chain 4: 5.13379 seconds (Total)
## Chain 4:Do you have any divergent transitions?
samps_hier_NB <- rstan::extract(fitted_model_NB_hier)print(fitted_model_NB_hier, pars = c('sigma_mu','beta','alpha','phi','mu'))## Inference for Stan model: hier_NB_regression.
## 4 chains, each with iter=4000; warmup=2000; thin=1;
## post-warmup draws per chain=2000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## sigma_mu 0.26 0.01 0.17 0.04 0.14 0.23 0.34 0.67 143 1.03
## beta -0.23 0.00 0.06 -0.35 -0.27 -0.23 -0.19 -0.12 816 1.00
## alpha 1.29 0.01 0.42 0.46 1.00 1.30 1.55 2.16 819 1.00
## phi 1.58 0.01 0.36 1.02 1.33 1.52 1.77 2.44 2248 1.00
## mu[1] 1.31 0.02 0.54 0.22 0.96 1.33 1.65 2.39 937 1.00
## mu[2] 1.26 0.02 0.52 0.23 0.93 1.25 1.59 2.33 959 1.00
## mu[3] 1.44 0.02 0.48 0.49 1.13 1.43 1.74 2.43 948 1.01
## mu[4] 1.47 0.02 0.47 0.57 1.17 1.45 1.77 2.44 829 1.01
## mu[5] 1.10 0.01 0.42 0.28 0.82 1.11 1.37 1.94 1007 1.00
## mu[6] 1.19 0.02 0.48 0.25 0.88 1.21 1.50 2.13 850 1.00
## mu[7] 1.48 0.02 0.51 0.52 1.15 1.46 1.81 2.55 891 1.01
## mu[8] 1.29 0.01 0.42 0.47 1.01 1.32 1.57 2.12 881 1.00
## mu[9] 1.44 0.02 0.56 0.30 1.07 1.47 1.81 2.52 819 1.00
## mu[10] 0.88 0.01 0.36 0.19 0.64 0.87 1.12 1.60 1078 1.00
##
## Samples were drawn using NUTS(diag_e) at Sun May 10 22:16:46 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).mcmc_trace(
as.array(fitted_model_NB_hier,pars = 'sigma_mu'),
np = nuts_params(fitted_model_NB_hier),
window = c(500,1000)
)
# assign to object so we can compare to another plot later
scatter_with_divs <- mcmc_scatter(
as.array(fitted_model_NB_hier),
pars = c("mu[4]", 'sigma_mu'),
transform = list('sigma_mu' = "log"),
np = nuts_params(fitted_model_NB_hier)
)
scatter_with_divs
N_sims <- 1000
log_sigma <- rep(NA, N_sims)
theta <- rep(NA, N_sims)
for (j in 1:N_sims) {
log_sigma[j] <- rnorm(1, mean = 0, sd = 1)
theta[j] <- rnorm(1, mean = 0, sd = exp(log_sigma[j]))
}
draws <- cbind("mu" = theta, "log(sigma_mu)" = log_sigma)
mcmc_scatter(draws)
parcoord_with_divs <- mcmc_parcoord(
as.array(fitted_model_NB_hier, pars = c("sigma_mu", "mu")),
np = nuts_params(fitted_model_NB_hier)
)
parcoord_with_divs
We see evidence that our problems concentrate when \(\texttt{sigma_mu}\) is small. Reparametrize the model using the non-centered parametrization and check diagnostics. The distribution of \(\texttt{mu}\) does not change: \[ \mu_b \sim \text{Normal}(\alpha + \texttt{building_data} \, \zeta, \sigma_{\mu}) \] We add an auxiliary variable in the parameters block \(\texttt{mu_raw}\sim \text{Normal}(0, 1)\) and we move \(\texttt{mu}\) in the transformed parameters block.
transformed parameters {
vector[J] mu;
mu = alpha + building_data * zeta + sigma_mu * mu_raw;
}comp_model_NB_hier_ncp <- stan_model('hier_NB_regression_ncp.stan')fitted_model_NB_hier_ncp <- sampling(comp_model_NB_hier_ncp, data = stan_dat_hier, #cores=4,
chains = 4)##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 8.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.85 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.76579 seconds (Warm-up)
## Chain 1: 1.05181 seconds (Sampling)
## Chain 1: 2.8176 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.35 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.77977 seconds (Warm-up)
## Chain 2: 1.03206 seconds (Sampling)
## Chain 2: 2.81183 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.90741 seconds (Warm-up)
## Chain 3: 0.992669 seconds (Sampling)
## Chain 3: 2.90008 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.72963 seconds (Warm-up)
## Chain 4: 0.980096 seconds (Sampling)
## Chain 4: 2.70973 seconds (Total)
## Chain 4:print(fitted_model_NB_hier_ncp, pars = c('sigma_mu','beta','alpha','phi','mu'))## Inference for Stan model: hier_NB_regression_ncp.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## sigma_mu 0.22 0.00 0.17 0.01 0.09 0.19 0.31 0.64 1610 1
## beta -0.23 0.00 0.06 -0.35 -0.27 -0.23 -0.19 -0.12 2678 1
## alpha 1.27 0.01 0.43 0.44 0.99 1.26 1.55 2.14 2653 1
## phi 1.57 0.00 0.35 1.01 1.34 1.53 1.77 2.37 4922 1
## mu[1] 1.30 0.01 0.54 0.26 0.94 1.30 1.67 2.40 2666 1
## mu[2] 1.24 0.01 0.52 0.24 0.89 1.23 1.58 2.28 2749 1
## mu[3] 1.41 0.01 0.49 0.48 1.09 1.40 1.72 2.40 3030 1
## mu[4] 1.45 0.01 0.48 0.51 1.12 1.44 1.76 2.42 2946 1
## mu[5] 1.09 0.01 0.41 0.28 0.83 1.09 1.35 1.92 3281 1
## mu[6] 1.20 0.01 0.48 0.25 0.88 1.20 1.52 2.14 2663 1
## mu[7] 1.47 0.01 0.52 0.46 1.12 1.46 1.82 2.53 3035 1
## mu[8] 1.26 0.01 0.43 0.43 0.97 1.25 1.53 2.13 3498 1
## mu[9] 1.44 0.01 0.56 0.31 1.08 1.43 1.80 2.55 2872 1
## mu[10] 0.88 0.01 0.37 0.19 0.63 0.87 1.12 1.65 3682 1
##
## Samples were drawn using NUTS(diag_e) at Sun May 10 22:17:58 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).scatter_no_divs <- mcmc_scatter(
as.array(fitted_model_NB_hier_ncp),
pars = c("mu[4]", 'sigma_mu'),
transform = list('sigma_mu' = "log"),
np = nuts_params(fitted_model_NB_hier_ncp)
)
bayesplot_grid(scatter_with_divs, scatter_no_divs,
grid_args = list(ncol = 2), ylim = c(-11, 1))
parcoord_no_divs <- mcmc_parcoord(
as.array(fitted_model_NB_hier_ncp, pars = c("sigma_mu", "mu")),
np = nuts_params(fitted_model_NB_hier_ncp)
)
bayesplot_grid(parcoord_with_divs, parcoord_no_divs,
ylim = c(-3, 3))
samps_NB_hier_ncp <- rstan::extract(fitted_model_NB_hier_ncp, pars = c('y_rep','inv_phi'))y_rep <- as.matrix(fitted_model_NB_hier_ncp, pars = "y_rep")
ppc_dens_overlay(stan_dat_hier$complaints, y_rep[1:200,])
ppc_stat_grouped(
y = stan_dat_hier$complaints,
yrep = y_rep,
group = pest_data$building_id,
stat = 'mean',
binwidth = 0.5
)
ppc_stat_grouped(
y = stan_dat_hier$complaints,
yrep = y_rep,
group = pest_data$building_id,
stat = 'sd',
binwidth = 0.5
)
ppc_intervals(
y = stan_dat_hier$complaints,
yrep = y_rep,
x = stan_dat_hier$traps
) +
labs(x = "Number of traps", y = "Number of complaints")
mean_y_rep <- colMeans(y_rep)
mean_inv_phi <- mean(as.matrix(fitted_model_NB_hier_ncp, pars = "inv_phi"))
std_resid <- (stan_dat_hier$complaints - mean_y_rep) / sqrt(mean_y_rep + mean_y_rep^2*mean_inv_phi)
qplot(mean_y_rep, std_resid) + hline_at(2) + hline_at(-2)
Try to expand the previous model with varying slope.
\[ \text{complaints}_{ib} \sim \text{Neg-Binomial}(\lambda_{ib}, \phi) \\ \lambda_{ib} = \exp{(\eta_{ib})}\\ \eta_{ib} = \mu_{b(i)} + \kappa_{b(i)} \, \texttt{traps}_{ib} + \text{log_sq_foot}_i \\ \mu_b \sim \text{Normal}(\alpha + \texttt{building_data} \, \zeta, \sigma_{\mu}) \\ \kappa_b \sim \text{Normal}(\beta + \texttt{building_data} \, \gamma, \sigma_{\kappa}) \]
stan_dat_hier <- readRDS('pest_data_longer_stan_dat.RDS')comp_model_NB_hier_slopes <- stan_model('hier_NB_regression_ncp_slopes_mod.stan')fitted_model_NB_hier_slopes <-
sampling(
comp_model_NB_hier_slopes,
data = stan_dat_hier,
chains = 4, #cores = 4,
control = list(adapt_delta = 0.95)
)##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp_slopes_mod' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000177 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.77 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 15.1196 seconds (Warm-up)
## Chain 1: 9.07921 seconds (Sampling)
## Chain 1: 24.1988 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp_slopes_mod' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 15.4089 seconds (Warm-up)
## Chain 2: 10.1335 seconds (Sampling)
## Chain 2: 25.5425 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp_slopes_mod' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 8.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.81 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 15.7828 seconds (Warm-up)
## Chain 3: 5.99432 seconds (Sampling)
## Chain 3: 21.7771 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'hier_NB_regression_ncp_slopes_mod' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 15.3129 seconds (Warm-up)
## Chain 4: 9.9742 seconds (Sampling)
## Chain 4: 25.2871 seconds (Total)
## Chain 4:mcmc_hist(
as.matrix(fitted_model_NB_hier_slopes, pars = "sigma_kappa"),
binwidth = 0.005
)
print(fitted_model_NB_hier_slopes, pars = c('kappa','beta','alpha','phi','sigma_mu','sigma_kappa','mu'))## Inference for Stan model: hier_NB_regression_ncp_slopes_mod.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## kappa[1] -0.02 0.00 0.08 -0.14 -0.07 -0.03 0.02 0.16 589 1.00
## kappa[2] -0.42 0.00 0.09 -0.60 -0.48 -0.41 -0.35 -0.24 1538 1.00
## kappa[3] -0.59 0.00 0.11 -0.80 -0.66 -0.59 -0.52 -0.39 3602 1.00
## kappa[4] -0.22 0.00 0.07 -0.35 -0.26 -0.22 -0.18 -0.09 2094 1.00
## kappa[5] -0.60 0.00 0.09 -0.79 -0.66 -0.60 -0.54 -0.43 4156 1.00
## kappa[6] -0.44 0.00 0.10 -0.67 -0.50 -0.43 -0.37 -0.25 2875 1.00
## kappa[7] -0.31 0.00 0.07 -0.45 -0.36 -0.31 -0.26 -0.18 4802 1.00
## kappa[8] -0.23 0.01 0.16 -0.58 -0.33 -0.22 -0.13 0.05 779 1.00
## kappa[9] 0.08 0.00 0.06 -0.03 0.04 0.08 0.12 0.19 4217 1.00
## kappa[10] -0.72 0.01 0.16 -1.02 -0.83 -0.73 -0.62 -0.37 672 1.00
## beta -0.35 0.00 0.06 -0.48 -0.39 -0.35 -0.31 -0.22 1497 1.00
## alpha 1.41 0.01 0.31 0.74 1.21 1.42 1.62 1.99 2408 1.00
## phi 1.62 0.01 0.21 1.26 1.48 1.60 1.74 2.07 439 1.01
## sigma_mu 0.50 0.03 0.43 0.02 0.17 0.38 0.71 1.61 262 1.01
## sigma_kappa 0.12 0.00 0.08 0.03 0.07 0.10 0.16 0.34 475 1.01
## mu[1] 0.30 0.03 0.73 -1.40 -0.10 0.41 0.81 1.43 593 1.00
## mu[2] 1.64 0.01 0.51 0.67 1.30 1.62 1.96 2.68 1407 1.00
## mu[3] 2.14 0.01 0.33 1.53 1.91 2.13 2.35 2.81 4150 1.00
## mu[4] 1.49 0.01 0.51 0.48 1.18 1.50 1.81 2.47 2237 1.00
## mu[5] 2.39 0.01 0.42 1.59 2.12 2.37 2.66 3.23 4258 1.00
## mu[6] 1.91 0.01 0.39 1.21 1.66 1.88 2.14 2.81 2102 1.00
## mu[7] 2.68 0.00 0.26 2.20 2.51 2.68 2.85 3.21 4351 1.00
## mu[8] -0.52 0.03 0.96 -2.33 -1.15 -0.54 0.07 1.58 917 1.00
## mu[9] 0.23 0.01 0.57 -0.86 -0.17 0.25 0.64 1.33 4200 1.00
## mu[10] 1.82 0.05 1.11 -0.75 1.22 1.99 2.60 3.59 539 1.00
##
## Samples were drawn using NUTS(diag_e) at Sun May 10 22:20:44 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).mcmc_hist(
as.matrix(fitted_model_NB_hier_slopes, pars = "beta"),
binwidth = 0.005
)
y_rep <- as.matrix(fitted_model_NB_hier_slopes, pars = "y_rep")
ppc_dens_overlay(
y = stan_dat_hier$complaints,
yrep = y_rep[1:200,]
)