rstanarm
rstanarm copied to clipboard
posterior_predict fails for binomial models response as cbind(success, failure)
Summary:
posterior_predict errors when supplying a stan_gamm4 binomial model with response specified as cbind(success, failure).
The error:
Error in cbind(incidence, size - incidence) :
object 'incidence' not found
Reproducible Steps:
dat <- structure(list(Group = c("b_t", "c_b", "crp", "eqm", "frs", "gro",
"mil", "min", "mmi", "mvo", "nmm", "oap", "ocr", "ofd", "omf",
"omt", "osd", "p_c", "pcr", "pfb", "sgr", "twl", "vol", "wht",
"wol", "wpp"), incidence = c(2, 0, 161, 427, 0, 0, 0, 1, 171,
204, 94, 0, 6, 4, 50, 0, 0, 8, 0, 0, 0, 18, 0, 0, 0, 35), size = c(52814,
17661, 224917, 328780, 1798, 1584, 1686, 7102, 190150, 267264,
60457, 208, 22534, 270418, 236736, 2177, 5830, 1236, 12196, 308,
3665, 105708, 5946, 112, 284, 79895), var1 = c(5245.67753079493,
0.953564060294559, 82343.9289806345, 197659.209119104, 63.2736859428056,
22.4404105569735, 1351.14684575416, 2769.92896596882, 44712.1423586455,
105171.025085931, 7071.26242155988, 119.660910329725, 430.228890875555,
14324.9500484617, 31603.4704954223, 1934.99222081046, 82.3389473082238,
42959.0387724196, 134.618019333429, 9.01855742076216, 349.877180671831,
39542.0596111563, 2196.30701684163, 6.25454990709854, 112.068551483099,
14049.389438344)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA,
-26L))
mod <-rstanarm::stan_gamm4(cbind(incidence, size - incidence) ~ s(var1),
random = ~
(1|Group),
data = dat,
family = binomial)
rstanarm::posterior_predict(mod)
Error message
Error in cbind(incidence, size - incidence) :
object 'incidence' not found
I think the issue is that the variables in cbind aren't correct in the model data object
names(mod$data)
[1] "cbind(incidence, size - incidence)" "var1" "Group"
[4] "temp_y"
It is important to note that this error does not occur when random effects are not included in the model.
mod2 <-rstanarm::stan_gamm4(cbind(incidence, size - incidence) ~ s(var1),
data = dat,
family = binomial)
rstanarm::posterior_predict(mod2)
Also note that the model data object in the non-random effect model is correctly specified
names(mod2$data)
[1] "Group" "incidence" "size" "var1" "temp_y"
RStanARM Version:
packageVersion("rstanarm")
[1] ‘2.21.3’
R Version:
getRversion()
[1] ‘4.2.0’
Operating System:
MacOS 12.5.1