[Lme4-commits] r1685 - pkg/lme4/R

[noreply at r-forge.r-project.org]

Author: bbolker
Date: 2012-03-20 14:18:30 +0100 (Tue, 20 Mar 2012)
New Revision: 1685

Modified:
   pkg/lme4/R/lmer.R
Log:

   rearrange ML/REML deviance stuff
   clean up terms.merMod
   various tweaks to comments and FIXMEs



Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-03-20 04:20:36 UTC (rev 1684)
+++ pkg/lme4/R/lmer.R	2012-03-20 13:18:30 UTC (rev 1685)
@@ -677,8 +677,6 @@
 	ss <- as.vector(object at pp$RX() %*% object at beta)^2
 	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 <- attr(terms, "term.labels")
 	if ("(Intercept)" %in% names(ss))
 	    nmeffects <- c("(Intercept)", nmeffects)
@@ -743,28 +741,34 @@
 ##' @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"]+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));
-}
-
+## FIXME: should these values (i.e. ML criterion for REML models
+##  and vice versa) be computed and stored in the object in the first place?
 ##' @importFrom stats deviance
 ##' @S3method deviance merMod
 deviance.merMod <- function(object, REML = NULL, ...) {
-    if (isTRUE(REML) && !isLMM(object)) stop("can't compute REML deviance for a non-LMM")
+    ## cf. (1) lmerResp::Laplace in respModule.cpp
+    ##     (2) section 5.6 of lMMwR, listing lines 34-42
+    if (isTRUE(REML) && !isLMM(object))
+        stop("can't compute REML deviance for a non-LMM")
+    cmp <- object at devcomp$cmp
     if (is.null(REML) || is.na(REML[1]))
         REML <- isREML(object)
-    if (REML) REMLdev(object) else {
+    if (REML) {
+        if (isREML(object)) {
+            cmp["REML"]
+        } else {
+            lnum <- log(2*pi*(cmp["pwrss"]))
+            n <- object at devcomp$dims["n"]
+            nmp <- n-length(object at beta)
+            unname(cmp["ldL2"]+cmp["ldRX2"]+nmp*(1.+lnum-log(nmp)))
+        }
+    } else {
         if (!isREML(object)) {
-            object at devcomp$cmp[["dev"]]
+            cmp[["dev"]]
         } else {
-            ## need to compute ML deviance from scratch here ...
+            n <- object at devcomp$dims["n"]
+            lnum <- log(2*pi*(cmp["pwrss"]))
+            unname(cmp["ldL2"]+n*(1+lnum-log(n)))
         }
     }
 }
@@ -1062,8 +1066,6 @@
         } else newresp <- newresp[-na.act]
       }
 
-
-
     if (isGLMM(object) && rr$family$family=="binomial") {
         if (is.matrix(newresp) && ncol(newresp)==2) {
             ntot <- rowSums(newresp)
@@ -1257,6 +1259,7 @@
 	      if(is.null(family$family)) stop("'family' not recognized")
 	      musim <- family$linkinv(etasim)
 	      ntot <- length(musim) ## FIXME: or could be dims["n"]?
+              ## FIXME: is it possible to leverage family$simulate ... ???
               val <- switch(family$family,
 			    poisson=rpois(ntot,lambda=musim),
 			    binomial={
@@ -1303,12 +1306,12 @@
 ##' @importFrom stats terms
 ##' @S3method terms merMod
 terms.merMod <- function(x, fixed.only=TRUE, ...) {
-  tt <- attr(x at frame, "terms")
   if (fixed.only) {
       ## do need to drop LHS (e.g. binomial example)
-      mm <- model.frame(formula(x,fixed.only=TRUE)[-2],data=model.frame(x))
-      terms(mm)
-  } else tt
+      ## mm <- model.frame(formula(x,fixed.only=TRUE)[-2],data=model.frame(x))
+      ## terms(mm)
+      terms.formula(formula(x,fixed.only=TRUE))
+  } else attr(x at frame,"terms")
 }
 
 ##' @importFrom stats update
@@ -1485,7 +1488,7 @@
 ##'     \item{devcomp}{a list consisting of a named numeric vector, \dQuote{cmp}, and
 ##'                    a named integer vector, \dQuote{dims}, describing the fitted model}
 ##'     \item{offset}{model offset}
-##'     \item{lower}{lower bounds on model parameters} ## FIXME: theta only?
+##'     \item{lower}{lower bounds on model parameters (random effects parameters only)}
 ##' }
 ##' @return Unspecified, as very much depending on the \code{\link{name}}.
 ##' @seealso \code{\link{getCall}()},
@@ -1566,7 +1569,8 @@
 
            "devcomp" = dc,
            "offset" = rsp$offset,
-           "lower" = object at lower,  ## FIXME: should this include -Inf values for beta?
+           "lower" = object at lower,
+            ## FIXME: current version gives lower bounds for theta parameters only -- these must be extended for [GN]LMMs -- give extended value including -Inf values for beta values?
 	   "..foo.." =# placeholder!
 	   stop(gettextf("'%s' is not implemented yet",
 			 sprintf("getME(*, \"%s\")", name))),
@@ -1676,7 +1680,7 @@
     if (is.null(cnms <- x at cnms))
 	stop("VarCorr methods require reTrms, not just reModule")
     if(missing(sigma)) # "bug": fails via default 'sigma=sigma(x)'
-	sigma <- lme4::sigma(x)  ## FIXME: do we need lme4:: ?
+	sigma <- lme4::sigma(x)  ## FIXME: do we still need lme4:: ?
     nc <- sapply(cnms, length)	  # no. of columns per term
     m <- mkVarCorr(sigma, cnms=cnms, nc=nc, theta = x at theta,
 	      nms = {fl <- x at flist; names(fl)[attr(fl, "assign")]})