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