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

Ben Bolker bbolker at gmail.com
Mon Mar 9 02:17:57 CET 2015

```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.
>
> 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")
}
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="")
mu.eta = mu.eta, valideta = valideta,
}

##Read in data, called 'mydata'

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:
> > >
> > >
> > > 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]]
>
>

```