Binomial regression

[sidravi1]

Hi,

Thank you for putting this package together.
I'm trying to fit a binomial model but unsure how to specify the number of failures (or total trials).

data:

> d_sample
  y2                        id      lat      lon   N
1 25 IND0070940639070004_Del_U 28.83111 77.16795 100
2 37 IND0070940639100003_Del_U 28.80293 77.12073 120
3 25 IND0070940639110007_Del_U 28.78427 77.20626 130
4 39 IND0070940639130005_Del_U 28.77438 77.16210 100
5 31 IND0070940639140018_Del_U 28.76010 77.14313 110
6 23 IND0070940639200027_Del_U 28.76169 77.01645  90

models:

I'm not specifying N here.

m1 <- glmmfields(y2 ~ as.factor(id), 
           data = d_sample,  family = binomial(), lat = "lat", lon = "lon", 
           nknots = 2, iter = 500,  
           prior_sigma = half_t(3, 0, 3),
           prior_gp_theta = half_t(3, 0, 10),
           prior_gp_sigma = half_t(3, 0, 3),)

gives

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: bernoulli_logit_lpmf: n[1] is 25, but must be in the interval [0, 1]  (in 'model_glmmfields' at line 224)

Surprised that is using the bernoulli_logit_lpmf instead of binomial_logit_lpmf.

Here's another model

> m2 <- glmmfields(cbind(y2, N - y2) ~ as.factor(id), 
              data = d_sample, family = binomial(), lat = "lat", lon = "lon", 
              nknots = 2, iter = 500,  
              prior_sigma = half_t(3, 0, 3),
              prior_gp_theta = half_t(3, 0, 10),
              prior_gp_sigma = half_t(3, 0, 3),)

gives

Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) : 
  Exception: mismatch in dimension declared and found in context; processing stage=data initialization; variable name=stationID; position=0; dims declared=(12); dims found=(6)  (in 'model_glmmfields' at line 6)

Session info

> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] yaml_2.2.1       pool_0.1.4.3     DBI_1.1.0        lubridate_1.7.9  readxl_1.3.1     lme4_1.1-23      Matrix_1.2-18    here_0.1        
 [9] styler_1.3.2     cowplot_1.0.0    broom_0.7.0      geojsonio_0.9.2  dplyr_1.0.0      ggplot2_3.3.2    glmmfields_0.1.4 Rcpp_1.0.5  

[seananderson]

Yes, it looks like we used bernoulli_logit_lpmf instead of binomial_logit_lpmf and never set up a way of passing in the number of trials. That should be a relatively easy fix, although I don't have time at this very moment. Feel free to make a pull request and if not I will try to get to this shortly.

[sidravi1]

Thanks! I'll let you know if I find time to put together a PR. It probably won't be till mid-Sept.

[ericward-noaa]

Thanks for the interest in the binomial N. I ended up tweaking the code to include this now, switching from bernoulli_lpmf to binomial_lpmf. I added a new argument letting you specify the name of the N parameter in the data frame, as a character string. When it's NULL, it defaults to the bernoulli case (N = 1). m <- glmmfields(y ~ 0, time = "time",binomial_N=NULL, lat = "lat", lon = "lon", data = s$dat, family=binomial(link="logit"), nknots = 12, iter = 100, chains = 2, seed = 1) But if N is stored in your data frame as something called "bin_N", you can include it with m <- glmmfields(y ~ 0, time = "time",binomial_N="bin_N", lat = "lat", lon = "lon", data = s$dat, family=binomial(link="logit"), nknots = 12, iter = 100, chains = 2, seed = 1)

…

On Wed, Aug 26, 2020 at 7:04 AM Sid Ravinutala ***@***.***> wrote: Thanks! I'll let you know if I find time to put together a PR. It probably won't be till mid-Sept. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#9 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABGQEZAMXINOOF74X2JKSTLSCUI5HANCNFSM4QAUIN5Q> .

-- Eric Ward Northwest Fisheries Science Center NOAA Fisheries | U.S. Department of Commerce Office: 206-302-1745 Mobile: 206-650-7401 https://faculty.washington.edu/warde/

[ericward-noaa]

Leaving this issue open for the time being, because we might want to set this up slightly differently. Two more common options are to (1) rename the 'binomial_N' argument to glmmfields to be 'weights' -- and potentially allow for weights to be used elsewhere in the GLMMs, and (2) pass in arguments in terms of a matrix where each row is cbind(successes,failures)~ ...
Looking into glmmTMB, I think they allow both approaches.

[m-moinuddin]

Hi, I have tried with binomial_N but it giving me an error "Error: passing unknown arguments: binomial_N.".

[ericward-noaa]

Sorry you're having trouble. Here's a simple example that hopefully shows how to set up the data. We should add something like this to the vignette to make it clear,

library(glmmfields)
set.seed(1)
# simulate data -- spatial, because just 1 draw
s <- sim_glmmfields(n_draws = 1, n_knots = 12, gp_theta = 1.5,
gp_sigma = 0.2, sd_obs = 0.2)

simulate some binomial variables, and let N vary a bit by observation

s$dat$p = plogis(s$dat$y)
s$dat$N = sample(10:14, size=nrow(s$dat), replace=TRUE)
s$dat$n = rbinom(n = nrow(s$dat), size=s$dat$N, prob = s$dat$p)

fit the model

m <- glmmfields(n ~ 0, time = "time",binomial_N = "N",
  family = binomial(),
 lat = "lat", lon = "lon", data = s$dat,
 nknots = 12, iter = 1000, chains = 2, seed = 1)