1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
params <-
list(EVAL = TRUE)
## ----SETTINGS-knitr, include=FALSE--------------------------------------------
stopifnot(require("knitr"))
library("bayesplot")
knitr::opts_chunk$set(
dev = "png",
dpi = 150,
fig.asp = 0.618,
fig.width = 5,
out.width = "60%",
fig.align = "center",
comment = NA,
eval = if (isTRUE(exists("params"))) params$EVAL else FALSE
)
## ----pkgs, include=FALSE------------------------------------------------------
library("ggplot2")
library("rstanarm")
set.seed(840)
## ----eval=FALSE---------------------------------------------------------------
# library("bayesplot")
# library("ggplot2")
# library("rstanarm")
## ----roaches-data-------------------------------------------------------------
head(roaches) # see help("rstanarm-datasets")
roaches$roach100 <- roaches$roach1 / 100 # pre-treatment number of roaches (in 100s)
## ----roaches-model-pois, message=FALSE----------------------------------------
# using rstanarm's default priors. For details see the section on default
# weakly informative priors at https://mc-stan.org/rstanarm/articles/priors.html
fit_poisson <- stan_glm(
y ~ roach100 + treatment + senior,
offset = log(exposure2),
family = poisson(link = "log"),
data = roaches,
seed = 1111,
refresh = 0 # suppresses all output as of v2.18.1 of rstan
)
## ----print-pois---------------------------------------------------------------
print(fit_poisson)
## ----roaches-model-nb, message=FALSE------------------------------------------
fit_nb <- update(fit_poisson, family = "neg_binomial_2")
## ----print-nb-----------------------------------------------------------------
print(fit_nb)
## ----y------------------------------------------------------------------------
y <- roaches$y
## ----yrep---------------------------------------------------------------------
yrep_poisson <- posterior_predict(fit_poisson, draws = 500)
yrep_nb <- posterior_predict(fit_nb, draws = 500)
dim(yrep_poisson)
dim(yrep_nb)
## ----ppc_dens_overlay---------------------------------------------------------
color_scheme_set("brightblue")
ppc_dens_overlay(y, yrep_poisson[1:50, ])
## ----ppc_dens_overlay-2, message=FALSE, warning=FALSE-------------------------
ppc_dens_overlay(y, yrep_poisson[1:50, ]) + xlim(0, 150)
## ----ppc_hist, message=FALSE--------------------------------------------------
ppc_hist(y, yrep_poisson[1:5, ])
## ----ppc_hist-nb, message=FALSE-----------------------------------------------
ppc_hist(y, yrep_nb[1:5, ])
## ----ppc_hist-nb-2, message=FALSE---------------------------------------------
ppc_hist(y, yrep_nb[1:5, ], binwidth = 20) +
coord_cartesian(xlim = c(-1, 300))
## ----prop_zero----------------------------------------------------------------
prop_zero <- function(x) mean(x == 0)
prop_zero(y) # check proportion of zeros in y
## ----ppc_stat, message=FALSE--------------------------------------------------
ppc_stat(y, yrep_poisson, stat = "prop_zero", binwidth = 0.005)
## ----ppc_stat-nb, message=FALSE-----------------------------------------------
ppc_stat(y, yrep_nb, stat = "prop_zero")
## ----ppc_stat-max, message=FALSE----------------------------------------------
ppc_stat(y, yrep_poisson, stat = "max")
ppc_stat(y, yrep_nb, stat = "max")
ppc_stat(y, yrep_nb, stat = "max", binwidth = 100) +
coord_cartesian(xlim = c(-1, 5000))
## ----available_ppc------------------------------------------------------------
available_ppc()
## ----available_ppc-grouped----------------------------------------------------
available_ppc(pattern = "_grouped")
## ----ppc_stat_grouped, message=FALSE------------------------------------------
ppc_stat_grouped(y, yrep_nb, group = roaches$treatment, stat = "prop_zero")
## ----pp_check.foo-------------------------------------------------------------
# @param object An object of class "foo".
# @param type The type of plot.
# @param ... Optional arguments passed on to the bayesplot plotting function.
pp_check.foo <- function(object, type = c("multiple", "overlaid"), ...) {
type <- match.arg(type)
y <- object[["y"]]
yrep <- object[["yrep"]]
stopifnot(nrow(yrep) >= 50)
samp <- sample(nrow(yrep), size = ifelse(type == "overlaid", 50, 5))
yrep <- yrep[samp, ]
if (type == "overlaid") {
ppc_dens_overlay(y, yrep, ...)
} else {
ppc_hist(y, yrep, ...)
}
}
## ----foo-object---------------------------------------------------------------
x <- list(y = rnorm(200), yrep = matrix(rnorm(1e5), nrow = 500, ncol = 200))
class(x) <- "foo"
## ----pp_check-1, message=FALSE------------------------------------------------
color_scheme_set("purple")
pp_check(x, type = "multiple", binwidth = 0.3)
## ----pp_check-2---------------------------------------------------------------
color_scheme_set("darkgray")
pp_check(x, type = "overlaid")
|
