Nothing
library("lme4")
L <- load(system.file("testdata","crabs_randdata2.Rda",package="lme4"))
## randdata0: simulated data, in form suitable for plotting
## randdata: simulated data, in form suitable for analysis
## fr ## alive/dead formula
## fr2 ## proportion alive formula (use with weights=initial.snail.density)
## FIXME: there are still bigger differences than I'd like between the approaches
## (mostly in the random-effects correlation). It's not clear who's right;
## lme4 thinks its parameters are better, but ?? Could be explored further.
if (FALSE) {
## library(ggplot2) ## commented to avoid triggering Suggests: requirement
library(grid)
zmargin <- theme(panel.margin=unit(0,"lines"))
theme_set(theme_bw())
g1 <- ggplot(randdata0,aes(x=snail.size,y=surv,colour=snail.size,fill=snail.size))+
geom_hline(yintercept=1,colour="black")+
stat_sum(aes(size=factor(..n..)),alpha=0.6)+
facet_grid(.~ttt)+zmargin+
geom_boxplot(fill=NA,outlier.colour=NULL,outlier.shape=3)+ ## set outliers to same colour as points
## (hard to see which are outliers, but it doesn't really matter in this case)
scale_size_discrete("# obs",range=c(2,5))
}
if (.Platform$OS.type != "windows") {
t1 <- system.time(glmer1 <- glmer(fr2,weights=initial.snail.density,
family ="binomial", data=randdata))
t1B <- system.time(glmer1B <- glmer(fr,family ="binomial", data=randdata))
res1 <- c(fixef(glmer1),c(VarCorr(glmer1)$plot))
res1B <- c(fixef(glmer1B),c(VarCorr(glmer1B)$plot))
p1 <- unlist(getME(glmer1,c("theta","beta")))
stopifnot(all.equal(res1,res1B))
dfun <- update(glmer1,devFunOnly=TRUE)
stopifnot(all.equal(dfun(p1),c(-2*logLik(glmer1))))
##
## library(lme4.0) ## version 0.999999.2 results
## t1_lme4.0 <- system.time(glmer1X <-
## glmer(fr2,weights=initial.snail.density,
## family ="binomial", data=randdata))
## dput(c(fixef(glmer1X),c(VarCorr(glmer1X)$plot)))
## p1X <- c(getME(glmer1X,"theta"),getME(glmer1X,"beta"))
p1X <- c(0.681301656652347, -1.14775239687404, 0.436143018123226,
2.77730476938968, 0.609023583738824, -1.60055813739844, 2.0324468778545,
0.624173873057839, -1.7908793509579, -2.44540201631615, -1.42365990002708,
-2.26780929006268, 0.700928084600075, -1.26220238391029, 0.369024582097804,
3.44325347343035, 2.26400391093108)
stopifnot(all.equal(unname(p1),p1X,tolerance=0.03))
dfun(p1X)
dfun(p1)
## ~ 1.8 seconds elapsed time
lme4.0_res <-
structure(c(2.77730476938968, 0.609023583738824, -1.60055813739844,
2.0324468778545, 0.624173873057839, -1.7908793509579, -2.44540201631615,
-1.42365990002708, -2.26780929006268, 0.700928084600075, -1.26220238391029,
0.369024582097804, 3.44325347343035, 2.26400391093108, 0.464171947357232,
-0.532754465140956, -0.532754465140956, 0.801690946568518),
.Names = c("(Intercept)",
"crab.speciesS", "crab.speciesW", "crab.sizeS", "crab.sizeM",
"snail.sizeS", "crab.speciesS:crab.sizeS", "crab.speciesS:crab.sizeM",
"crab.speciesS:snail.sizeS", "crab.speciesW:snail.sizeS",
"crab.sizeS:snail.sizeS", "crab.sizeM:snail.sizeS",
"crab.speciesS:crab.sizeS:snail.sizeS",
"crab.speciesS:crab.sizeM:snail.sizeS", "", "", "", ""))
stopifnot(all.equal(res1,lme4.0_res,tolerance=0.015))
## library("glmmADMB")
## prop/weights formulation: ~ 7 seconds
## t1_glmmadmb <- system.time(glmer1B <- glmmadmb(fr,family ="binomial",
## corStruct="full",data=randdata))
## dput(c(fixef(glmer1B),c(VarCorr(glmer1B)$plot)))
glmmADMB_res <- structure(c(2.7773101267224, 0.609026276823218, -1.60055704634712,
2.03244174458562, 0.624171008585953, -1.79088398816641, -2.44540300134182,
-1.42366043619683, -2.26780858382505, 0.700927141726545, -1.26219964572264,
0.369029052442189, 3.44326297908383, 2.26403738918967, 0.46417,
-0.53253, -0.53253, 0.80169), .Names = c("(Intercept)", "crab.speciesS",
"crab.speciesW", "crab.sizeS", "crab.sizeM", "snail.sizeS", "crab.speciesS:crab.sizeS",
"crab.speciesS:crab.sizeM", "crab.speciesS:snail.sizeS", "crab.speciesW:snail.sizeS",
"crab.sizeS:snail.sizeS", "crab.sizeM:snail.sizeS", "crab.speciesS:crab.sizeS:snail.sizeS",
"crab.speciesS:crab.sizeM:snail.sizeS", "", "", "", ""))
stopifnot(all.equal(res1B,glmmADMB_res,tolerance=0.015))
} ## skip on windows (for speed)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.