[Lme4-commits] r1679 - in pkg/lme4: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 19 19:24:07 CET 2012
Author: bbolker
Date: 2012-03-19 19:24:07 +0100 (Mon, 19 Mar 2012)
New Revision: 1679
Added:
pkg/lme4/tests/is.R
Modified:
pkg/lme4/NAMESPACE
pkg/lme4/R/AllGeneric.R
pkg/lme4/R/lmer.R
pkg/lme4/R/profile.R
pkg/lme4/man/isREML.Rd
pkg/lme4/man/profile-methods.Rd
pkg/lme4/tests/glmmExt.R
pkg/lme4/tests/profile.R
Log:
exposed is[GN]*LMM, added documentation to isREML.Rd, added tests
added 'which' to profiling
Modified: pkg/lme4/NAMESPACE
===================================================================
--- pkg/lme4/NAMESPACE 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/NAMESPACE 2012-03-19 18:24:07 UTC (rev 1679)
@@ -18,6 +18,9 @@
export(GQdk)
export(isNested)
export(isREML)
+export(isLMM)
+export(isGLMM)
+export(isNLMM)
export(lmer)
export(lmerResp)
export(lmList)
@@ -91,7 +94,6 @@
importMethodsFrom(Matrix,t)
importMethodsFrom(Matrix,tcrossprod)
S3method(anova,merMod)
-S3method(as.data.frame,thpr)
S3method(as.function,merMod)
S3method(coef,lmList)
S3method(coef,merMod)
@@ -109,10 +111,12 @@
S3method(fitted,merMod)
S3method(fixef,merMod)
S3method(formula,lmList)
+S3method(isLMM,merMod)
+S3method(isGLMM,merMod)
+S3method(isNLMM,merMod)
S3method(isREML,merMod)
S3method(logLik,merMod)
S3method(log,thpr)
-S3method(mcmcsamp,merMod)
S3method(model.frame,merMod)
S3method(model.matrix,merMod)
S3method(nobs,merMod)
Modified: pkg/lme4/R/AllGeneric.R
===================================================================
--- pkg/lme4/R/AllGeneric.R 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/R/AllGeneric.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -21,7 +21,7 @@
## Create a Markov chain Monte Carlo sample from the posterior
## distribution of the parameters
##
-##
+##
## @title Create an MCMC sample
## @param object a fitted model object
## @param n number of samples to generate. Defaults to 1; for real use values of 200-1000 are more typical
@@ -41,21 +41,37 @@
##' @export
sigma <- function(object, ...) UseMethod("sigma")
-##' Check if a model has been fit according to the REML criterion
+##' Check characteristics of models: whether a model fit corresponds to a linear (LMM), generalized linear (GLMM), or
+##' nonlinear (NLMM) mixed model, and whether a linear mixed model has been fitted by REML or not (\code{isREML(x)}
+##' is always \code{FALSE} for GLMMs and NLMMs).
##'
##' This is a generic function. At present the only methods are for mixed-effects
##' models of class \code{\linkS4class{merMod}}.
-##' @title Check for a REML fit
+##' @title Check characteristics of models
##' @param x a fitted model.
-##' @param ... additional, optional arguments. (None are used in the merMod method)
-##' @return \code{TRUE} if \code{x} has been fit by REML, otherwise FALSE
+##' @param ... additional, optional arguments. (None are used in the merMod methods)
+##' @return a logical value
+##' @seealso getME
##' @export
isREML <- function(x, ...) UseMethod("isREML")
+##' @rdname isREML
+##' @export
+isLMM <- function(x, ...) UseMethod("isLMM")
+
+##' @rdname isREML
+##' @export
+isNLMM <- function(x, ...) UseMethod("isNLMM")
+
+##' @rdname isREML
+##' @export
+isGLMM <- function(x, ...) UseMethod("isGLMM")
+
##' Refit a model using the maximum likelihood criterion
##'
##' This function is primarily used to get a maximum likelihood fit of
##' a linear mixed-effects model for an \code{\link{anova}} comparison.
+##'
##' @title Refit a model by maximum likelihood criterion
##' @param x a fitted model, usually of class \code{"\linkS4class{lmerMod}"},
##' to be refit according to the maximum likelihood criterion
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/R/lmer.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -872,6 +872,24 @@
##' @S3method isREML merMod
isREML.merMod <- function(x, ...) as.logical(x at devcomp$dims["REML"])
+##' @S3method isGLMM merMod
+isGLMM.merMod <- function(x,...) {
+ as.logical(x at devcomp$dims["GLMM"])
+ ## or: is(x at resp,"glmResp")
+}
+
+##' @S3method isNLMM merMod
+isNLMM.merMod <- function(x,...) {
+ as.logical(x at devcomp$dims["NLMM"])
+ ## or: is(x at resp,"nlsResp")
+}
+
+##' @S3method isLMM merMod
+isLMM.merMod <- function(x,...) {
+ !isGLMM(x) && !isNLMM(x)
+ ## or: is(x at resp,"lmerResp") ?
+}
+
##' @importFrom stats logLik
##' @S3method logLik merMod
logLik.merMod <- function(object, REML = NULL, ...) {
@@ -1312,6 +1330,11 @@
}
if (!is.null(cc <- so$call$formula))
cat("Formula:", deparse(cc),"\n")
+ if (!is.null(so$family)) {
+ cat("Family: ",so$family,
+ " (link=",so$link,")\n",
+ sep="")
+ }
if (!is.null(cc <- so$call$data))
cat(" Data:", deparse(cc), "\n")
if (!is.null(cc <- so$call$subset))
@@ -1711,8 +1734,8 @@
link <- fam <- NULL
if(is(resp, "glmerResp")) {
- fam <- resp$fam()
- link <- resp$link()
+ fam <- resp$family$family
+ link <- resp$family$link
}
coefs <- cbind("Estimate" = fixef(object),
"Std. Error" = sig * sqrt(diag(object at pp$unsc())))
@@ -1887,21 +1910,7 @@
object at resp$weights
}
-isGLMM <- function(object) {
- as.logical(object at devcomp$dims["GLMM"])
- ## or: is(object at resp,"glmResp")
-}
-isNLMM <- function(object) {
- as.logical(object at devcomp$dims["NLMM"])
- ## or: is(object at resp,"nlsResp")
-}
-
-isLMM <- function(object) {
- !isGLMM(object) && !isNLMM(object)
- ## **not** is(object at resp,"lmerResp") ?
-}
-
getOptfun <- function(optimizer) {
if (is.character(optimizer)) {
optfun <- try(get(optimizer),silent=TRUE)
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/R/profile.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -8,7 +8,7 @@
##' @aliases profile-methods profile.merMod
##' @docType methods
##' @param fitted a fitted model, e.g., the result of \code{\link{lmer}(..)}.
-##' @param which which parameters to profile (UNDER CONSTRUCTION): default is all parameters
+##' @param which (integer) which parameters to profile: default is all parameters. The parameters are ordered as follows: (1) random effects (theta) parameters; (2) residual standard deviation (or scale parameter for GLMMs where appropriate); (3) fixed effect parameters. Random effects parameters are ordered as in \code{getME(.,"theta")}, i.e. as the lower triangle of a matrix with standard deviations on the diagonal and correlations off the diagonal.
##' @param alphamax maximum alpha value for likelihood ratio confidence regions; used to establish the range of values to be profiled
##' @param maxpts maximum number of points (in each direction, for each parameter) to evaluate in attempting to construct the profile
##' @param delta stepping scale for deciding on next point to profile
@@ -18,7 +18,7 @@
##' @param optimizer (character or function) optimizer to use (see \code{\link{lmer}} for details)
##' @param \dots potential further arguments for \code{profile} methods.
##' @section Methods: FIXME: These (signatures) will change soon --- document
-##' \bold{after} change!
+##' \bold{after} change!
##' \describe{
##' \item{signature(fitted = \"merMod\")}{ ... } }
##' @seealso For (more expensive) alternative confidence intervals:
@@ -30,6 +30,15 @@
##' system.time( tpr <- profile(fm01ML) )
##' (confint(tpr) -> CIpr)
##' print(xyplot(tpr))
+##' tpr2 <- profile(fm01ML, which=1:2) ## Batch and residual variance only
+##' ## GLMM example (running time ~8 seconds on a modern machine)
+##' \dontrun{
+##' gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+##' data = cbpp, family = binomial)
+##' system.time(pr4 <- profile(gm1)) ## ~ 7 seconds
+##' xyplot(pr4,layout=c(5,1),as.table=TRUE)
+##' splom(pr4)
+##' }
##' @importFrom splines backSpline interpSpline periodicSpline
##' @importFrom stats profile
##' @method profile merMod
@@ -42,7 +51,7 @@
## 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 selection of individual variables to profile by name?
## FIXME: allow for failure of bounds (non-pos-definite correlation matrices) when >1 cor parameter
## FIXME: generalize to GLMMs
## (use different devfun;
@@ -63,10 +72,13 @@
n <- environment(dd)$n
p <- length(pp$beta0)
- ans <- lapply(opt <- attr(dd, "optimum"), function(el) NULL)
+ opt <- attr(dd, "optimum")
+ nptot <- length(opt)
+
+ ans <- lapply(opt[which], function(el) NULL)
bakspl <- forspl <- ans
- nptot <- length(opt)
+
nvp <- nptot - p # number of variance-covariance pars
fe.orig <- opt[-seq_len(nvp)]
res <- c(.zeta = 0, opt)
@@ -109,8 +121,8 @@
nextstart <- function(mat, pind, r, method="global") {
## FIXME: indexing may need to be checked (e.g. for fixed-effect parameters)
switch(method,
- global=opt[seqnvp][-pind], ## address opt, no zeta column
- prev=mat[r,1+seqnvp][-pind],
+ global=opt[seqpar1][-pind], ## address opt, no zeta column
+ prev=mat[r,1+seqpar1][-pind],
extrap=stop("stub")) ## do something with mat[r-(1:0),1+seqnvp])[-pind] ...
}
@@ -166,14 +178,15 @@
stopifnot(all.equal(unname(dd(opt[seq(npar1)])),base,tol=1e-5))
seqnvp <- intersect(seq_len(npar1),which)
- lowvp <- lower[seqnvp]
- upvp <- upper[seqnvp]
+ seqpar1 <- seq_len(npar1)
+ lowvp <- lower[seqpar1]
+ upvp <- upper[seqpar1]
form <- .zeta ~ foo # pattern for interpSpline formula
for (w in seqnvp) {
if (verbose) cat(if (isLMM(fitted)) "var-cov " else "", "parameter ",w,":\n",sep="")
wp1 <- w + 1L
- start <- opt[seqnvp][-w]
+ start <- opt[seqpar1][-w]
pw <- opt[w]
lowcut <- lower[w]
upcut <- upper[w]
@@ -208,18 +221,19 @@
pres <- nres <- res
## assign one row, determined by inc. sign, from a small shift
## FIXME:: do something if pw==0 ???
- nres[1, ] <- pres[2, ] <- zeta(pw * 1.01, start=opt[seqnvp][-w])
+ nres[1, ] <- pres[2, ] <- zeta(pw * 1.01, start=opt[seqpar1][-w])
## fill in the rest of the arrays and collapse them
bres <-
as.data.frame(unique(rbind2(fillmat(pres,lowcut, upcut, zeta, wp1),
fillmat(nres,lowcut, upcut, zeta, wp1))))
- bres$.par <- names(opt)[w]
- ans[[w]] <- bres[order(bres[, wp1]), ]
- form[[3]] <- as.name(names(opt)[w])
+ pname <- names(opt)[w]
+ bres$.par <- pname
+ ans[[pname]] <- bres[order(bres[, wp1]), ]
+ form[[3]] <- as.name(pname)
## FIXME: test for bad things here??
- bakspl[[w]] <- try(backSpline(forspl[[w]] <- interpSpline(form, bres)),silent=TRUE)
- if (inherits(bakspl[[w]],"try-error")) {
+ bakspl[[pname]] <- try(backSpline(forspl[[pname]] <- interpSpline(form, bres)),silent=TRUE)
+ if (inherits(bakspl[[pname]],"try-error")) {
warning("non-monotonic profile")
}
@@ -229,6 +243,7 @@
if (isLMM(fitted)) {
offset.orig <- fitted at resp$offset
fp <- seq_len(p)
+ fp <- fp[(fp+nvp) %in% which]
for (j in fp) {
if (verbose) cat("fixed-effect parameter ",j,":\n",sep="")
pres <- # intermediate results for pos. incr.
Modified: pkg/lme4/man/isREML.Rd
===================================================================
--- pkg/lme4/man/isREML.Rd 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/man/isREML.Rd 2012-03-19 18:24:07 UTC (rev 1679)
@@ -1,26 +1,41 @@
\name{isREML}
+\alias{isGLMM}
+\alias{isLMM}
+\alias{isNLMM}
\alias{isREML}
-\title{Check for a REML fit}
+\title{Check characteristics of models}
\usage{
isREML(x, ...)
+
+ isLMM(x, ...)
+
+ isNLMM(x, ...)
+
+ isGLMM(x, ...)
}
\arguments{
\item{x}{a fitted model.}
\item{...}{additional, optional arguments. (None are
- used in the merMod method)}
+ used in the merMod methods)}
}
\value{
- \code{TRUE} if \code{x} has been fit by REML, otherwise
- FALSE
+ a logical value
}
\description{
- Check if a model has been fit according to the REML
- criterion
+ Check characteristics of models: whether a model fit
+ corresponds to a linear (LMM), generalized linear (GLMM),
+ or nonlinear (NLMM) mixed model, and whether a linear
+ mixed model has been fitted by REML or not
+ (\code{isREML(x)} is always \code{FALSE} for GLMMs and
+ NLMMs).
}
\details{
This is a generic function. At present the only methods
are for mixed-effects models of class
\code{\linkS4class{merMod}}.
}
+\seealso{
+ getME
+}
Modified: pkg/lme4/man/profile-methods.Rd
===================================================================
--- pkg/lme4/man/profile-methods.Rd 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/man/profile-methods.Rd 2012-03-19 18:24:07 UTC (rev 1679)
@@ -4,15 +4,25 @@
\alias{profile-methods}
\title{Methods for profile() of [ng]lmer fitted models}
\usage{
- \method{profile}{merMod} (fitted, alphamax = 0.01,
- maxpts = 100, delta = cutoff/8, verbose = 0,
- devtol = 1e-09, startval = "prev",
+ \method{profile}{merMod} (fitted, which = 1:nptot,
+ alphamax = 0.01, maxpts = 100, delta = cutoff/8,
+ verbose = 0, devtol = 1e-09, startmethod = "prev",
optimizer = "bobyqa", ...)
}
\arguments{
\item{fitted}{a fitted model, e.g., the result of
\code{\link{lmer}(..)}.}
+ \item{which}{(integer) which parameters to profile:
+ default is all parameters. The parameters are ordered as
+ follows: (1) random effects (theta) parameters; (2)
+ residual standard deviation (or scale parameter for GLMMs
+ where appropriate); (3) fixed effect parameters. Random
+ effects parameters are ordered as in
+ \code{getME(.,"theta")}, i.e. as the lower triangle of a
+ matrix with standard deviations on the diagonal and
+ correlations off the diagonal.}
+
\item{alphamax}{maximum alpha value for likelihood ratio
confidence regions; used to establish the range of values
to be profiled}
@@ -30,7 +40,7 @@
\item{devtol}{tolerance for fitted deviances less than
baseline (supposedly minimum) deviance}
- \item{startval}{method for picking starting conditions
+ \item{startmethod}{method for picking starting conditions
for optimization (STUB)}
\item{optimizer}{(character or function) optimizer to use
@@ -55,7 +65,16 @@
system.time( tpr <- profile(fm01ML) )
(confint(tpr) -> CIpr)
print(xyplot(tpr))
+tpr2 <- profile(fm01ML, which=1:2) ## Batch and residual variance only
+## GLMM example (running time ~8 seconds on a modern machine)
+\dontrun{
+gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ data = cbpp, family = binomial)
+system.time(pr4 <- profile(gm1)) ## ~ 7 seconds
+xyplot(pr4,layout=c(5,1),as.table=TRUE)
+splom(pr4)
}
+}
\seealso{
For (more expensive) alternative confidence intervals:
\code{\link{bootMer}}.
Modified: pkg/lme4/tests/glmmExt.R
===================================================================
--- pkg/lme4/tests/glmmExt.R 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/tests/glmmExt.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -61,9 +61,12 @@
gG1 <- glmer(y ~ 1 + (1|block), data=dG, family=gaussian(link="log"), verbose=TRUE)
gG2 <- glmer(y ~ x + (1|block), data=dG, family=gaussian(link="log"), verbose=TRUE)
+
+## 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))
## sets variance to zero, converges back to GLM solution
+## FIXME: cloglog link
Added: pkg/lme4/tests/is.R
===================================================================
--- pkg/lme4/tests/is.R (rev 0)
+++ pkg/lme4/tests/is.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -0,0 +1,28 @@
+library(lme4)
+fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+stopifnot(isREML(fm1),
+ isLMM(fm1),
+ !isGLMM(fm1),
+ !isNLMM(fm1))
+
+fm1ML <- refitML(fm1)
+stopifnot(!isREML(fm1ML),
+ isLMM(fm1ML),
+ !isGLMM(fm1ML),
+ !isNLMM(fm1ML))
+
+gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ data = cbpp, family = binomial)
+stopifnot(!isREML(gm1),
+ !isLMM(gm1),
+ isGLMM(gm1),
+ !isNLMM(gm1))
+
+nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
+ Orange, start = c(Asym = 200, xmid = 725, scal = 350))
+stopifnot(!isREML(nm1),
+ !isLMM(nm1),
+ !isGLMM(nm1),
+ isNLMM(nm1))
+
+
Modified: pkg/lme4/tests/profile.R
===================================================================
--- pkg/lme4/tests/profile.R 2012-03-18 19:07:33 UTC (rev 1678)
+++ pkg/lme4/tests/profile.R 2012-03-19 18:24:07 UTC (rev 1679)
@@ -6,6 +6,15 @@
##
system.time( tpr <- profile(fm01ML) )
+## test all combinations of 'which'
+wlist <- list(1:3,1:2,1,2:3,2,3,c(1,3))
+invisible(lapply(wlist,function(w) xyplot(profile(fm01ML,which=w))))
+tpr2 <- profile(fm01ML,which=1:3)
+tpr3 <- profile(fm01ML,which=2:3) ## can't plot
+tpr4 <- profile(fm01ML,which=3)
+tpr5 <- profile(fm01ML,which=2) ## can't plot
+tpr6 <- profile(fm01ML,which=1)
+
(confint(tpr) -> CIpr)
print(xyplot(tpr))
## comparing against lme4a reference values -- but lme4 returns sigma
More information about the Lme4-commits
mailing list