[R-sig-ME] [FORGED] Working formula for glmm model in R (bernoulli response)
Rolf Turner
r@turner @end|ng |rom @uck|@nd@@c@nz
Fri Mar 20 00:33:44 CET 2020
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
More information about the R-sig-mixed-models
mailing list