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

Ben Bolker bbolker at gmail.com
Wed Mar 11 13:05:41 CET 2015 

  Huh.  I saw that zero-exposure case, but clamping the inverse-link
function seemed to make the GLM work, so I didn't think that it would
still be screwing up the GLMM ...

On Wed, Mar 11, 2015 at 6:43 AM, Elwyn Sharps <e.sharps at gmail.com> wrote:
>
> 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
>>
>> > > 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]]
>>
>> _______________________________________________
>> R-sig-mixed-models at 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