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

Thank you very much, Rolf!
Your advise helped much!
lme4 works well, although is produces some warninngs:
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,? :
 ? Model failed to converge with max|grad| = 0.527726 (tol = 0.001, 
component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,? :
 ? Model is nearly unidentifiable: very large eigenvalue
 ?- Rescale variables?;Model is nearly unidentifiable: large eigenvalue 
ratio
 ?- Rescale variables?
Nevertheless, it says that there?s a significant time trend in occupancy:
Fixed effects:
 ????????????? Estimate Std. Error z value Pr(>|z|)
(Intercept) 94.7762510? 1.0176729?? 93.13?? <2e-16 ***
year??????? -0.0466590? 0.0005135? -90.86?? <2e-16 ***
======================================
glmm also works well, though quite slow. What I thought to be an 
infinite loop, turned to be very long loop. It is running about an hour, 
then produces some output. It also warns about lack of convergence.
In glmmTMB package the model doesn?t converge at all, so it is useless 
for me since I?m interested more in p-values of time trend than in 
slope. I do not see how it is possible to rescale since they are binary.

I?m still working on it. Probably there will be more questions if you 
don?t mind.

Best regards,
Michael

    ???????, 20 ????? 2020, 10:33 +11:00 ?? Rolf Turner
    <r.turner at auckland.ac.nz>:

    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