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

[Mark Herzog]

   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

> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

> >

>

>       [[alternative HTML version deleted]]

>

>

_______________________________________________

R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

	[[alternative HTML version deleted]]