[R-sig-ME] glmmADMB function maximizer failure
Ben Bolker
bbolker at gmail.com
Sat Feb 16 18:19:35 CET 2013
[cc'ing back to r-sig-mixed]
On 13-02-16 08:33 AM, Adriaan De Jong wrote:
> Hi,
>
> Some time ago I posted a question to the r-sig-mixed-models forum (copy
> below), but, much to my surprise, received no answer yet (although the
> R-sig-mixed-models Archives
> <https://stat.ethz.ch/pipermail/r-sig-mixed-models/> stored an answer to
> another question under this question’s header).
I saw the response from Billy Requena at <
http://article.gmane.org/gmane.comp.lang.r.lme4.devel/9591 >
recommending that you try admbControl(shess=FALSE,noinit=FALSE) ; I had
hoped that might work for you in some cases.
> Since I posted my question, I found that others have encountered this
> error and that the issue has been discussed by you and others on the
> r-sig-mixed-models forum. From what I’ve seen in the archives and from
> the results of further web-searches, I’ve come to the conclusion that
> this matter is not finally settled. Am I wrong? If not, I hope this
> e-mail can contribute to the solution.
This matter is certainly not finally settled.
Apologies in advance if I sound a bit defensive about this, but:
nonlinear modeling is *hard*, and especially hard on small, noisy data
sets that are so common in ecology. (I'm tempted to go back
and italicize/boldface/uppercase that a few times.)
Just because more sophisticated models exist and we'd like to use
them (of course I would like to allow for all of the realistic
complications that I know are out there), doesn't mean it's
practical to fit them to data, no matter how slick the computational
tools are.
I'm not trying to blame the data or make light of the difficulty
of gathering it -- just saying that it's really hard to
fit a complicated model to a sparse data set.
glmmADMB attempts to be a black-box, one-size-fits-all solution.
It is often possible to do a bit better by hand-tuning the model in
AD Model Builder, but this requires quite a bit of expertise and
knowledge about how the model-fitting procedure works.
Furthermore, glmmADMB errs a tiny bit on the side of caution by
refusing to give an answer when the estimated Hessian matrix
(the matrix of second derivatives of the log-likelihood with respect
to the parameters, which when inverted gives an estimate of the
variances/standard) is not well-behaved (positive definite). *Sometimes*
the answers might be reasonable even when the Hessian is bad, but
it's not possible to know without checking carefully, and it would
be very easy to get nonsensical answers this way. I've considered
adding an option to allow fits to be retrieved even when the Hessian
is bad, but I'm afraid it would be badly misused ...
It's sometimes the case that exactly the same glmmADMB fit behaves
differently on different platforms (Windows, MacOS, Linux, 32 vs
64-bit). This is frustrating, but it generally means that the model
is very close to the edge of stability in any case ...
More generally, I'm not entirely enamored of the AIC-inspired
"run all possible cases" scenario. I can appreciate wanting to
make a reasonable job estimating an appropriate place in the bias-variance
tradeoff curve, or get good model-averaged predictions, but I
wonder what your goal is here. Are you trying to test hypotheses
about the effect of the treatment? Are you trying to figure out
the best predictive model for management purposes? Admittedly it
would be nice if there were good, flexible penalized regression models
(IMO a better way to do good estimation when the appropriate variables
are unknown than AIC-based model averaging) that could handle the
things ecologists like to do (negative binomial models,
random effects, etc.)).
The bottom line is that I'm sorry that glmmADMB doesn't work
better than it does, but it's very hard to give any general prescriptions
other than (1) accept that some data are too sparse to fit some models;
(2) work up your own models from scratch using WinBUGS/JAGS or
AD Model Builder, so you have more control to tweak settings
Depending on what I was doing with this problem, I would try:
* check whether the among-site variances seem large in general;
if not (I bet they're not), simplify life by dropping the random
effects and just fitting Poisson or NB GLMs
* use permutation tests to allow for hypothesis testing with
fewer distributional assumptions
* rather than fitting each species separately, try modeling all
the data together, with a random effect of species (and
species-by-treatment).
[This may sound contradictory since I'm encouraging you to simplify
your models, but this would have the advantage of allowing the
rarer species to "share strength" from the commoner species ...]
> I also continued my tests of glmmADMB. The data I use are of 13 species
> from the very same investigation (nine years of territory mapping of
> breeding birds in 19 different sites), but because I’ve excluded sites
> without any observation during the nine year’s period from the analysis
> of that species
This seems reasonable, as those territories can't contribute
any information about the dynamics of that bird through time
or with respect to treatment,s although I would be very careful with it.
> the numbers of cases varies among species, and of
> course, species range from rare to common (all data files attached).
>
> The pass/fail pattern of the six models over the species was as follows:
>
> Spec. m100 m200 m300 m400 m500 m600 (models
> described in the copy of my forum posting below)
> Aa x + + +
> x x
>
> Ap + x + x
> x x
>
> Cd + x + x
> + +
>
> Ce + x + +
> + x
>
> Gg x x + x
> x x
>
> Hr + x + x
> + +
>
> Lc + x + x
> + +
>
> Mf x x x +
> x x
>
> Na + x + x
> + x
>
> Sr x x + x
> + x
>
> Sv x x + +
> x x
>
> To + x + x
> x +
>
> Vv x + + x
> + +
>
> Clearly, the negative binomial model without zeroinflation correction
> (m300) worked fine for most species. Unfortunately, all the failing
> models made the use of AICctab and the graphs presented in the “Getting
> started with the glmmADMB package” rather meaningless.
>
> After scrutinizing the data further, I concluded that the Cd, Ce, Gg, Lc
> and To data had few >1 observations, and thus were better treated with
> binomial glmer models.
>
> The distribution of the data also “explained” why the Mf data could be
> fitted with zeroInflation = TRUE but not with zeroInflation = FALSE.
> This was the only (non-binomial) dataset with obvious zero inflation.
>
> From there I ran “nbinom” models for the eight non-binomial species
> (with zeroInflation = TRUE for the Mf data) with the following
> combinations of fixed and random factors:
>
> treatment * year + (year|site)
> treatment + year + (year|site)
> year + (year|site)
> treatment + (year|site)
> 1 + (year|site)
> treatment * year + (1|site)
> treatment + year + (1|site)
> year + (1|site)
> treatment + (1|site)
> 1 + (1|site)
>
> This worked fine in all species except Ap. Here, the 1 + year|site model
> produced the “The function maximizer failed (couldn't find STD file)”
> error, while the 1|site, the model ran fine. This seems to contradict
> the suggestion by Joshua Wiley (“What happens if you try a simpler
> model”) when answering Melanie Browne (Jun 24, 2012). Here the next
> simplest model produced the error while the more complex models didn’t.
>
> Basically, I’m fairly satisfied with the model selection process I’ve
> been able to conduct, but I would have preferred to test each treatment
> * year + (year|site) model under m100 – m600 conditions and choose the
> optimal one before further comparisons (with simpler model versions).
> The current ad hoc solution is reasonable given the distribution of the
> data, but a true AICc based selection would be easier to “sell” to
> editors, reviewers and readers.
>
> Thanks in advance for any comments or suggestions.
>
> Cheers,
>
> Adjan
>
> PS. I’ve been unable to find a full documentation of the glmmADMB
> package. Only found “Getting started with the glmmADMB package” dated
> Jan. 2, 2012.
I guess I'm wondering what other kind of documentation you
would like other than the help pages in the package (which should
provide a complete, if terse, description of all available options)
and the vignette you describe (which also comes with the package,
and can be accessed via vignette("glmmadmb") when the package is loaded).
I'm open to suggestions about what is missing, and especially open
to contributions!
> Copy of the 8 Feb. posting:
>
> Hi,
>
> I try to evaluate models with various families, with and without zero
> inflation.
>
> This is the set of models I’ve tried so far:
[snip]
> For species A, m100, m300 and m600 run smoothly, but the others all
> return the same error message:
>
>
>
> Error in glmmadmb(resp ~ treatment * year + (year | site) +
> offset(log(bot$area)), :
>
> The function maximizer failed (couldn't find STD file)
>
>
>
> although I try to follow the “Getting started with the glmmADMB package”
> as closely as possible.
> The error message always comes with an additional warning:
>
> In addition: Warning message:
>
> running command 'C:\Windows\system32\cmd.exe /c "C:/Program
> Files/R/R-2.15.2/library/glmmADMB/bin/windows64/glmmadmb.exe" -maxfn 500
> -maxph 5 -noinit -shess' had status 1
>
>
>
> For other species, other models produce this error while models not
> working for species A work fine (each species has the same
> data-structure, but with different amounts and distributions of the
> zeros, and different number of cases).
>
>
>
> The glmmADMB package was downloaded today from
> http://glmmadmb.r-forge.r-project.org/ and I run R x64 2.15.2
> (2012-10-26) under Windows 7 in a university network.
>
> Grateful for any hints or comments that can lead me out of this trench.
>
> Have a nice weekend!
>
> Adjan
>
>
Here's my R code:
## read in files and merge them:
dfiles <- list.files(pattern="[A-Za-z]{2}.txt")
dfiles <- setdiff(dfiles,"Eh.txt") ## leave out Eh (too rare)
spnames <- gsub("\\.txt","",dfiles)
datlist <- lapply(dfiles,read.table,header=TRUE)
datlist2 <- lapply(datlist,
function(x) {
names(x)[ncol(x)] <- "number"
x
})
alldat <- do.call(rbind,datlist2)
alldat$sp <- rep(spnames,sapply(datlist,nrow))
alldat <- transform(alldat,sp=reorder(sp,number/area))
library("plyr")
library("reshape2")
ddply(alldat,"sp",summarise,tot=sum(number))
dcast(ddply(alldat,c("sp","treatment"),summarise,tot=sum(number)),
sp~treatment,value=tot)
library("ggplot2")
theme_set(theme_bw())
library("grid")
zm <- theme(panel.margin=unit(0,"line"))
ggplot(alldat,aes(x=year,y=number/area,colour=treatment))+
geom_line(aes(group=site))+
facet_wrap(~sp)+zm +
scale_y_sqrt() ## expand axes a bit: would like to use log,
## but too many zeroes to make it worth it
## collapse to boxplots (ignore sites/years since it
## sort of looks like noise anyway
ggplot(alldat,aes(x=treatment,y=number/area,colour=treatment))+
geom_boxplot(outlier.colour=NULL)+
facet_wrap(~sp)+zm+scale_y_sqrt()
## "beeswarm" plot for more detail
ggplot(alldat,aes(x=treatment,y=number/area,
fill=treatment,color=treatment))+
geom_dotplot(stackdir="center",binaxis="y",alpha=0.8)+
facet_wrap(~sp)+zm
## look at raw numbers rather than densities, to get a better
## sense of the Poisson/NB issue
ggplot(alldat,aes(x=treatment,y=number,
fill=treatment,color=treatment))+
geom_dotplot(stackdir="center",binaxis="y",alpha=0.8)+
facet_wrap(~sp)+zm
## Based on these pictures and summaries, (a) I see very little
## pattern overall (b) some of the data sets are pretty clearly too
## small to do much with (e.g. in the first 7 species, there are
## at most 7 non-zero observations in the "During" phase); but
## let's see how far we can get anyway.
## Try each species with all three families, with/without zero inflation,
## with/without the shess/noninit trick
form <- number ~ treatment * year + (year|site)+ offset(log(area))
library("glmmADMB")
fvec <- c("poisson","nbinom","nbinom1")
params <- expand.grid(family=fvec,zi=c(FALSE,TRUE),sp=levels(alldat$sp),
shess=c(TRUE,FALSE))
nmodels <- nrow(params)
mresults <- vector(nmodels,mode="list")
i <- 1
for (sp in levels(alldat$sp)) {
for (zi in c(FALSE,TRUE)) {
for (fam in fvec) {
for (shess in c(TRUE,FALSE)) {
cat(i,fam,zi,sp,shess,"\n")
mresults[[i]] <-
try(glmmadmb(form,family=fam,zeroInflation=zi,
data=alldat[alldat$sp==sp,], ## watch out for
subset()!
admb.opts=if (shess) admbControl() else
admbControl(shess=FALSE,noinit=FALSE)),
silent=TRUE)
i <- i+1
}
}
}
}
save("mresults",file="mresults.RData")
params$OK <- sapply(mresults,inherits,"try-error")
## whether we succeed
ggplot(params,aes(x=family,y=sp,fill=OK))+geom_tile()+
facet_grid(zi~shess,
labeller=label_both)
## whether they can be made to work with either 'shess/noinit' setting
params2 <- dcast(params,family+zi+sp~.,fun.aggregate=any,value.var="OK")
names(params2)[4] <- "OK"
ggplot(params,aes(x=family,y=sp,fill=OK))+geom_tile()+
facet_grid(zi~.,
labeller=label_both)
More information about the R-sig-mixed-models
mailing list