[Lme4-commits] r1652 - in pkg/lme4Eigen: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 14 14:24:11 CET 2012
Author: bbolker
Date: 2012-03-14 14:24:11 +0100 (Wed, 14 Mar 2012)
New Revision: 1652
Added:
pkg/lme4Eigen/tests/profile.R
Modified:
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/profile.R
pkg/lme4Eigen/man/mcmcsamp.Rd
Log:
add 'upper' to optwrap(), adjust lmer.R calls accordingly
make profiling use optimizer choice/optwrap (default back to bobyqa)
expose plot.merMod
remove bogus second optimization call
documentation tweak to mcmcsamp
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-03-13 22:38:25 UTC (rev 1651)
+++ pkg/lme4Eigen/NAMESPACE 2012-03-14 13:24:11 UTC (rev 1652)
@@ -117,6 +117,7 @@
S3method(nobs,merMod)
S3method(plot,coef.mer)
S3method(plot,lmList.confint)
+S3method(plot,merMod)
S3method(plot,ranef.mer)
S3method(predict,merMod)
S3method(print,merMod)
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-03-13 22:38:25 UTC (rev 1651)
+++ pkg/lme4Eigen/R/lmer.R 2012-03-14 13:24:11 UTC (rev 1652)
@@ -136,7 +136,8 @@
if (devFunOnly) return(devfun)
opt <- optwrap(optimizer,
- devfun, rho$pp$theta, lower=reTrms$lower, control=control, rho, adj=FALSE)
+ devfun, rho$pp$theta, lower=reTrms$lower, control=control,
+ rho=rho, adj=FALSE)
mkMerMod(environment(devfun), opt, reTrms, fr, mc)
}## { lmer }
@@ -347,7 +348,8 @@
if (length(optimizer)==1) {
optimizer <- replicate(2,optimizer)
}
- opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower, control, rho,
+ opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower,
+ control=control, rho=rho,
adj=FALSE, verbose=verbose)
rho$nAGQ <- nAGQ
@@ -364,7 +366,8 @@
devfun <- mkdevfun(rho, nAGQ)
if (devFunOnly) return(devfun)
- opt <- optwrap(optimizer[[2]],devfun,c(rho$pp$theta, rho$beta0), rho$lower, control, rho,
+ opt <- optwrap(optimizer[[2]],devfun,c(rho$pp$theta, rho$beta0),
+ rho$lower, control=control, rho=rho,
adj=TRUE, verbose=verbose)
}
@@ -430,7 +433,8 @@
optimizer <- replicate(2,optimizer)
}
- opt <- optwrap(optimizer[[1]],devfun, rho$pp$theta, rho$lower, control, rho, adj=FALSE)
+ opt <- optwrap(optimizer[[1]],devfun, rho$pp$theta, rho$lower,
+ control=control, rho=rho, adj=FALSE)
if (nAGQ > 0L) {
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
@@ -445,7 +449,8 @@
devfun <- mkdevfun(rho, nAGQ)
if (devFunOnly) return(devfun)
- opt <- optwrap(optimizer[[2]],devfun, par=c(rho$pp$theta, rho$beta0), lower=rho$lower, control=control, rho,
+ opt <- optwrap(optimizer[[2]],devfun, par=c(rho$pp$theta, rho$beta0),
+ lower=rho$lower, control=control, rho=rho,
adj=TRUE, verbose=verbose)
@@ -1908,9 +1913,16 @@
optfun
}
-optwrap <- function(optimizer,fn,par,lower,control,rho,adj, verbose=0L) {
+optwrap <- function(optimizer, fn, par, lower=-Inf, upper=Inf,
+ control=list(), rho, adj=FALSE, verbose=0L) {
+ ## control and rho must be specified if adj==TRUE;
+ ## otherwise this is a fairly simple wrapper
optfun <- getOptfun(optimizer)
## special-purpose control parameter tweaks: only for second round in nlmer, glmer
+
+ lower <- rep(lower, length.out=length(par))
+ upper <- rep(upper, length.out=length(par))
+
if (adj && is.character(optimizer))
switch(optimizer,
bobyqa = {
@@ -1927,12 +1939,13 @@
if (is.null(control$xt)) control$xt <- control$xst*5e-4
rho$control <- control
})
- arglist <- list(fn=fn, par=par, lower=lower, control=control)
+ arglist <- list(fn=fn, par=par, lower=lower, upper=upper, control=control)
## optimx: must pass method in control (?) because 'method' was previously
## used in lme4 to specify REML vs ML
## FIXME: test -- does deparse(substitute(...)) clause work?
if (optimizer=="optimx" || deparse(substitute(optimizer))=="optimx") {
- if (is.null(method <- control$method)) stop("must specify 'method' explicitly for optimx")
+ if (is.null(method <- control$method))
+ stop("must specify 'method' explicitly for optimx")
arglist$control$method <- NULL
arglist <- c(arglist,list(method=method))
}
Modified: pkg/lme4Eigen/R/profile.R
===================================================================
--- pkg/lme4Eigen/R/profile.R 2012-03-13 22:38:25 UTC (rev 1651)
+++ pkg/lme4Eigen/R/profile.R 2012-03-14 13:24:11 UTC (rev 1652)
@@ -34,15 +34,18 @@
profile.merMod <- function(fitted, alphamax = 0.01, maxpts = 100, delta = cutoff/8,
## tr = 0, ## FIXME: remove if not doing anything ...
verbose=0, devtol=1e-9,
- startval = "prev", ...) {
+ startval = "prev",
+ optimizer="bobyqa", ...) {
- ## FIXME: allow choice of optimizer? choice of nextstep/nextstart algorithm?
+ ## FIXME: allow choice of nextstep/nextstart algorithm?
+ ## FIXME: by default, get optimizer from within fitted object
## FIXME: allow selection of individual variables to profile (by location and?? name)
## FIXME: allow for failure of bounds (non-pos-definite correlation matrices) when >1 cor parameter
## FIXME: generalize to GLMMs
## (use different devfun;
## be careful with scale parameter;
## profile all parameters at once rather than RE first and then fixed)
+
dd <- devfun2(fitted)
base <- attr(dd, "basedev")
@@ -141,15 +144,10 @@
lowcut <- lower[w]
upcut <- upper[w]
zeta <- function(xx,start) {
- ## FIXME: use Nelder_Mead, or (for GLMMs) optimizer choice??
- ores <- bobyqa(start,
- function(x) dd(mkpar(nvp, w, xx, x)),
- lower = lowvp[-w],
- upper = upvp[-w])
- ores2 <- Nelder_Mead(par=start,
- fn=function(x) dd(mkpar(nvp, w, xx, x)),
- lower = lowvp[-w],
- upper = upvp[-w])
+ ores <- optwrap(optimizer, par=start,
+ fn=function(x) dd(mkpar(nvp, w, xx, x)),
+ lower = lowvp[-w],
+ upper = upvp[-w])
devdiff <- ores$fval-base
if (is.na(ores$fval)) {
@@ -216,8 +214,10 @@
rr$setOffset(Xw * fw + offset.orig)
rho <- as.environment(list(pp=pp1, resp=rr))
parent.env(rho) <- parent.frame()
- ores <- bobyqa(thopt, mkdevfun(rho, 0L), lower = pmax(fitted at lower, -1.0),
- upper = ifelse(fitted at lower==0,Inf,1.0))
+ ores <- optwrap(optimizer,
+ par=thopt, fn=mkdevfun(rho, 0L),
+ lower = pmax(fitted at lower, -1.0),
+ upper = ifelse(fitted at lower==0,Inf,1.0))
fv <- ores$fval
sig <- sqrt((rr$wrss() + pp1$sqrL(1))/n)
c(sign(fw - est) * sqrt(fv - base),
Modified: pkg/lme4Eigen/man/mcmcsamp.Rd
===================================================================
--- pkg/lme4Eigen/man/mcmcsamp.Rd 2012-03-13 22:38:25 UTC (rev 1651)
+++ pkg/lme4Eigen/man/mcmcsamp.Rd 2012-03-14 13:24:11 UTC (rev 1652)
@@ -12,7 +12,7 @@
\item{verbose}{should verbose output be given?}
- \item{...}{}
+ \item{...}{additional arguments}
}
\value{
a Markov chain Monte Carlo sample as a matrix
Added: pkg/lme4Eigen/tests/profile.R
===================================================================
--- pkg/lme4Eigen/tests/profile.R (rev 0)
+++ pkg/lme4Eigen/tests/profile.R 2012-03-14 13:24:11 UTC (rev 1652)
@@ -0,0 +1,39 @@
+library(lme4Eigen)
+
+fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
+
+## 0.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
+##
+system.time( tpr <- profile(fm01ML) )
+
+(confint(tpr) -> CIpr)
+print(xyplot(tpr))
+## comparing against lme4a reference values -- but lme4Eigen returns sigma
+## rather than log(sigma)
+stopifnot(dim(CIpr) == c(3,2),
+ all.equal(unname(CIpr[".sigma",]),exp(c(3.64362, 4.21446)), tol=1e-6),
+ all.equal(unname(CIpr["(Intercept)",]),c(1486.451500,1568.548494)))
+
+## 2D profiles
+fm2ML <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin, REML=0)
+pr2 <- profile(fm2ML)
+(confint(pr2) -> CIpr2)
+
+lme4a_CIpr2 <-
+structure(c(0.633565787613112, 1.09578224011285, -0.721864513060904,
+21.2666273835452, 1.1821039843372, 3.55631937954106, -0.462903300019305,
+24.6778176174587), .Dim = c(4L, 2L), .Dimnames = list(c(".sig01",
+".sig02", ".lsig", "(Intercept)"), c("2.5 %", "97.5 %")))
+lme4a_CIpr2[".lsig",] <- exp(lme4a_CIpr2[".lsig",])
+
+stopifnot(all.equal(unname(CIpr2),unname(lme4a_CIpr2)))
+
+print(xyplot(pr2, absVal=0, aspect=1.3, layout=c(4,1)))
+print(splom(pr2))
+
+## NOT RUN: takes ~ 30 seconds on my machine ...
+fm3ML <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=FALSE)
+## system.time(pr3 <- profile(fm3ML))
+## xyplot(pr3)
+## print(splom(pr3))
+
More information about the Lme4-commits
mailing list