[R-sig-ME] Logistic exposure glmer polynomial covariate issue
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Mon May 17 03:13:11 CEST 2021
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
>
More information about the R-sig-mixed-models
mailing list