rstanarm icon indicating copy to clipboard operation
rstanarm copied to clipboard

posterior_predict fails for binomial models response as cbind(success, failure)

Open jscamac opened this issue 2 years ago • 0 comments

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

jscamac avatar Mar 15 '23 01:03 jscamac