[R-sig-ME] [FORGED] Working formula for glmm model in R (bernoulli response)

[Rolf Turner]

On 20/03/20 2:14 am, Michael Romanov via R-sig-mixed-models wrote:

> Dear all,
>   
> Sorry for the silly question, but I got stuck on it.

The only silly question is the one that is not asked.
>   
> My dataset is occupancy of eagle territories in different years (see
> the example below). I’m trying to test my data for time trend, taking
> into account territory quality (as a random effect).
>   
> territory year occupied
> CHV-1 2006 0
> CHV-12 2006 0
> CHV-120 2006 1
> CHV-13 2009 0
> CHV-14 2009 1
> CHV-15 2010 1
> CHV-16 2010 1
>   
> My thoughts are following. ‘year’ is a fixed effect variable and
> ‘territory’ is a random effect variable. For example, in lme4 package
> the formula could be following: occupied = year + (1|territory).
> However, lme4 doesn’t support Bernoulli family,

This is simply NOT TRUE!!!  The Bernoulli family is just a special case 
of the binomial family. See below.

OTOH, the glmer() function from lme4 seems reluctant to handle the 
simple model that you appear to have in mind.  (Again, see below.)

This, in my experience, is a recurring problem with lme4 --- it cannot 
deal with the simple models that constitute edge cases.

(It is also quite possible that I've got the syntax wrong.  Perhaps a 
younger and wiser head will chime in.)

> so I chose glmm package. According to its specifications, the glmm
> function accepts two different formulae, separately for fixed and random
> effects: >
> glmm(fixed, random, varcomps.names, data, family.glmm, m,
> varcomps.equal, doPQL = TRUE,debug=FALSE, p1=1/3,p2=1/3, p3=1/3,
> rmax=1000,iterlim=1000, par.init, zeta=5, cluster=NULL)
> 
> So now I need two different formulae. I tried ‘occupied ~ year’,
> ‘occupied ~ territory’, but it doesn’t work. Namely, the RStudio gets
> into an infinite loop and doesn’t produce any output.

Please DO NOT confuse RStudio with R.  RStudio is simply an interface 
that allows those mentally handicapped people who can only use GUIs to 
access R.  :-)

>   
> Can anyone help me to write a working formula (or formulae) to properly
> run the glmm model?

I have no experience with glmm and don't really understand what you are 
trying to model.  Moreover you did not supply a reproducible example, as 
you are instructed to do by the Posting Guide.  This makes it very 
difficult to help you.  However I persisted, and constructed a simulated 
toy example data set, and fitted a model to it:

library(glmm)
set.seed(42)
X <- data.frame(territory=factor(paste0("CHV-",
                                  sample(6:10,200,TRUE))),
                 year=sample(2001:2020,200,TRUE),
                 occupied=sample(0:1,200,TRUE))
fit1 <- glmm(fixed=occupied~year,random=occupied~territory,
             family=bernoulli.glmm,varcomps.names="territory",
             data=X,m=1e3)

No infinite loops seem to have arisen.  Note that coef(fit1) gives:

> (Intercept)        year 
> 23.17057405 -0.01164075 

The "variance" component estimate is:
>            Estimate ... 
> territory 1.533e-05 ...

Now back to lme4:

X$unoccupied <- 1-X$occupied
library(lme4)
fit2 <- glmer(cbind(occupied,unoccupied) ~ year + (1|territory),
               family=binomial(),data=X)

This seems not to work properly; it says:

> boundary (singular) fit: see ?isSingular

However the fixed effects in fit2 are effectively the same as those in fit1:
> Fixed Effects:
> (Intercept)         year  
>    23.18216     -0.01164  

The variance component is said to be 0.

I then tried the glmmTMB() package:

library(glmmTMB)
fit3 <- glmmTMB(cbind(occupied,unoccupied) ~ year + (1|territory),
               family=binomial(),data=X)

This gave a mild squawk about a "Model convergence problem" but produced 
an answer with essentially the same fixed effect estimates as the other 
two methods:

>             Estimate ...
> (Intercept) 23.18216 ...
> year        -0.01164

The variance component estimate was vanishingly small: 1.062e-91.

Bottom line:  I think I'd go with the glmmTMB package for your problem,
but glmm seems to work, if you want to stick with that.

HTH

cheers,

Rolf Turner

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276