[Lme4-commits] r1683 - in pkg/lme4: . R inst/tests tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 20 04:50:54 CET 2012
Author: bbolker
Date: 2012-03-20 04:50:54 +0100 (Tue, 20 Mar 2012)
New Revision: 1683
Added:
pkg/lme4/tests/REMLdev.R
Modified:
pkg/lme4/NAMESPACE
pkg/lme4/R/lmer.R
pkg/lme4/inst/tests/test-lmer.R
pkg/lme4/tests/drop.R
Log:
too much at once
updated terms: use nobars() properly rather than trying to hack
it by removing items from @flist.
- Fixes problems with empty models in drop1, anova, and
problem with SASmixed ?Cultivation example
- fix bug in drop1() for empty scope (seq(length()) vs seq_along())
- add empty models to tests in drop.R
started on (deviance|logLik)(x,REML) implementation
- commented is.na(deviance()) test
- (deviance|loglik)(REMLfit,REML=FALSE) not yet working
- fix (major) typo in REMLdev
- added (incomplete) REMLdev.R tests
added formula.merMod (PROBABLY needs fix for eval() issues)
Modified: pkg/lme4/NAMESPACE
===================================================================
--- pkg/lme4/NAMESPACE 2012-03-20 03:50:20 UTC (rev 1682)
+++ pkg/lme4/NAMESPACE 2012-03-20 03:50:54 UTC (rev 1683)
@@ -111,6 +111,7 @@
S3method(fitted,merMod)
S3method(fixef,merMod)
S3method(formula,lmList)
+S3method(formula,merMod)
S3method(isLMM,merMod)
S3method(isGLMM,merMod)
S3method(isNLMM,merMod)
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-03-20 03:50:20 UTC (rev 1682)
+++ pkg/lme4/R/lmer.R 2012-03-20 03:50:54 UTC (rev 1683)
@@ -653,7 +653,7 @@
val <- data.frame(Df = Df,
AIC = .sapply(llks, AIC),
BIC = .sapply(llks, BIC),
- logLik = llk,
+ logLik = llk,
deviance = -2*llk,
Chisq = chisq,
"Chi Df" = dfChisq,
@@ -678,7 +678,8 @@
names(ss) <- colnames(X)
terms <- terms(object)
## FIXME: this setdiff() should be obsolete since terms now keeps only fixed effects by default
- nmeffects <- setdiff(attr(terms, "term.labels"), names(object at flist))
+ ## nmeffects <- setdiff(attr(terms, "term.labels"), names(object at flist))
+ nmeffects <- attr(terms, "term.labels")
if ("(Intercept)" %in% names(ss))
nmeffects <- c("(Intercept)", nmeffects)
ss <- unlist(lapply(split(ss, asgn), sum))
@@ -742,11 +743,13 @@
##' @S3method coef merMod
coef.merMod <- coefMer
+## FIXME: should this be computed and stored?
REMLdev <- function(object) {
n <- object at devcomp$dims["N"] ## FIXME: N or n ??
+ nmp <- n-length(object at beta)
cmp <- object at devcomp$cmp
- cmp["ldL2"]+cmp["ldRX2"]+n*(1+log(2*pi*cmp["pwrss"]/n))
- ## from lme4:
+ cmp["ldL2"]+cmp["ldRX2"]+nmp*(1+log(2*pi*cmp["pwrss"]/nmp))
+ ## from lme4 (old):
## d[REML_POS] = d[ldL2_POS] + d[ldRX2_POS] +
## dnmp * (1. + log(d[pwrss_POS]) + log(2. * PI / dnmp));
}
@@ -754,9 +757,16 @@
##' @importFrom stats deviance
##' @S3method deviance merMod
deviance.merMod <- function(object, REML = NULL, ...) {
- if (!missing(REML)) stop("REML argument not supported")
- ##
- object at devcomp$cmp[["dev"]]
+ if (isTRUE(REML) && !isLMM(object)) stop("can't compute REML deviance for a non-LMM")
+ if (is.null(REML) || is.na(REML[1]))
+ REML <- isREML(object)
+ if (REML) REMLdev(object) else {
+ if (!isREML(object)) {
+ object at devcomp$cmp[["dev"]]
+ } else {
+ ## need to compute ML deviance from scratch here ...
+ }
+ }
}
##' @importFrom stats drop1
@@ -778,7 +788,7 @@
ans[1, ] <- extractAIC(object, scale, k = k, ...)
n0 <- nobs(object, use.fallback = TRUE)
env <- environment(formula(object))
- for(i in seq(ns)) {
+ for(i in seq_along(scope)) { ## was seq(ns), failed on empty scope
tt <- scope[i]
if(trace > 1) {
cat("trying -", tt, "\n", sep='')
@@ -893,11 +903,12 @@
##' @importFrom stats logLik
##' @S3method logLik merMod
logLik.merMod <- function(object, REML = NULL, ...) {
- if (!missing(REML)) stop("REML argument not supported")
+ if (is.null(REML) || is.na(REML[1]))
+ REML <- isREML(object)
+ val <- -deviance(object, REML = REML)/2
dc <- object at devcomp
dims <- dc$dims
- val <- - dc$cmp[["dev"]]/2
- attr(val, "nall") <- attr(val, "nobs") <- nrow(object at frame)
+ attr(val, "nall") <- attr(val, "nobs") <- nrow(object at frame) ## FIXME use nobs() ?
attr(val, "df") <- length(object at beta) + length(object at theta) + dims[["useSc"]]
class(val) <- "logLik"
val
@@ -1290,7 +1301,11 @@
terms.merMod <- function(x, fixed.only=TRUE, ...) {
tt <- attr(x at frame, "terms")
if (fixed.only) {
- drop.terms(tt,match(names(x at flist),attr(tt,"term.labels")))
+ ## FIXME: e.g. ~ Inoc * Cult + (1|Block) + (1|Cult)
+ ## need to use nobars() to get the right answer
+ ## drop.terms(tt,match(names(x at flist),attr(tt,"term.labels")))
+ mm <- model.frame(as.formula(nobars(formula(x)[-2])),data=model.frame(x))
+ terms(mm)
} else tt
}
@@ -1987,3 +2002,7 @@
opt
}
+formula.merMod <- function(x,...) {
+ x at call$formula
+}
+
Modified: pkg/lme4/inst/tests/test-lmer.R
===================================================================
--- pkg/lme4/inst/tests/test-lmer.R 2012-03-20 03:50:20 UTC (rev 1682)
+++ pkg/lme4/inst/tests/test-lmer.R 2012-03-20 03:50:54 UTC (rev 1683)
@@ -13,7 +13,7 @@
expect_that(REMLfun(0), equals(326.023232155879))
expect_that(family(fm1), equals(gaussian()))
expect_that(isREML(fm1ML <- refitML(fm1)), equals(FALSE))
- expect_that(is.na(deviance(fm1)), equals(TRUE))
+ ## expect_that(is.na(deviance(fm1)), equals(TRUE))
expect_that(deviance(fm1ML), equals(327.327059881135))
expect_that(sigma(fm1), equals(49.5100503990048))
expect_that(sigma(fm1ML), equals(49.5100999308089))
Added: pkg/lme4/tests/REMLdev.R
===================================================================
--- pkg/lme4/tests/REMLdev.R (rev 0)
+++ pkg/lme4/tests/REMLdev.R 2012-03-20 03:50:54 UTC (rev 1683)
@@ -0,0 +1,25 @@
+library(lme4)
+fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+fm1ML <- refitML(fm1)
+deviance(fm1)
+deviance(fm1ML)
+deviance(fm1,REML=FALSE) ## FIXME: not working yet (NA)
+deviance(fm1,REML=TRUE)
+
+## from lme4.0
+oldvals <- c(REML=1743.6282722424, ML=1751.98581103058)
+## leave out ML values for REML fits for now ...
+stopifnot(
+ all.equal(deviance(fm1),deviance(fm1,REML=TRUE),deviance(fm1ML,REML=TRUE),oldvals["REML"]),
+ all.equal(deviance(fm1ML),deviance(fm1ML,REML=FALSE),oldvals["ML"]),
+ all.equal(deviance(fm1)/-2,c(logLik(fm1)),c(logLik(fm1ML,REML=TRUE)),c(logLik(fm1,REML=TRUE))),
+ all.equal(deviance(fm1ML)/-2,c(logLik(fm1ML,REML=FALSE)),
+ c(logLik(fm1ML,REML=FALSE))))
+
+## should be:
+## stopifnot(
+## all.equal(deviance(fm1),deviance(fm1,REML=TRUE),deviance(fm1ML,REML=TRUE),oldvals["REML"]),
+## all.equal(deviance(fm1ML),deviance(fm1,REML=FALSE),deviance(fm1ML,REML=FALSE),oldvals["ML"]),
+## all.equal(deviance(fm1)/2,c(logLik(fm1)),c(logLik(fm1ML,REML=TRUE)),c(logLik(fm1,REML=TRUE))),
+## all.equal(deviance(fm1ML)/2,c(logLik(fm1,REML=FALSE)),c(logLik(fm1ML,REML=FALSE)),
+## c(logLik(fm1ML,REML=FALSE))))
Modified: pkg/lme4/tests/drop.R
===================================================================
--- pkg/lme4/tests/drop.R 2012-03-20 03:50:20 UTC (rev 1682)
+++ pkg/lme4/tests/drop.R 2012-03-20 03:50:54 UTC (rev 1683)
@@ -4,8 +4,9 @@
## slightly weird model but plausible --- not that
## one would want to try drop1() on this model ...
fm2 <- lmer(Reaction ~ 1+ (Days|Subject), sleepstudy)
+drop1(fm2) ## empty
update(fm1, . ~ . - Days)
-try(anova(fm2)) ## fails because there is nothing to test
+anova(fm2) ## empty
terms(fm1)
terms(fm1,fixed.only=FALSE)
More information about the Lme4-commits
mailing list