[R-sig-ME] Logistic exposure glmer polynomial covariate issue

[Ben Bolker]

On 5/15/21 9:21 AM, Kirsti Carr wrote:
> I am having issues with a glmer logistic regression with a logistic
> exposure link function (method popularized in 2004 paper by Terry Shaffer),
> to analyze juvenile survival data in respect to land cover data around
> nest sites.

   I'll start by saying that while I appreciate why people use the 
logistic exposure link function, the approach that uses the 
complementary log-log link (described here: 
https://www.rpubs.com/bbolker/logregexp ) seems very much more sensible 
to me and may work more smoothly.
> 
> Had no issues with rescaling warnings when I ran data as glm, but when I
> added the random effect I started getting these warnings for my polynomial
> models. I attempted to rescale by log transforming the covariates. This
> helped with some, but not all. I am also still getting an error message for
> some models - both models included get the warning or error. I am using
> lme4 version 1.1.25.
> 
> Any help would be much appreciated! Thanks in advance.
> 
> Example data here
> <https://www.dropbox.com/s/0w3csm2vfcj47yq/data_email.csv?dl=0>

   I have looked at this a bit and may post more tomorrow.  Things that 
jump out at me:

- `with(FL1,table(fate))` shows 237 surviving nests, only 6 being 
predated.  This means your effective sample size is 6 (!), which means 
it will only work to fit very simple models (the rule of thumb from 
Frank Harrell's book and elsewhere is that you should aim for 10 data 
points *per parameter you're going to try to estimate*). Thus a cubic 
model is probably overkill.

- the whole log(1+G250) thing looks odd to me, unless there is some 
particular reason that G250 makes more sense on a log scale (log(1+x) is 
always a bit of a caution that you're trying to treat something as 
logarithmic that might not want to be ...). The values of G250 are 
pretty reasonably spaced as it is, there's no obvious *computational* 
reason to want to log-transform them.

  How far I got so far:

gm1 <- glmer(fate~(1|brood) + G250 + offset(log(Interval)),
              family=binomial(link="cloglog"),
              data=FL1)

gm2 <- update(gm1, . ~ . -G250 + poly(G250,2))
gm3 <- update(gm1, . ~ . -G250 + poly(G250,3)) ## warning
library(bbmle)
AICctab(gm1, gm2, gm3, mnames=c("linear","quad","cubic"))

        dAICc df
quad   0.0   4
cubic  2.1   5
linear 2.8   3

The cubic model gives a warning about large eigenvalues, and the AICc 
calculation says that it's not worth it.

If you want to go nuts, here are the quadratic and linear models:

library(sjPlot)
library(ggeffects)

g1 <- ggpredict(gm1, terms="G250 [0:8.5 by=0.2]")
plot_model(gm2, type="pred", terms="G250 [0:8.5 by=0.2]") +
   geom_line(data=g1, colour="red") +
   geom_ribbon(data=g1, colour=NA, fill="red", alpha=0.2,
               aes(ymin=conf.low, ymax=conf.high)) +
   stat_sum(data=FL1, alpha=0.5, aes(x=G250, y=fate))


   This is all with cloglog - it could probably be done with logexp() 
too, but would all be a little less stable / more fussy.

> 
> My code (logexp function used from posts by Ben Bolker online):
> 
> library(MASS)
> library(lme4)
> packages(bbmle)
> 
> logexp <- function(expos = 1)
> {
>    linkfun <- function(mu) qlogis(mu^(1/expos))
>    ## FIXME: is there some trick we can play here to allow
>    ##   evaluation in the context of the 'data' argument?
>    linkinv <- function(eta)  plogis(eta)^expos
>    mu.eta <- function(eta) expos * plogis(eta)^(expos-1) *
>      .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
>    valideta <- function(eta) TRUE
>    link <- paste("logexp(", deparse(substitute(expos)), ")",
>                  sep="")
>    structure(list(linkfun = linkfun, linkinv = linkinv,
>                   mu.eta = mu.eta, valideta = valideta,
>                   name = link),
>              class = "link-glm")
> }
> 
> FL1 <- read.csv("data_email.csv")
> 
> g250.3
> <-glmer(fate~(1|brood)+log(I(1+G250))+(log(I(1+G250^2)))+(log(I(1+G250^3))),
> family=binomial(link=logexp(FL1$Interval)), data=FL1, start=NULL)
> 
> #error I received:
> 
> # Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, maxit
> = 100L,  :
> #(maxstephalfit) PIRLS step-halvings failed to reduce deviance in
> pwrssUpdate
> 
> g500.2 <-glmer(fate~(1|brood)+log(I(1+G500))+(log(I(1+G500^2))),
> family=binomial(link=logexp(FL1$Interval)), data=FL1, start=NULL)
> 
> # Warning for g500.2:
> 
> # Warning message:
>    # In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>                  # Model is nearly unidentifiable: large eigenvalue ratio
>                 # Rescale variables?
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>