[R-sig-ME] Nest survival: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

[Ben Bolker]

Elwyn Sharps <e.sharps at ...> writes:

> 
> Hi Ben
> 
> Thank you very much for your reply. If you click on this link, it should
> give you the data in a CSV file.
> https://sites.google.com/site/es263datahlp/habitat-type_example.csv
> 
> Many thanks
> 
> Elwyn

  I worked on this for a while, without complete success. The main
issue is that the inverse-link function and derivative functions
need some clamping so that they don't hit 0/1 ... this still doesn't
solve the lme4 problem, but at least it allows the GLM to work.

  Have you considered a cloglog link + offset(log(exposure))
model?  That *might* be a little more stable ...


library(lme4)
library(MASS)

logexp <- function(exposure = 1, eps=1e-8, maxlink=Inf)
{
  linkfun <- function(mu) {
      r <- qlogis(mu^(1/exposure))
      ## clamp link function: not actually necessary?
      ##  maxlink set to Inf
      if (any(toobig <- abs(r)>maxlink)) {
          ## cat("max threshold hit") 
          r[toobig] <- sign(r[toobig])*maxlink
      }
      return(r)
  }
  ## utility for clamping inverse-link, derivative function
  clamp <- function(x) {
      x <- pmax(eps,x)
      if (upr) x <- pmin(1-eps,x)
      return(x)
  }
  linkinv <- function(eta)  clamp(plogis(eta)^exposure)
  mu.eta <- function(eta) {
      r <- exposure * clamp(plogis(eta)^(exposure-1)) *
          .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
      return(r)
  }
  valideta <- function(eta) TRUE
  link <- paste("logexp(", deparse(substitute(exposure)), ")",
                sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
                 mu.eta = mu.eta, valideta = valideta,
                 name = link),
            class = "link-glm")
}

##Read in data, called 'mydata'

mydata <- read.csv("habitat-type_example.csv")

library("ggplot2")
with(mydata,table(survive,trials))
with(mydata,table(survive,habitat))
ggplot(mydata,aes(log(1+expos),survive,colour=habitat))+
    geom_point()+
    geom_smooth(method="glm",family="binomial")
ggplot(subset(mydata,habitat=="Conregrowth"),
       aes(expos,survive))+
           stat_sum(aes(size=..n..))+
    geom_smooth(method="glm",family="binomial")+
        scale_size_area()

## trials is always == 1 in this data set

## the fact that glm() fails means that the problem is more
## basic than a GLMM problem
glm1 <- glm(survive~habitat,
            family=binomial(logexp(exposure=mydata$expos)),
            data=mydata)

Mod1 <- glmer(survive~habitat + (1|site)+(1|year),
              family=binomial(logexp(exposure=mydata$expos)),data=mydata,
              nAGQ=1,
              devFunOnly=TRUE,
              control=glmerControl(nAGQ0initStep=FALSE),
              start=list(beta=coef(glm1),theta=1e-5),
              verbose=100)

Mod2 <- glmer(survive~habitat + (1|year),
              family=binomial(logexp(exposure=mydata$expos)),data=mydata,
              start=list(theta=c(1e-6,1e-6)),
              nAGQ=0,
              devFunOnly=TRUE)

Mod3 <- glmer(survive~habitat + (1|site),
              family=binomial(logexp(exposure=mydata$expos)),data=mydata,
              start=list(theta=c(1e-6,1e-6)),
              nAGQ=0,
              devFunOnly=TRUE)

mydata3 <- droplevels(subset(mydata,habitat!="Conregrowth"))
Mod4 <- glmer(survive~habitat + (1|year),
              family=binomial(logexp(exposure=mydata3$expos)),data=mydata3)

Mod5 <- glmer(survive~habitat + (1|site),
              family=binomial(logexp(exposure=mydata3$expos)),data=mydata3,
              nAGQ=1,
              devFunOnly=TRUE,
              control=glmerControl(nAGQ0initStep=FALSE),
              start=list(beta=coef(glm1),theta=1e-5),
              verbose=100)

with(mydata3,table(site,habitat,survive))
with(mydata,table(year,habitat,survive))

> 
> On 5 March 2015 at 03:11, Ben Bolker <bbolker at ...> wrote:
> 
> > Elwyn Sharps <e.sharps <at> ...> writes:
> >
> > >
> >
> > [snip]
> >
> > > I am using a nest survival model (glmer) with random effects and a
> > logistic
> > > exposure link function, as described here:
> > >
> > >
> > http://stackoverflow.com/questions/19012128/user-defined-link-function-for-
> > > glmer-for-known-fate-survival-modelling
> > >
> > > I am running a number of different models, with varying fixed effects.
> > Some
> > > of them are running well, with no error or warning messages, however
> > > for other models, I am getting the following message:
> > >
> > > *Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in
> > > pwrssUpdate*
> > >
> > > I'm not sure what is causing this error. I have tried to check the data
> > for
> > > simple problems, however can't see anything that could be causing
> > trouble.
> > >
> > > I've also tried running the model without the random effects. This
> > results
> > > in a different error message:
> > >
> > > *Error: cannot find valid starting values: please specify some*
> >
> >   Example data doesn't seem to be attached: it may have been
> > stripped by the mailing list software.  Can you post it somewhere
> > public and provide a URL?
> >
> >   My guess it that there is something rather wonky about the data
> > for this example, e.g. complete separation (for example, no individuals
> > die for some combination of predictor variables).  Hard to say
> > without the data though.
> >
> >   Ben Bolker
> >
> > _______________________________________________
> > R-sig-mixed-models at ... mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> 
> 	[[alternative HTML version deleted]]
> 
>