[Lme4-commits] r1710 - in pkg/lme4: R man tests testsx
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 24 15:59:34 CEST 2012
Author: mmaechler
Date: 2012-04-24 15:59:33 +0200 (Tue, 24 Apr 2012)
New Revision: 1710
Modified:
pkg/lme4/R/bootMer.R
pkg/lme4/R/lmer.R
pkg/lme4/R/predict.R
pkg/lme4/man/bootMer.Rd
pkg/lme4/man/cbpp.Rd
pkg/lme4/man/glmer.Rd
pkg/lme4/man/nlmer.Rd
pkg/lme4/tests/glmer-1.R
pkg/lme4/tests/glmmExt.R
pkg/lme4/testsx/testcolonizer.R
Log:
for now, skip the "nAGQ > 1" tests; some \pkg{.} some \link{} in the Roxygen doc strings --> *Rd
Modified: pkg/lme4/R/bootMer.R
===================================================================
--- pkg/lme4/R/bootMer.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/R/bootMer.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -49,7 +49,7 @@
##' display. Default is \code{"none"}.
##' @param PBargs a list of additional arguments to the progress bar function.
##' @return an object of S3 \code{\link{class}} \code{"boot"}, compatible with
-##' \pkg{boot} package's \code{boot()} result.
+##' \pkg{boot} package's \code{\link[boot]{boot}()} result.
##' @seealso For inference, including confidence intervals,
##' \code{\link{profile-methods}}.
##'
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/R/lmer.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -199,9 +199,9 @@
##' for the preliminary (random effects parameters only) optimization, while
##' the second will be used for the final (random effects plus
##' fixed effect parameters) phase. The built-in optimizers are
-##' \code{\link{Nelder_Mead}} and \code{\link{bobyqa}} (from
-##' the \code{minqa} package; the default
-##' is to use \code{\link{bobyqa}} for the first and
+##' \code{\link{Nelder_Mead}} and \code{\link[minqa]{bobyqa}} (from
+##' the \pkg{minqa} package; the default
+##' is to use \code{\link[minqa]{bobyqa}} for the first and
##' \code{\link{Nelder_Mead}} for the final phase.
##' (FIXME: simplify if possible!). For difficult model fits we have found
##' \code{\link{Nelder_Mead}} to be more reliable but occasionally slower than
@@ -217,9 +217,9 @@
##' message, or explanation of convergence failure).
##' Special provisions are made for \code{\link{bobyqa}},
##' \code{\link{Nelder_Mead}}, and optimizers wrapped in
-##' the \code{optimx} package; to use \code{optimx} optimizers
-##' (including \code{L-BFGS-B} from base \code{optim} and
-##' \code{nlminb}), pass the \code{method} argument to \code{optim}
+##' the \pkg{optimx} package; to use \pkg{optimx} optimizers
+##' (including \code{L-BFGS-B} from base \code{\link{optim}} and
+##' \code{\link{nlminb}}), pass the \code{method} argument to \code{optim}
##' in the \code{control} argument.
##'
##' @param mustart optional starting values on the scale of the conditional mean,
@@ -243,15 +243,14 @@
##' xyplot(incidence/size ~ period, group=herd, type="a", data=cbpp)
##' (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
##' data = cbpp, family = binomial))
-##' ## using nAGQ=0L only gets close to the optimum
+##' ## using nAGQ=0 only gets close to the optimum
##' (gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
-##' cbpp, binomial, nAGQ = 0L))
-##' ## using nAGQ=9L provides a better evaluation of the deviance
-##' (gm1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
-##' cbpp, binomial, nAGQ = 9L))
-##' ## check with nAGQ=25L
-##' (gm1c <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
-##' cbpp, binomial, nAGQ = 25L))
+##' cbpp, binomial, nAGQ = 0))
+##' if(FALSE) { ##________ FIXME _______
+##' ## using nAGQ = 9 provides a better evaluation of the deviance
+##' (gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+##' cbpp, binomial, nAGQ = 9))
+##' }#__ end{FIXME} __
##'
##' ## GLMM with individual-level variability (accounting for overdispersion)
##' cbpp$obs <- 1:nrow(cbpp)
@@ -332,7 +331,7 @@
devfun <- function(theta) {
pp$setTheta(theta)
## for consistency start from known mu and weights
- resp$updateMu(lp0)
+ resp$updateMu(lp0)
pwrssUpdate(pp, resp, tol=1e-7)
}
environment(devfun) <- rho
@@ -1118,7 +1117,8 @@
### FIXME: Probably should save the control settings and the optimizer name in the merMod object
opt <- Nelder_Mead(ff, x0, lower=lower, control=control)
mkMerMod(environment(ff), opt, list(flist=object at flist, cnms=object at cnms, Gp=object at Gp,
- lower=object at lower), object at frame, getCall(object))
+ lower=object at lower),
+ object at frame, getCall(object))
}
##' @S3method refitML merMod
Modified: pkg/lme4/R/predict.R
===================================================================
--- pkg/lme4/R/predict.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/R/predict.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -30,7 +30,7 @@
terms=NULL, type=c("link","response"),
allow.new.levels=FALSE, ...) {
## FIXME: appropriate names for result vector?
- if (any(getME(object,"offset")!=0)) stop("offsets not handled yet") ## FIXME
+ if (any(getME(object,"offset")!=0)) stop("offsets not handled yet") ## FIXME for glmer()
type <- match.arg(type)
if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
X_orig <- getME(object, "X")
@@ -57,7 +57,7 @@
## what's the appropriate test?
if (is.language(REform)) {
re <- ranef(object)
- ##
+ ##
ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],levels)
re_x <- mapply(function(x,n) {
@@ -93,5 +93,5 @@
}
return(pred)
}
-
-
+
+
Modified: pkg/lme4/man/bootMer.Rd
===================================================================
--- pkg/lme4/man/bootMer.Rd 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/man/bootMer.Rd 2012-04-24 13:59:33 UTC (rev 1710)
@@ -45,8 +45,8 @@
}
\value{
an object of S3 \code{\link{class}} \code{"boot"},
- compatible with \pkg{boot} package's \code{boot()}
- result.
+ compatible with \pkg{boot} package's
+ \code{\link[boot]{boot}()} result.
}
\description{
Perform model-based (Semi-)parametric bootstrap for mixed
Modified: pkg/lme4/man/cbpp.Rd
===================================================================
--- pkg/lme4/man/cbpp.Rd 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/man/cbpp.Rd 2012-04-24 13:59:33 UTC (rev 1710)
@@ -39,28 +39,18 @@
\examples{
## response as a matrix
(m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- cbpp, binomial, nAGQ=25L))
+ family = binomial, data = cbpp))
## response as a vector of probabilities and usage of argument "weights"
m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
- cbpp, binomial, nAGQ=25L)
+ family = binomial, data = cbpp)
## Confirm that these are equivalent:
stopifnot(all.equal(fixef(m1), fixef(m1p), tol = 1e-5),
- all.equal(ranef(m1), ranef(m1p), tol = 1e-5),
- TRUE)
-## Can this section be moved to a test file? I don't think it belongs in an example. DB
-for(m in c(m1, m1p)) {
- cat("-------\\n\\nCall: ",
- paste(format(getCall(m)), collapse="\\n"), "\\n")
- print(logLik(m)); cat("AIC:", AIC(m), "\\n") ; cat("BIC:", BIC(m),"\\n")
-}
-stopifnot(all.equal(logLik(m1), logLik(m1p), tol = 1e-5),
- all.equal(AIC(m1), AIC(m1p), tol = 1e-5),
- all.equal(BIC(m1), BIC(m1p), tol = 1e-5))
+ all.equal(ranef(m1), ranef(m1p), tol = 1e-5))
+%% more extensive variations of the above --> ../tests/glmer-1.R
## GLMM with individual-level variability (accounting for overdispersion)
cbpp$obs <- 1:nrow(cbpp)
-(m2 <- glmer(cbind(incidence, size - incidence) ~ period +
- (1 | herd) + (1|obs),
+(m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd) + (1|obs),
family = binomial, data = cbpp))
}
\keyword{datasets}
Modified: pkg/lme4/man/glmer.Rd
===================================================================
--- pkg/lme4/man/glmer.Rd 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/man/glmer.Rd 2012-04-24 13:59:33 UTC (rev 1710)
@@ -6,7 +6,7 @@
control = list(), start = NULL, verbose = 0L,
nAGQ = 1L, compDev = TRUE, subset, weights, na.action,
offset, contrasts = NULL, mustart, etastart,
- devFunOnly = FALSE, tolPwrss = 1e-10,
+ devFunOnly = FALSE, tolPwrss = 1e-07,
optimizer = c("bobyqa", "Nelder_Mead"), ...)
}
\arguments{
@@ -43,9 +43,9 @@
parameters only) optimization, while the second will be
used for the final (random effects plus fixed effect
parameters) phase. The built-in optimizers are
- \code{\link{Nelder_Mead}} and \code{\link{bobyqa}} (from
- the \code{minqa} package; the default is to use
- \code{\link{bobyqa}} for the first and
+ \code{\link{Nelder_Mead}} and \code{\link[minqa]{bobyqa}}
+ (from the \pkg{minqa} package; the default is to use
+ \code{\link[minqa]{bobyqa}} for the first and
\code{\link{Nelder_Mead}} for the final phase. (FIXME:
simplify if possible!). For difficult model fits we have
found \code{\link{Nelder_Mead}} to be more reliable but
@@ -62,10 +62,10 @@
\code{message} (informational message, or explanation of
convergence failure). Special provisions are made for
\code{\link{bobyqa}}, \code{\link{Nelder_Mead}}, and
- optimizers wrapped in the \code{optimx} package; to use
- \code{optimx} optimizers (including \code{L-BFGS-B} from
- base \code{optim} and \code{nlminb}), pass the
- \code{method} argument to \code{optim} in the
+ optimizers wrapped in the \pkg{optimx} package; to use
+ \pkg{optimx} optimizers (including \code{L-BFGS-B} from
+ base \code{\link{optim}} and \code{\link{nlminb}}), pass
+ the \code{method} argument to \code{optim} in the
\code{control} argument.}
\item{mustart}{optional starting values on the scale of
@@ -202,15 +202,14 @@
xyplot(incidence/size ~ period, group=herd, type="a", data=cbpp)
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial))
-## using nAGQ=0L only gets close to the optimum
+## using nAGQ=0 only gets close to the optimum
(gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- cbpp, binomial, nAGQ = 0L))
-## using nAGQ=9L provides a better evaluation of the deviance
-(gm1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- cbpp, binomial, nAGQ = 9L))
-## check with nAGQ=25L
-(gm1c <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- cbpp, binomial, nAGQ = 25L))
+ cbpp, binomial, nAGQ = 0))
+if(FALSE) { ##________ FIXME _______
+## using nAGQ = 9 provides a better evaluation of the deviance
+(gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ cbpp, binomial, nAGQ = 9))
+}#__ end{FIXME} __
## GLMM with individual-level variability (accounting for overdispersion)
cbpp$obs <- 1:nrow(cbpp)
Modified: pkg/lme4/man/nlmer.Rd
===================================================================
--- pkg/lme4/man/nlmer.Rd 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/man/nlmer.Rd 2012-04-24 13:59:33 UTC (rev 1710)
@@ -107,9 +107,9 @@
parameters only) optimization, while the second will be
used for the final (random effects plus fixed effect
parameters) phase. The built-in optimizers are
- \code{\link{Nelder_Mead}} and \code{\link{bobyqa}} (from
- the \code{minqa} package; the default is to use
- \code{\link{bobyqa}} for the first and
+ \code{\link{Nelder_Mead}} and \code{\link[minqa]{bobyqa}}
+ (from the \pkg{minqa} package; the default is to use
+ \code{\link[minqa]{bobyqa}} for the first and
\code{\link{Nelder_Mead}} for the final phase. (FIXME:
simplify if possible!). For difficult model fits we have
found \code{\link{Nelder_Mead}} to be more reliable but
@@ -126,10 +126,10 @@
\code{message} (informational message, or explanation of
convergence failure). Special provisions are made for
\code{\link{bobyqa}}, \code{\link{Nelder_Mead}}, and
- optimizers wrapped in the \code{optimx} package; to use
- \code{optimx} optimizers (including \code{L-BFGS-B} from
- base \code{optim} and \code{nlminb}), pass the
- \code{method} argument to \code{optim} in the
+ optimizers wrapped in the \pkg{optimx} package; to use
+ \pkg{optimx} optimizers (including \code{L-BFGS-B} from
+ base \code{\link{optim}} and \code{\link{nlminb}}), pass
+ the \code{method} argument to \code{optim} in the
\code{control} argument.}
}
\description{
Modified: pkg/lme4/tests/glmer-1.R
===================================================================
--- pkg/lme4/tests/glmer-1.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/tests/glmer-1.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -37,20 +37,43 @@
## now
#bobyqa(m1e, control = list(iprint = 2L))
-m0 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- family = binomial, data = cbpp, verbose = 2L)
+(m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, verbose = 2L)
+## response as a vector of probabilities and usage of argument "weights"
+m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
+ family = binomial, data = cbpp, verbose = 2L)
+## Confirm that these are equivalent:
+stopifnot(all.equal(fixef(m1), fixef(m1p), tol = 1e-5),
+ all.equal(ranef(m1), ranef(m1p), tol = 1e-5),
+ TRUE)
+for(m in c(m1, m1p)) {
+ cat("-------\\n\\nCall: ",
+ paste(format(getCall(m)), collapse="\\n"), "\\n")
+ print(logLik(m)); cat("AIC:", AIC(m), "\\n") ; cat("BIC:", BIC(m),"\\n")
+}
+stopifnot(all.equal(logLik(m1), logLik(m1p), tol = 1e-5),
+ all.equal(AIC(m1), AIC(m1p), tol = 1e-5),
+ all.equal(BIC(m1), BIC(m1p), tol = 1e-5))
-m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- optimizer="bobyqa",
- family = binomial, data = cbpp, verbose = 2L,
- control = list(rhobeg=0.2, rhoend=2e-7), tolPwrss=1e-8)
+m1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ optimizer="bobyqa",
+ family = binomial, data = cbpp, verbose = 2L,
+ control = list(rhobeg=0.2, rhoend=2e-7), tolPwrss=1e-8)
+
+if(FALSE) { ##_____________ FIXME _____________ not yet nAGQ > 1 ______________
+
+## using nAGQ=9L provides a better evaluation of the deviance
+m.9 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 9)
+
+## check with nAGQ = 25
m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- family = binomial, data = cbpp, nAGQ=25L)
+ family = binomial, data = cbpp, nAGQ = 25)
## loosened tolerance on parameters
-stopifnot(is((cm1 <- coef(m2)), "coef.mer"),
- dim(cm1$herd) == c(15,4),
+stopifnot(is((cm2 <- coef(m2)), "coef.mer"),
+ dim(cm2$herd) == c(15,4),
all.equal(fixef(m2),
### lme4a [from an Ubuntu 11.10 amd64 system]
### c(-1.39922533406847, -0.991407294757321,
@@ -63,19 +86,22 @@
## with bobyqa first (AGQ=0), then
all.equal(deviance(m2), 101.119749563, tol=1e-9)
)
+}##_____________ end{FIXME} _____________ not yet nAGQ > 1 ______________
-stopifnot(is((cm1 <- coef(m1)), "coef.mer"),
+
+stopifnot(is((cm1 <- coef(m1b)), "coef.mer"),
dim(cm1$herd) == c(15,4),
- all.equal(fixef(m1),
+ all.equal(fixef(m1b),
## these values are those of "old-lme4":
## c(-1.39853504914, -0.992334711,
## -1.12867541477, -1.58037390498),
## lme4[r 1636], 64-bit ubuntu 11.10:
c(-1.3788385, -1.0589543,
-1.1936382, -1.6306271),
- tol = 1.e-3,
+ tol = 1e-3,
check.attr=FALSE)
)
+## FIXME --- compare m1b with m1 and m0 ---
## Deviance for the new algorithm is lower, eventually we should change the previous test
Modified: pkg/lme4/tests/glmmExt.R
===================================================================
--- pkg/lme4/tests/glmmExt.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/tests/glmmExt.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -1,7 +1,7 @@
library(lme4)
## generate a basic Gamma/random effects sim
set.seed(101)
-d <- expand.grid(block=LETTERS[1:26],rep=1:100)
+d <- expand.grid(block=LETTERS[1:26], rep=1:100, KEEP.OUT.ATTRS = FALSE)
d$x <- runif(nrow(d)) ## sd=1
reff_f <- rnorm(length(levels(d$block)),sd=1)
## need intercept large enough to avoid negative values
@@ -10,6 +10,7 @@
## inverse link
d$mu <- 1/d$eta
d$y <- rgamma(nrow(d),scale=d$mu/2,shape=2)
+str(d)## 2600 obs .. 'block' with 26 levels
## version with Poisson errors instead (to check GLM in general)
dP <- d
@@ -30,30 +31,28 @@
dGi$y <- rnorm(nrow(d),dGi$mu,sd=0.01)
############
-gm0 <- glm(y~1, data=d, family=Gamma)
-gm1 <- glm(y~block-1, data=d, family=Gamma)
-sd(coef(gm1))
+gm0 <- glm(y ~ 1, data=d, family=Gamma)
+gm1 <- glm(y ~ block-1, data=d, family=Gamma)
+sd(coef(gm1)) # 1.007539
-gm2 <- glmer(y ~ 1 + (1|block), d, Gamma,
- verbose=TRUE)
+gm2 <- glmer(y ~ 1 + (1|block), d, Gamma, verbose = 4)
## do we do any better with a correctly specified model??
## no.
-gm3 <- glmer(y ~ x + (1|block), d, Gamma,
- verbose=TRUE)
+gm3 <- glmer(y ~ x + (1|block), d, Gamma, verbose = 4)
## correctly specified model with "true" parameters as starting values
gm3B <- glmer(y ~ x + (1|block), d, Gamma,
start=list(fixef=c(4,3),ST=list(matrix(1))),
- verbose=TRUE)
+ verbose = 4)
-stopifnot(all.equal(fixef(gm3),fixef(gm3B)),
+stopifnot(all.equal(fixef (gm3),fixef (gm3B)),
all.equal(VarCorr(gm3),VarCorr(gm3B)))
###
## Poisson
-gP1 <- glmer(y ~ 1 + (1|block), data=dP, family=poisson, verbose=TRUE)
-gP2 <- glmer(y ~ x + (1|block), data=dP, family=poisson, verbose=TRUE)
+gP1 <- glmer(y ~ 1 + (1|block), data=dP, family=poisson, verbose= 2)
+gP2 <- glmer(y ~ x + (1|block), data=dP, family=poisson, verbose= 2)
## works just fine.
@@ -62,11 +61,14 @@
gG2 <- glmer(y ~ x + (1|block), data=dG, family=gaussian(link="log"), verbose=TRUE)
+if(Sys.info()["user"] != "maechler") { # <- seg.faults (MM)
+
## FIXME: get these working
## Gaussian with inverse link ... FIXME: not working yet
-try(gGi1 <- glmer(y ~ 1 + (1|block), data=dGi, family=gaussian(link="inverse"), verbose=TRUE))
-try(gGi2 <- glmer(y ~ x + (1|block), data=dGi, family=gaussian(link="inverse"), verbose=TRUE))
+try(gGi1 <- glmer(y ~ 1 + (1|block), data=dGi, family=gaussian(link="inverse"), verbose= 3))
+try(gGi2 <- glmer(y ~ x + (1|block), data=dGi, family=gaussian(link="inverse"), verbose= 3))
## sets variance to zero, converges back to GLM solution
+}
## FIXME: cloglog link
Modified: pkg/lme4/testsx/testcolonizer.R
===================================================================
--- pkg/lme4/testsx/testcolonizer.R 2012-04-23 14:52:20 UTC (rev 1709)
+++ pkg/lme4/testsx/testcolonizer.R 2012-04-24 13:59:33 UTC (rev 1710)
@@ -1,15 +1,15 @@
library(lme4.0)
-load("colonizer_rand.RData")
-m1 <- glmer(form1,data=randdat,
- family=poisson)
-m2 <- glmer(form2,data=randdat,
- family=poisson)
-detach("package:lme4.0")
+## Emacs M-<Enter> --> setwd() correctly
+(load("colonizer_rand.RData"))
+summary(m1.0 <- glmer(form1, data=randdat, family=poisson))
+summary(m2.0 <- glmer(form2, data=randdat, family=poisson))
+
+detach("package:lme4.0", unload=TRUE)
library(lme4)
-try(m1 <- glmer(form1,data=randdat,
- family=poisson,verbose=10L))
-try(m1 <- glmer(form1,data=randdat,
- family=poisson,verbose=10L,tolPwrss=1e-13))
-try(m2 <- glmer(form2,data=randdat,
- family=poisson,verbose=10L))
+packageDescription("lme4")
+## currently (r 1704), all give "Downdated VtV is not positive definite"
+try(m1 <- glmer(form1,data=randdat, family=poisson, verbose=10L))
+try(m1 <- glmer(form1,data=randdat, family=poisson, verbose=10L, tolPwrss=1e-13))
+try(m2 <- glmer(form2,data=randdat, family=poisson, verbose=10L))
+
More information about the Lme4-commits
mailing list