[R-sig-ME] Nest survival: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
Elwyn Sharps
e.sharps at gmail.com
Wed Mar 11 11:43:19 CET 2015
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
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list