Nest survival: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

Thank you very much for working on this Ben and your very useful comment
Mark! I really appreciate all of your help. It looks like it was the zero in
the exposure days that was causing the problem. As the data were recently
sent on to me from a variety of sources, I hadn't spotted this! It looks
like it is due to an observer witnessing a nest mortality on the same day as
he found the active nest, which then resulted in zero exposure days.

Very grateful for all your help

Best wishes

Elwyn

On 9 March 2015 at 20:30, Mark Herzog <mherzog at usgs.gov> wrote:

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

I think it might be a result of you having an exposure period of zero days
in your data... That won't work very well given the link function....

Sorry this is based using my R package, so the codes a little different,
but here's the example:

library(nestsurvival)

mydata<-read.csv("c:/users/mherzog/Downloads/habitat-type_example.csv")

glm1<-glm(survive/trials~habitat,family=binomial(logexp(days=mydata$expos)),data=mydata)

Error: cannot find valid starting values: please specify some

mydata<-subset(mydata,expos>0)

glm1<-glm(survive/trials~habitat,family=binomial(logexp(days=mydata$expos)),data=mydata)

summary(glm1)

Call:

glm(formula = survive/trials ~ habitat, family = binomial(logexp(days =
mydata$expos)),

    data = mydata)

Deviance Residuals:

    Min       1Q   Median       3Q      Max

-2.0963  -0.8808   0.5694   0.8052   2.2708

Coefficients:

              Estimate Std. Error z value Pr(>|z|)

(Intercept)     1.9629     0.4036   4.863 1.16e-06 ***

habitatForest  -0.0342     0.4327  -0.079   0.9370

habitatHeath    0.5090     0.4126   1.234   0.2173

habitatScrub    0.8395     0.4497   1.867   0.0619 .

---

Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 787.11  on 568  degrees of freedom

Residual deviance: 621.21  on 565  degrees of freedom

AIC: 629.21

Number of Fisher Scoring iterations: 5

 -----Original Message-----

From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org
<r-sig-mixed-models-bounces at r-project.org>] On Behalf Of Ben Bolker

Sent: Sunday, March 08, 2015 6:18 PM

To: r-sig-mixed-models at r-project.org

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

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-functi

on-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

      [[alternative HTML version deleted]]