[Lme4-commits] r1498 - in pkg/lme4Eigen: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jan 8 00:35:27 CET 2012
Author: dmbates
Date: 2012-01-08 00:35:27 +0100 (Sun, 08 Jan 2012)
New Revision: 1498
Modified:
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/profile.R
pkg/lme4Eigen/tests/lmer-1.R
Log:
Nlmer now working (for nAGQ=1). Updated tests and got profile working again.
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-01-07 17:29:59 UTC (rev 1497)
+++ pkg/lme4Eigen/R/lmer.R 2012-01-07 23:35:27 UTC (rev 1498)
@@ -1,3 +1,31 @@
+##' Fit a linear mixed model or a generalized linear mixed model or a nonlinear mixed model.
+##'
+##'
+##'
+##' @title Fit a linear mixed-effects model
+##' @param formula a two-sided linear formula object describing the
+##' fixed-effects part of the model, with the response on the left of a
+##' \code{~} operator and the terms, separated by \code{+} operators, on
+##' the right. The vertical bar character \code{"|"} separates an
+##' expression for a model matrix and a grouping factor.
+##' @param data an optional data frame containing the variables named in
+##' \code{formula}. By default the variables are taken from the
+##' environment from which \code{lmer} is called.
+##' @param REML logical - Should the estimates be chosen to optimize the
+##' REML criterion (as opposed to the log-likelihood)? Defaults to \code{TRUE}.
+##' @param sparseX logical - should a sparse model matrix be used for the
+##' fixed-effects terms? Defaults to \code{FALSE}.
+##' @param control
+##' @param start
+##' @param verbose
+##' @param subset
+##' @param weights
+##' @param na.action
+##' @param offset
+##' @param contrasts
+##' @param devFunOnly
+##' @param ...
+##' @return
lmer <- function(formula, data, REML = TRUE, sparseX = FALSE,
control = list(), start = NULL,
verbose = 0L, subset, weights, na.action, offset,
@@ -154,7 +182,6 @@
if (family$family %in% c("quasibinomial", "quasipoisson", "quasi"))
stop('"quasi" families cannot be used in glmer')
- if (nAGQ < 0L || nAGQ > 25L) warning("nAGQ must be in the range [0, 25]")
stopifnot(length(nAGQ <- as.integer(nAGQ)) == 1L,
nAGQ >= 0L,
nAGQ <= 25L,
@@ -482,31 +509,6 @@
}
}
-## setMethod("show", signature("lmerResp"), function(object)
-## {
-## with(object,
-## print(head(cbind(weights, offset, mu, y, sqrtrwt, wtres, sqrtXwt))))
-## })
-
-## setMethod("show", signature("glmerResp"), function(object)
-## {
-## with(object,
-## print(head(cbind(weights, offset, eta, mu, y, muEta,
-## variance, sqrtrwt, wtres, sqrtXwt,
-## sqrtWrkWt, wrkResids, wrkResp))))
-## })
-
-## setMethod("show", signature("Rcpp_reModule"), function(object)
-## {
-## with(object, print(head(cbind(u0, incr, u))))
-## })
-
-
-## setMethod("show", signature("Rcpp_deFeMod"), function(object)
-## {
-## with(object, print(cbind(coef0, incr, coef, Vtr)))
-## })
-
## create a deviance evaluation function that uses the sigma parameters
## df2 <- function(dd) {
## stopifnot(is.function(dd),
@@ -667,38 +669,31 @@
##' @param data an optional data frame containing the variables named in
##' \code{formula}. By default the variables are taken from the
##' environment from which \code{nlmer} is called.
-##' @param family
##' @param start starting estimates for the nonlinear model
##' parameters, as a named numeric vector
##' @param verbose integer scalar passed to nlminb. If negative then
##' diagnostic output from the PIRLS (penalized iteratively
##' reweighted least squares) step is also provided.
##' @param nAGQ number of adaptive Gauss-Hermite quadrature points to use
-##' @param doFit logical scalar. If FALSE the optimization
-##' environment is returned. Otherwise the parameters are estimated
-##' and an object of S4 class "mer" is returned.
##' @param subset further model specifications as in
##' \code{\link[stats]{lm}}; see there for details.
-##' @param weights further model specifications as in
+##' @param weights further model specifications as in
##' \code{\link[stats]{lm}}; see there for details.
-##' @param na.action further model specifications as in
+##' @param na.action further model specifications as in
##' \code{\link[stats]{lm}}; see there for details.
-##' @param mustart
-##' @param etastart
-##' @param sparseX
-##' @param contrasts further model specifications as in
+##' @param offset
+##' @param contrasts further model specifications as in
##' \code{\link[stats]{lm}}; see there for details.
-##' @param control a list of control parameters passed to bobyqa.
-##' @param ...
-
-##' @return an object of S4 class "merMod"
-nlmer <- function(formula, data, family = gaussian, start = NULL,
- verbose = 0L, nAGQ = 1L, doFit = TRUE,
- subset, weights, na.action, offset,
- contrasts = NULL, devFunOnly=0L,
- tolPwrss = 1e-10, optimizer=c("bobyqa","NelderMead"), ...)
+##' @param devFunOnly
+##' @param tolPwrss
+##' @param optimizer
+##' @param ...
+##' @param control a list of control parameters passed to the optimizer.
+nlmer <- function(formula, data, control = list(), start = NULL, verbose = 0L,
+ nAGQ = 1L, subset, weights, na.action, offset,
+ contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-10,
+ optimizer=c("NelderMead","bobyqa"), ...)
{
- if (!missing(family)) stop("code not yet written")
mf <- mc <- match.call()
m <- match(c("data", "subset", "weights", "na.action",
"offset", "etastart", "mustart"),
@@ -763,10 +758,74 @@
rho$beta0 <- rho$pp$beta0
rho$lower <- reTrms$lower
devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
- devfun
-}
+ if (devFunOnly && !nAGQ) return(devfun)
+ control$iprint <- min(verbose, 3L)
+ opt <- bobyqa(rho$pp$theta, devfun, rho$lower, control=control)
+ if (nAGQ > 0L) {
+ rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
+ rho$u0 <- rho$pp$u0
+ rho$beta0 <- rho$pp$beta0
+ rho$dpars <- seq_along(rho$pp$theta)
+ if (nAGQ > 1L) {
+ if (length(reTrms$flist) != 1L || length(reTrms$cnms[[1]]) != 1L)
+ stop("nAGQ > 1 is only available for models with a single, scalar random-effects term")
+ rho$fac <- reTrms$flist[[1]]
+ }
+ devfun <- mkdevfun(rho, nAGQ)
+ if (devFunOnly) return(devfun)
+ opt <- switch(match.arg(optimizer),
+ bobyqa = {
+ if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
+ if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
+ rho$control <- control
+ bobyqa(c(rho$pp$theta, rho$beta0), devfun, rho$lower, control=control)
+ },
+ NelderMead = {
+ xst <- c(rep.int(0.1, length(rho$dpars)), sqrt(diag(environment(devfun)$pp$unsc())))
+ nM <- NelderMead$new(lower=rho$lower, upper=rep.int(Inf, length(rho$lower)), xst=0.2*xst,
+ x0=with(environment(devfun), c(pp$theta, pp$beta0)),
+ xt=xst*0.0001)
+ cc <- do.call(function(iprint=0L, maxfun=10000L, FtolAbs=1e-5,
+ FtolRel=1e-15, XtolRel=1e-7,
+ MinfMax=.Machine$double.xmin) {
+ if (length(list(...))>0) warning("unused control arguments ignored")
+ list(iprint=iprint, maxfun=maxfun, FtolAbs=FtolAbs, FtolRel=FtolRel,
+ XtolRel=XtolRel, MinfMax=MinfMax)
+ }, control)
+ nM$setMaxeval(cc$maxfun)
+ nM$setFtolAbs(cc$FtolAbs)
+ nM$setFtolRel(cc$FtolRel)
+ nM$setMinfMax(cc$MinfMax)
+ nM$setIprint(cc$iprint)
+ while ((nMres <- nM$newf(devfun(nM$xeval()))) == 0L) {}
+ if (nMres < 0L) {
+ if (nMres > -4L)
+ stop("convergence failure, code ", nMres, " in NelderMead")
+ else
+ warning("failure to converge in 1000 evaluations")
+ }
+ list(fval=nM$value(), pars=nM$xpos(), code=nMres)
+ })
+ }
+ sqrLenU <- rho$pp$sqrL(0.)
+ wrss <- rho$resp$wrss()
+ pwrss <- wrss + sqrLenU
+ n <- nrow(fr)
+ dims <- c(N=n, n=n, nmp=n-p, nth=length(rho$pp$theta), p=p, q=nrow(reTrms$Zt),
+ nAGQ=nAGQ, useSc=1L, reTrms=length(reTrms$cnms),
+ spFe=0L, REML=0L, GLMM=0L, NLMM=1L)
+ cmp <- c(ldL2=rho$pp$ldL2(), ldRX2=rho$pp$ldRX2(), wrss=wrss,
+ ussq=sqrLenU, pwrss=pwrss, drsum=NA,
+ drsum=wrss, dev=opt$fval, REML=NA,
+ sigmaML=sqrt(pwrss/n), sigmaREML=NA)
+ new("nlmerMod", call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms,
+ Gp=reTrms$Gp, theta=rho$pp$theta, beta=rho$pp$beta0, u=rho$pp$u0,
+ lower=reTrms$lower, devcomp=list(cmp=cmp, dims=dims), pp=rho$pp, resp=rho$resp)
+}## {nlmer}
+
+
## Methods for the merMod class
fixef.merMod <- function(object, ...)
structure(object at beta, names = dimnames(object at pp$X)[[2]])
@@ -1076,20 +1135,20 @@
refitML.merMod <- function (x, ...) {
if (!isREML(x)) return(x)
stopifnot(is(rr <- x at resp, "lmerResp"))
- resp <- new(class(rr), y=rr$y, offset=rr$offset, weights=rr$weights, REML=0L)
-# resp$offset <- rr$offset
-# resp$weights <- rr$weights
-# resp$REML <- 0L
+ rho <- new.env(parent=parent.env(environment()))
+ rho$resp <- new(class(rr), y=rr$y, offset=rr$offset, weights=rr$weights, REML=0L)
xpp <- x at pp
- pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat,
- Lind=xpp$Lind, theta=xpp$theta)
- opt <- bobyqa(x at theta, mkdevfun(pp, resp, emptyenv()), x at lower)
+ rho$pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat,
+ Lind=xpp$Lind, theta=xpp$theta, S=1L)
+ devfun <- mkdevfun(rho, 0L)
+ opt <- bobyqa(x at theta, devfun, x at lower)
n <- length(rr$y)
+ pp <- rho$pp
p <- ncol(pp$X)
dims <- c(N=n, n=n, nmp=n-p, nth=length(pp$theta), p=p, q=nrow(pp$Zt),
nAGQ=NA_integer_, useSc=1L, reTrms=length(x at cnms),
spFe=0L, REML=0L, GLMM=0L, NLMM=0L)
- wrss <- resp$wrss()
+ wrss <- rho$resp$wrss()
ussq <- pp$sqrL(1)
pwrss <- wrss + ussq
cmp <- c(ldL2=pp$ldL2(), ldRX2=pp$ldRX2(), wrss=wrss, ussq=ussq,
@@ -1100,7 +1159,7 @@
### tricky to do so without causing the call to be evaluated
new("lmerMod", call=x at call, frame=x at frame, flist=x at flist,
cnms=x at cnms, theta=pp$theta, beta=pp$delb, u=pp$delu,
- lower=x at lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=resp)
+ lower=x at lower, devcomp=list(cmp=cmp, dims=dims), pp=pp, resp=rho$resp)
}
getCall.merMod <- function(x, ...) x at call
Modified: pkg/lme4Eigen/R/profile.R
===================================================================
--- pkg/lme4Eigen/R/profile.R 2012-01-07 17:29:59 UTC (rev 1497)
+++ pkg/lme4Eigen/R/profile.R 2012-01-07 23:35:27 UTC (rev 1498)
@@ -179,7 +179,9 @@
rr$weights <- fitted at resp$weights
fe.zeta <- function(fw) {
rr$offset <- Xw * fw + offset.orig
- ores <- bobyqa(thopt, mkdevfun(pp1, rr, emptyenv()), lower = fitted at lower)
+ rho <- as.environment(list(pp=pp1, resp=rr))
+ parent.env(rho) <- parent.frame()
+ ores <- bobyqa(thopt, mkdevfun(rho, 0L), lower = fitted at lower)
fv <- ores$fval
sig <- sqrt((rr$wrss() + pp1$sqrL(1))/n)
c(sign(fw - est) * sqrt(fv - base),
Modified: pkg/lme4Eigen/tests/lmer-1.R
===================================================================
--- pkg/lme4Eigen/tests/lmer-1.R 2012-01-07 17:29:59 UTC (rev 1497)
+++ pkg/lme4Eigen/tests/lmer-1.R 2012-01-07 23:35:27 UTC (rev 1498)
@@ -82,12 +82,12 @@
FUN <- get(nm)
F.fmX1s <- FUN(fmX1s)
# F.fmX2s <- FUN(fmX2s)
- if(nm == "model.matrix") {
- F.fmX1s <- as(F.fmX1s, "denseMatrix")
+# if(nm == "model.matrix") {
+# F.fmX1s <- as(F.fmX1s, "denseMatrix")
# F.fmX2s <- as(F.fmX2s, "denseMatrix")
# FF <- function(.) {r <- FUN(.); row.names(r) <- NULL
# as(r, "generalMatrix") }
- } # else
+# } # else
FF <- FUN
stopifnot(
all.equal( FF(fmX1), F.fmX1s, tol = 1e-6)
@@ -122,18 +122,17 @@
# family = binomial, data = cbpp, doFit = FALSE)
## now
#bobyqa(m1e, control = list(iprint = 2L))
-m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
- family = binomial, data = cbpp, verbose = 2L)
+m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), cbpp, binomial, nAGQ=25L)
stopifnot(is((cm1 <- coef(m1)), "coef.mer"),
dim(cm1$herd) == c(15,4),
- all.equal(fixef(m1), ## these values are those of "old-lme4":
- c(-1.39853504914, -0.992334711,
- -1.12867541477, -1.58037390498),
- tol = 1.e-3,
- check.attr=FALSE)
+ all.equal(fixef(m1), ## these values are from an Ubuntu 11.10 amd64 system
+ c(-1.39922135307046, -0.991415396352428,
+ -1.12781521322006, -1.57947198508598),
+ tol = 1.e-7,
+ check.attr=FALSE),
+ all.equal(deviance(m1), 100.010030539916, tol=1e-7)
)
-## Deviance for the new algorithm is lower, eventually we should change the previous test
-#stopifnot(deviance(m1) <= deviance(m1e))
+
## Simple example by Andrew Gelman (2006-01-10) ----
n.groups <- 10 ; n.reps <- 2
n <- length(group.id <- gl(n.groups, n.reps))
@@ -151,7 +150,7 @@
all.equal(ranef(fit.1)[["group.id"]][,"(Intercept)"],
c(1.80469, -1.80977, 1.61465, 1.54083, -0.1332,
-3.33067, -1.82593, -0.873515, -0.359131, 3.37204),
- tol = 1e-4)
+ tol = 1e-5)
)
@@ -176,7 +175,7 @@
contrasts(bacteria$trt) <-
structure(contr.sdif(3),
dimnames = list(NULL, c("diag", "encourage")))
- print(fm5 <- glmer(y ~ trt + wk2 + (1|ID), bacteria, binomial))
+ print(fm5 <- glmer(y ~ trt + wk2 + (1|ID), bacteria, binomial, nAGQ=25L))
## used to fail with nlminb() : stuck at theta=1
showProc.time() #
@@ -184,16 +183,15 @@
stopifnot(
all.equal(logLik(fm5),
## was -96.127838
- structure(-96.13069, nobs = 220L, nall = 220L,
+ structure(-95.89706, nobs = 220L, nall = 220L,
df = 5L, REML = FALSE,
class = "logLik"),
tol = 1e-5, check.attributes = FALSE)
,
all.equal(fixef(fm5),
- ## was 2.834218798 -1.367099481
- c("(Intercept)"= 2.831609490, "trtdiag"= -1.366722631,
- ## now 0.5842291915, -1.599148773
- "trtencourage"=0.5840147802, "wk2TRUE"=-1.598591346), tol = 1e-4)
+ c("(Intercept)"= 2.85970407988865, "trtdiag"= -1.36896064623335,
+ "trtencourage"=0.579864265134738, "wk2TRUE"=-1.62687300090901),
+ tol = 1e-6)
)
}
More information about the Lme4-commits
mailing list