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