library(lme4)
n <- nrow(sleepstudy)
op <- options(warn = 1, # show as they happen ("false" convergence warnings)
useFancyQuotes = FALSE)
if (.Platform$OS.type != "windows") {
##' remove all attributes but names
dropA <- function(x) `attributes<-`(x, list(names = names(x)))
##' transform result of "numeric" all.equal.list() to a named vector
all.eqL <- function(x1, x2, ...) {
r <- sub("^Component ", '', all.equal(x1, x2, tolerance = 0, ...))
r <- strsplit(sub(": Mean relative difference:", "&&", r),
split="&&", fixed=TRUE)
setNames(as.numeric(vapply(r, `[`, "1.234", 2L)),
## drop surrounding "..."
nm = sub('"$', '', substring(vapply(r, `[`, "nam", 1L), first=2)))
}
seedF <- function(s) {
if(s %in% c(6, 39, 52, 57, 63, 74, 76, 86))
switch(as.character(s)
, "52"=, "63"=, "74" = 2
, "6"=, "39" = 3
, "86" = 8 # needs 4 on Lnx-64b
, "76" = 70 # needs 42 on Lnx-64b
, "57" = 90 # needs 52 on Lnx-64b
)
else if(s %in% c(1, 12, 15, 34, 36, 41, 42, 43, 49, 55, 59, 67, 80, 85)) ## seeds 41,59, .. 15
1.0
else ## seeds 22, 20, and better
0.25
}
## be fast, running only 10 seeds by default:
sMax <- if(lme4:::testLevel() > 1) 99L else 9L
mySeeds <- 0L:sMax
lapply(setNames(,mySeeds), function(seed) {
cat("\n------ random seed =", seed, "---------\n")
set.seed(seed)
v <- rpois(n,1) + 1
w <- 1/v
cat("weights w:\n")
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, REML=FALSE, weights = w); cat("..2:\n")
fm2 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML=FALSE, weights = w)
cat("weights w*10:\n")
fm1.10 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, REML=FALSE, weights = w*10);cat("..2:\n")
fm2.10 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML=FALSE, weights = w*10)
##
ano12... <- dropA(anova(fm1, fm2 ))
ano12.10 <- dropA(anova(fm1.10, fm2.10))
print(aEQ <- all.eqL(ano12..., ano12.10)) # showing differences
if(!exists("notChisq"))
notChisq <<-
local({ n <- names(ano12...)
grep("Chisq", n, value=TRUE, fixed=TRUE, invert=TRUE) })
stopifnot(
all.equal(ano12...$Chisq,
ano12.10$Chisq, tol = 1e-6 * seedF(seed))
,
all.equal(ano12...[notChisq],
ano12.10[notChisq], tol= 1.5e-8 * seedF(seed))
)
aEQ
}) -> rallEQ
cat("=====================================\n")
rallEQ <- t(simplify2array(rallEQ))
notChisq <- intersect(notChisq, colnames(rallEQ))
## sort according to "severity":
srallEQ <- rallEQ[with(as.data.frame(rallEQ), order(AIC, Chisq)), ]
round(log10(srallEQ), 2)
saveRDS(srallEQ, "priorWeightsMod_relerr.rds")
if(!dev.interactive(orNone=TRUE)) pdf("priorWeightsMod_relerr.pdf")
matplot(mySeeds, log10(srallEQ), type="l", xlab=NA) ; grid()
legend("topleft", ncol=3, bty="n",
paste(1:6, colnames(srallEQ), sep = ": "), col=1:6, lty=1:6)
tolD <- sqrt(.Machine$double.eps) # sqrt(eps_C)
abline(h = log10(tolD), col = "forest green", lty=3)
axis(4, at=log10(tolD), label=quote(sqrt(epsilon[c])), las=1)
LRG <- which(srallEQ[,"AIC"] > tolD)
if (length(LRG)>0) {
text(LRG, log10(srallEQ[LRG, "AIC"]), names(LRG), cex = .8)
}
## how close are we ..
str(tF <- sapply(mySeeds, seedF))
round(sort( rallEQ[, "Chisq"] / (tF * 1e-6 ), decreasing=TRUE), 1)
round(sort(apply(rallEQ[,notChisq] / (tF * 1.5e-8), 1, max), decreasing=TRUE), 1)
} ## skip on windows (for speed)
options(op)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.