[Lme4-commits] r1767 - in pkg/lme4: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 16 19:57:36 CEST 2012
Author: mmaechler
Date: 2012-06-16 19:57:36 +0200 (Sat, 16 Jun 2012)
New Revision: 1767
Added:
pkg/lme4/man/densityplot.thpr.Rd
pkg/lme4/man/splom.thpr.Rd
Modified:
pkg/lme4/NAMESPACE
pkg/lme4/R/AllGeneric.R
pkg/lme4/R/GHrule.R
pkg/lme4/R/lmer.R
pkg/lme4/R/mcmcsamp.R
pkg/lme4/R/plot.R
pkg/lme4/R/profile.R
pkg/lme4/man/VerbAgg.Rd
pkg/lme4/man/getME.Rd
pkg/lme4/man/isREML.Rd
pkg/lme4/man/mkdevfun.Rd
pkg/lme4/man/profile-methods.Rd
pkg/lme4/man/refit.Rd
pkg/lme4/man/refitML.Rd
pkg/lme4/tests/nlmer.Rout.save
Log:
isLMM() {etc] examples; doxy <-> Rd consistency; document densityplot() and splom() methods; other docs
Modified: pkg/lme4/NAMESPACE
===================================================================
--- pkg/lme4/NAMESPACE 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/NAMESPACE 2012-06-16 17:57:36 UTC (rev 1767)
@@ -99,6 +99,7 @@
S3method(coef,lmList)
S3method(coef,merMod)
S3method(confint,lmList)
+S3method(confint,merMod)
S3method(confint,thpr)
S3method(densityplot,thpr)
S3method(deviance,merMod)
Modified: pkg/lme4/R/AllGeneric.R
===================================================================
--- pkg/lme4/R/AllGeneric.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/AllGeneric.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -45,13 +45,26 @@
##' 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
+##' These are generic functions. At present the only methods are for mixed-effects
##' models of class \code{\linkS4class{merMod}}.
##' @title Check characteristics of models
##' @param x a fitted model.
##' @param ... additional, optional arguments. (None are used in the merMod methods)
##' @return a logical value
##' @seealso getME
+##' @examples
+##' fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+##' gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+##' data = cbpp, family = binomial)
+##' nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
+##' Orange, start = c(Asym = 200, xmid = 725, scal = 350))
+##'
+##' isLMM(fm1)
+##' isGLMM(gm1)
+##' ## check all :
+##' is.MM <- function(x) c(LMM = isLMM(x), GLMM= isGLMM(x), NLMM= isNLMM(x))
+##' stopifnot(cbind(is.MM(fm1), is.MM(gm1), is.MM(nm1))
+##' == diag(rep(TRUE,3)))
##' @export
isREML <- function(x, ...) UseMethod("isREML")
@@ -93,6 +106,15 @@
##' of the same length as the original response.
##' @param ... optional additional parameters. None are used at present.
##' @return an object like \code{x} but fit by maximum likelihood
+##' @examples
+##' ## using refit() to fit each column in a matrix of responses
+##' set.seed(101)
+##' Y <- matrix(rnorm(1000),ncol=10)
+##' res <- list()
+##' d <- data.frame(y=Y[,1],x=rnorm(100),f=rep(1:10,10))
+##' fit1 <- lmer(y~x+(1|f),data=d)
+##' res <- c(fit1,lapply(as.data.frame(Y[,-1]),
+##' refit,object=fit1))
##' @export
refit <- function(object, newresp, ...) UseMethod("refit")
Modified: pkg/lme4/R/GHrule.R
===================================================================
--- pkg/lme4/R/GHrule.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/GHrule.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -160,7 +160,7 @@
}
} else fr <- rbind(fr[rev(seq_len(nr)),], fr)
if (ord > 1L) fr[seq_len(ord %/% 2L), "z"] <- -fr[seq_len(ord %/% 2L), "z"]
- rownames(fr) <- NULL
+ rownames(fr) <- z <- NULL
fr <- within(fr, ldnorm <- dnorm(z, log=TRUE))
if (asMatrix) return(as.matrix(fr))
fr
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/lmer.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -254,7 +254,7 @@
##' ## which is not directly comparable with the nAGQ=0 or nAGQ=1 result.
##' (gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
##' cbpp, binomial, nAGQ = 9))
-##'
+##'
##' ## GLMM with individual-level variability (accounting for overdispersion)
##' ## For this data set the model is the same as one allowing for a period:herd
##' ## interaction, which the plot indicates could be needed.
@@ -324,12 +324,12 @@
form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
X <- model.matrix(form, fr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet
p <- ncol(X)
-
+
## Environment for deviance function. For the optimizers the
## deviance function must be a simple function of a numeric
## parameter. We put all the other information in the
## environment rho which is assigned as the environment of the
- ## deviance function.
+ ## deviance function.
rho <- as.environment(list(verbose=verbose, tolPwrss=tolPwrss))
parent.env(rho) <- parent.frame()
rho$pp <- do.call(merPredD$new,
@@ -354,7 +354,7 @@
control=control,
adj=FALSE, verbose=verbose)
rho$control <- attr(opt,"control")
-
+
rho$nAGQ <- nAGQ
if (nAGQ > 0L) {
rho$nAGQ <- nAGQ
@@ -473,6 +473,7 @@
## dummy
globalVariables <- function(...) {}
}
+if(FALSE)## not ok for roxygen2
globalVariables(c("pp","resp","lp0","pwrssUpdate","compDev",
"baseOffset","GQmat","fac","nlmerAGQ","tolPwrss",
"dpars","verbose"),
@@ -512,43 +513,51 @@
mkdevfun <- function(rho, nAGQ=1L, verbose=0) {
## FIXME: should nAGQ be automatically embedded in rho?
stopifnot(is.environment(rho), is(rho$resp, "lmResp"))
- ff <- NULL
+
+ ## R CMD check "no visible binding for global variable" ...
+ ## MM *preferred* to globalVariables()
+ fac <- pp <- resp <- lp0 <- compDev <- dpars <- baseOffset <- tolPwrss <-
+ pwrssUpdate <- ## <-- even though it's a function below
+ GQmat <- nlmerAGQ <- NULL
+
+ ## The deviance function (to be returned):
+ ff <-
if (is(rho$resp, "lmerResp")) {
- rho$lmer_Deviance <- lmer_Deviance
- ff <- function(theta) .Call(lmer_Deviance, pp$ptr(), resp$ptr(), as.double(theta))
+ rho$lmer_Deviance <- lmer_Deviance
+ function(theta) .Call(lmer_Deviance, pp$ptr(), resp$ptr(), as.double(theta))
} else if (is(rho$resp, "glmResp")) {
- ff <- if (nAGQ == 0L)
- function(theta) {
- resp$updateMu(lp0)
- pp$setTheta(theta)
- pwrssUpdate(pp, resp, 1e-7, GHrule(0L), compDev, verbose)
- }
- else
- function(pars) {
- resp$updateMu(lp0)
- pp$setTheta(as.double(pars[dpars])) # theta is first part of pars
- resp$setOffset(baseOffset + pp$X %*% as.numeric(pars[-dpars]))
- pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose)
- }
+ if (nAGQ == 0L)
+ function(theta) {
+ resp$updateMu(lp0)
+ pp$setTheta(theta)
+ pwrssUpdate(pp, resp, 1e-7, GHrule(0L), compDev, verbose)
+ }
+ else
+ function(pars) {
+ resp$updateMu(lp0)
+ pp$setTheta(as.double(pars[dpars])) # theta is first part of pars
+ resp$setOffset(baseOffset + pp$X %*% as.numeric(pars[-dpars]))
+ pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose)
+ }
} else if (is(rho$resp, "nlsResp")) {
- if (nAGQ < 2L) {
- rho$nlmerLaplace <- nlmerLaplace
- ff <- switch(nAGQ + 1L,
- function(theta)
- .Call(nlmerLaplace, pp$ptr(), resp$ptr(), as.double(theta),
- as.double(u0), beta0, verbose, FALSE, tolPwrss),
- function(pars)
- .Call(nlmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], u0,
- pars[-dpars], verbose, TRUE, tolPwrss))
- } else {
- rho$nlmerAGQ <- nlmerAGQ
- rho$GQmat <- GHrule(nAGQ)
- ff <- function(pars)
- .Call(nlmerAGQ, pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars],
- u0, pars[-dpars], tolPwrss)
- }
+ if (nAGQ < 2L) {
+ rho$nlmerLaplace <- nlmerLaplace
+ switch(nAGQ + 1L,
+ function(theta)
+ .Call(nlmerLaplace, pp$ptr(), resp$ptr(), as.double(theta),
+ as.double(u0), beta0, verbose, FALSE, tolPwrss),
+ function(pars)
+ .Call(nlmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], u0,
+ pars[-dpars], verbose, TRUE, tolPwrss))
+ } else {
+ rho$nlmerAGQ <- nlmerAGQ
+ rho$GQmat <- GHrule(nAGQ)
+ function(pars)
+ .Call(nlmerAGQ, pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars],
+ u0, pars[-dpars], tolPwrss)
+ }
}
- if (is.null(ff)) stop("code not yet written")
+ else stop("code not yet written")
environment(ff) <- rho
ff
}
@@ -1109,12 +1118,14 @@
class(ans) <- "ranef.mer"
}
ans
-}
+}## ranef.merMod
-##' @S3method refit merMod
+
+##' @method refit merMod
+##' @rdname refit
+##' @export
refit.merMod <- function(object, newresp=NULL, ...)
{
-
rr <- object at resp$copy()
if (!is.null(newresp)) {
@@ -1139,10 +1150,8 @@
newresp <- as.numeric(newresp)-1
}
}
-
stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == length(rr$y))
rr$setResp(newresp)
-
}
pp <- object at pp$copy()
@@ -1184,7 +1193,13 @@
object at frame, getCall(object))
}
-##' @S3method refitML merMod
+##-- BUG in roxygen2: If we use @S3method instead of @method,
+##-- the \usage{ ... } will have
+##-- refitML.merMod(..) instead of \method{refitML}{mermod}(..)
+##' @param optimizer a string indicating the optimizer to be used.
+##' @method refitML merMod
+##' @rdname refitML
+##' @export
refitML.merMod <- function (x, optimizer="bobyqa", ...) {
## FIXME: optimizer is set to 'bobyqa' for back-compatibility, but that's not
## consistent with lmer (default NM). Should be based on internally stored 'optimizer' value
@@ -1574,7 +1589,7 @@
##' @keywords utilities
##' @examples
##'
-##' ## shows many methods you should consider *before* getME():
+##' ## shows many methods you should consider *before* using getME():
##' methods(class = "merMod")
##'
##' (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
@@ -1586,6 +1601,12 @@
##'
##' ## All that can be accessed [potentially ..]:
##' (nmME <- eval(formals(getME)$name))
+##' \dontshow{
+##' ## internal consistency check ensuring that all work:
+##' ## "try(.)" because some are not yet implemented:
+##' str(parts <- sapply(nmME, function(nm) try(getME(fm1, nm)),
+##' simplify=FALSE))
+##' }% dont..
##'
##' @export
getME <- function(object,
Modified: pkg/lme4/R/mcmcsamp.R
===================================================================
--- pkg/lme4/R/mcmcsamp.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/mcmcsamp.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -1,3 +1,4 @@
+if(FALSE) ## C++ code in ../src/mcmcsamp.cpp -- is also #ifdef 0
# @S3method mcmcsamp merMod
mcmcsamp.merMod <- function(object, n=1L, verbose=FALSE, saveb=FALSE, ...) {
n <- max(1L, as.integer(n)[1])
Modified: pkg/lme4/R/plot.R
===================================================================
--- pkg/lme4/R/plot.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/plot.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -140,9 +140,8 @@
## Diagnostic plots based on residuals and/or fitted values
{
object <- x
- if (!inherits(form, "formula")) {
- stop("\"Form\" must be a formula")
- }
+ if (!inherits(form, "formula"))
+ stop("\"form\" must be a formula")
## constructing data
## can I get away with using object at frame???
allV <- all.vars(asOneFormula(form, id, idLabels))
@@ -232,7 +231,9 @@
stop("\"Id\" can only be a formula or numeric.")
)
if (is.null(idLabels)) {
- idLabels <- getGroups(object)
+ ## FIXME: Have no example yet, where this is triggered...
+ message("**** getGroups() case in plot.merMod() ****")
+ idLabels <- getGroupsFormula(object)
if (length(idLabels) == 0) idLabels <- 1:object$dims$N
idLabels <- as.character(idLabels)
} else {
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/R/profile.R 2012-06-16 17:57:36 UTC (rev 1767)
@@ -30,19 +30,25 @@
##' ## 0.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
##' ## This is slower because of the Nelder-Mead optimizer but using bobyqa, the default,
##' ## produces a warning.
-##' system.time( tpr <- profile(fm01ML, optimizer="Nelder_Mead") )
+##' system.time( tpr <- profile(fm01ML, optimizer="Nelder_Mead") )
+##' system.time( tpr2 <- profile(fm01ML, optimizer="bobyqa") )
##' (confint(tpr) -> CIpr)
-##' print(xyplot(tpr))
+##' stopifnot(all.equal(CIpr, confint(tpr2), tol= 1e-11))
+##' xyplot(tpr)
+##' densityplot(tpr, main="densityplot( profile(lmer(..)) )")
+##' splom(tpr)
+##'
##' tpr2 <- profile(fm01ML, which=1:2, optimizer="Nelder_Mead") ## 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)
+##' splom(pr4)% FIXME! currently fails!
##' }
-##' @importFrom splines backSpline interpSpline periodicSpline
+##' ##' @importFrom splines backSpline interpSpline periodicSpline
##' @importFrom stats profile
##' @method profile merMod
##' @export
@@ -308,10 +314,10 @@
} ## if isLMM
ans <- do.call(rbind, ans)
+ row.names(ans) <- NULL
ans$.par <- factor(ans$.par)
attr(ans, "forward") <- forspl
attr(ans, "backward") <- bakspl
- row.names(ans) <- NULL
class(ans) <- c("thpr", "data.frame")
ans
}
@@ -528,27 +534,33 @@
## FIXME: make bootMer more robust; make profiling more robust;
## more warnings; documentation ...
+## --> then @export (we have "surprising arguments")
+
##' @importFrom stats confint
##' @S3method confint merMod
-confint.merMod <- function(object, parm, method=c("profile","quadratic","bootMer"), level = 0.95, zeta, nsim=500, boot.type="perc", ...)
+confint.merMod <- function(object, parm, level = 0.95,
+ method=c("profile","quadratic","bootMer"),
+ zeta, nsim=500, boot.type="perc", ...)
{
method <- match.arg(method)
- if (method=="profile") {
- if (!missing(parm)) {
- pp <- profile(object,which=parm)
- } else {
- pp <- profile(object)
- }
- return(confint(pp,level=level,zeta=zeta,...))
- } else if (method=="quadratic") {
- stop("stub")
- ## only give confidence intervals on fixed effects?
- ## or warn??
- } else if (method=="bootMer") {
- message("Computing bootstrap confidence intervals ...")
- bb <- bootMer(object, fixef, nsim=nsim)
- boot::boot.ci(bb,type=boot.type,conf=level)
- }
+ switch(method,
+ "profile" =
+ {
+ pp <- if(missing(parm)) profile(object) else profile(object, which=parm)
+ confint(pp,level=level,zeta=zeta,...)
+ },
+ "quadratic" =
+ {
+ stop("stub")
+ ## only give confidence intervals on fixed effects?
+ ## or warn??
+ },
+ "bootMer" =
+ {
+ message("Computing bootstrap confidence intervals ...")
+ bb <- bootMer(object, fixef, nsim=nsim)
+ boot::boot.ci(bb,type=boot.type,conf=level)
+ })
}
## Convert x-cosine and y-cosine to average and difference.
@@ -562,7 +574,7 @@
{
a <- (xc + yc)/2
d <- (xc - yc)
- cbind(ifelse(d > 0, a, -a), abs(d))
+ cbind(sign(d)* a, abs(d))
}
## convert d versus a (as an xyVector) and level to a matrix of taui and tauj
@@ -602,14 +614,20 @@
}
-## Profile pairs plot
-## Contours are for the marginal two-dimensional regions (i.e. using
-## df = 2)
+##' Draws profile pairs plots. Contours are for the marginal
+##' two-dimensional regions (i.e. using df = 2).
+##'
+##' @title Profile pairs plot
+##' @param x the result of \code{\link{profile}} (or very similar structure)
+##' @param data unused - only for compatibility with generic
+##' @param levels the contour levels to be shown; usually derived from \code{conf}
+##' @param conf numeric vector of confidence levels to be shown as contours
+##' @param ... further arguments passed to \code{\link{splom}}
##' @importFrom grid gpar viewport
##' @importFrom lattice splom
-##' @S3method splom thpr
-splom.thpr <-
- function (x, data, ## unused - only for compatibility with generic
+##' @method splom thpr
+##' @export
+splom.thpr <- function (x, data,
levels = sqrt(qchisq(pmax.int(0, pmin.int(1, conf)), 2)),
conf = c(50, 80, 90, 95, 99)/100, ...)
{
@@ -801,6 +819,7 @@
x
}
+### FIXME: Not exported and nowhere used
## Transform a profile from the standard deviation parameters to the variance
##
## @title Transform to variance component scale
@@ -861,15 +880,18 @@
fr
}
-## Densityplot method for a mixed-effects model profile
+##' Densityplot method for a mixed-effects model profile
##
-## @title densityplot from a mixed-effects profile
-## @param x a mixed-effects profile
-## @param data not used - for compatibility with generic
-## @param ... optional arguments to densityplot
-## @return a density plot
+##' @title densityplot from a mixed-effects profile
+##' @param x a mixed-effects profile
+##' @param data not used - for compatibility with generic
+##' @param ... optional arguments to \code{\link[lattice]{densityplot}()}
+##' from package \pkg{lattice}.
+##' @return a density plot
+##' @examples ## see example("profile.merMod")
##' @importFrom lattice densityplot
-##' @S3method densityplot thpr
+##' @method densityplot thpr
+##' @export
densityplot.thpr <- function(x, data, ...) {
ll <- c(list(...),
list(x=density ~ pval|pnm,
@@ -907,16 +929,19 @@
## convert profile to data frame, adding a .focal parameter to simplify lattice/ggplot plotting
##' @method as.data.frame thpr
+##' @param x the result of \code{\link{profile}} (or very similar structure)
##' @export
+##' @rdname profile-methods
as.data.frame.thpr <- function(x,...) {
class(x) <- "data.frame"
m <- as.matrix(x[,seq(ncol(x))-1]) ## omit .par
- x[[".focal"]] <- m[cbind(seq(nrow(x)),match(x$.par,names(x)))]
- x[[".par"]] <- factor(x[[".par"]],levels=unique(as.character(x[[".par"]]))) ## restore order
+ x.p <- x[[".par"]]
+ x[[".focal"]] <- m[cbind(seq(nrow(x)),match(x.p,names(x)))]
+ x[[".par"]] <- factor(x.p, levels=unique(as.character(x.p))) ## restore order
x
}
-if (FALSE) {
+if (FALSE) { ## FIXME Ben , e.g. mv to ../testsx/ ?
source("~/R/misc/slice.R")
dd <- lme4:::devfun2(nm1,useSc=TRUE)
s1 <- slice(attr(dd,"optimum"),dd,dim=2)
Modified: pkg/lme4/man/VerbAgg.Rd
===================================================================
--- pkg/lme4/man/VerbAgg.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/VerbAgg.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -41,7 +41,7 @@
xtabs(~ btype + resp, VerbAgg)
round(100 * ftable(prop.table(xtabs(~ situ + mode + resp, VerbAgg), 1:2), 1))
person <- unique(subset(VerbAgg, select = c(id, Gender, Anger)))
-if (require(lattice)) { # is this necessary when the package depends on lattice?
+if(require(lattice)) { # not needed as long as the lme4 depends on lattice
densityplot(~ Anger, person, groups = Gender, auto.key = list(columns = 2),
xlab = "Trait Anger score (STAXI)")
}
Added: pkg/lme4/man/densityplot.thpr.Rd
===================================================================
--- pkg/lme4/man/densityplot.thpr.Rd (rev 0)
+++ pkg/lme4/man/densityplot.thpr.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -0,0 +1,25 @@
+\name{densityplot.thpr}
+\alias{densityplot.thpr}
+\title{densityplot from a mixed-effects profile}
+\usage{
+ \method{densityplot}{thpr} (x, data, ...)
+}
+\arguments{
+ \item{x}{a mixed-effects profile}
+
+ \item{data}{not used - for compatibility with generic}
+
+ \item{...}{optional arguments to
+ \code{\link[lattice]{densityplot}()} from package
+ \pkg{lattice}.}
+}
+\value{
+ a density plot
+}
+\description{
+ Densityplot method for a mixed-effects model profile
+}
+\examples{
+## see example("profile.merMod")
+}
+
Modified: pkg/lme4/man/getME.Rd
===================================================================
--- pkg/lme4/man/getME.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/getME.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -68,7 +68,7 @@
\code{\link{vcov}}, etc.
}
\examples{
-## shows many methods you should consider *before* getME():
+## shows many methods you should consider *before* using getME():
methods(class = "merMod")
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
@@ -80,6 +80,12 @@
## All that can be accessed [potentially ..]:
(nmME <- eval(formals(getME)$name))
+\dontshow{
+## internal consistency check ensuring that all work:
+## "try(.)" because some are not yet implemented:
+str(parts <- sapply(nmME, function(nm) try(getME(fm1, nm)),
+ simplify=FALSE))
+}\% dont..
}
\seealso{
\code{\link{getCall}()}, More standard methods for mer
Modified: pkg/lme4/man/isREML.Rd
===================================================================
--- pkg/lme4/man/isREML.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/isREML.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -31,10 +31,24 @@
NLMMs).
}
\details{
- This is a generic function. At present the only methods
+ These are generic functions. At present the only methods
are for mixed-effects models of class
\code{\linkS4class{merMod}}.
}
+\examples{
+fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ data = cbpp, family = binomial)
+nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
+ Orange, start = c(Asym = 200, xmid = 725, scal = 350))
+
+isLMM(fm1)
+isGLMM(gm1)
+## check all :
+is.MM <- function(x) c(LMM = isLMM(x), GLMM= isGLMM(x), NLMM= isNLMM(x))
+stopifnot(cbind(is.MM(fm1), is.MM(gm1), is.MM(nm1))
+ == diag(rep(TRUE,3)))
+}
\seealso{
getME
}
Modified: pkg/lme4/man/mkdevfun.Rd
===================================================================
--- pkg/lme4/man/mkdevfun.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/mkdevfun.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -2,7 +2,7 @@
\alias{mkdevfun}
\title{Create a deviance evaluation function from a predictor and a response module}
\usage{
- mkdevfun(rho, nAGQ = 1L)
+ mkdevfun(rho, nAGQ = 1L, verbose = 0)
}
\arguments{
\item{rho}{an environment containing \code{pp}, a
@@ -15,6 +15,8 @@
that both the fixed-effects parameters and the random
effects are optimized by the iteratively reweighted least
squares algorithm.}
+
+ \item{verbose}{Logical: print verbose output?}
}
\value{
A function of one numeric argument.
Modified: pkg/lme4/man/profile-methods.Rd
===================================================================
--- pkg/lme4/man/profile-methods.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/profile-methods.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -1,5 +1,6 @@
\docType{methods}
\name{profile-methods}
+\alias{as.data.frame.thpr}
\alias{profile.merMod}
\alias{profile-methods}
\title{Methods for profile() of [ng]lmer fitted models}
@@ -8,6 +9,8 @@
alphamax = 0.01, maxpts = 100, delta = cutoff/8,
verbose = 0, devtol = 1e-09, maxmult = 10,
startmethod = "prev", optimizer = "bobyqa", ...)
+
+ \method{as.data.frame}{thpr} (x, ...)
}
\arguments{
\item{fitted}{a fitted model, e.g., the result of
@@ -51,6 +54,9 @@
\item{\dots}{potential further arguments for
\code{profile} methods.}
+
+ \item{x}{the result of \code{\link{profile}} (or very
+ similar structure)}
}
\description{
Methods for function \code{\link{profile}} (package
@@ -67,18 +73,25 @@
## 0.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
## This is slower because of the Nelder-Mead optimizer but using bobyqa, the default,
## produces a warning.
-system.time( tpr <- profile(fm01ML, optimizer="Nelder_Mead") )
+system.time( tpr <- profile(fm01ML, optimizer="Nelder_Mead") )
+system.time( tpr2 <- profile(fm01ML, optimizer="bobyqa") )
(confint(tpr) -> CIpr)
-print(xyplot(tpr))
+stopifnot(all.equal(CIpr, confint(tpr2), tol= 1e-11))
+xyplot(tpr)
+densityplot(tpr, main="densityplot( profile(lmer(..)) )")
+splom(tpr)
+
tpr2 <- profile(fm01ML, which=1:2, optimizer="Nelder_Mead") ## 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)
+splom(pr4)\% FIXME! currently fails!
}
+##'
}
\seealso{
For (more expensive) alternative confidence intervals:
Modified: pkg/lme4/man/refit.Rd
===================================================================
--- pkg/lme4/man/refit.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/refit.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -1,8 +1,11 @@
\name{refit}
\alias{refit}
+\alias{refit.merMod}
\title{Refit a model by maximum likelihood criterion}
\usage{
refit(object, newresp, ...)
+
+ \method{refit}{merMod} (object, newresp = NULL, ...)
}
\arguments{
\item{object}{a fitted model, usually of class
@@ -29,4 +32,14 @@
creation of the model representation and goes directly to
the optimization step.
}
+\examples{
+## using refit() to fit each column in a matrix of responses
+set.seed(101)
+Y <- matrix(rnorm(1000),ncol=10)
+res <- list()
+d <- data.frame(y=Y[,1],x=rnorm(100),f=rep(1:10,10))
+fit1 <- lmer(y~x+(1|f),data=d)
+res <- c(fit1,lapply(as.data.frame(Y[,-1]),
+ refit,object=fit1))
+}
Modified: pkg/lme4/man/refitML.Rd
===================================================================
--- pkg/lme4/man/refitML.Rd 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/man/refitML.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -1,8 +1,11 @@
\name{refitML}
\alias{refitML}
+\alias{refitML.merMod}
\title{Refit a model by maximum likelihood criterion}
\usage{
refitML(x, ...)
+
+ \method{refitML}{merMod} (x, optimizer = "bobyqa", ...)
}
\arguments{
\item{x}{a fitted model, usually of class
@@ -11,6 +14,9 @@
\item{...}{optional additional parameters. None are used
at present.}
+
+ \item{optimizer}{a string indicating the optimizer to be
+ used.}
}
\value{
an object like \code{x} but fit by maximum likelihood
Added: pkg/lme4/man/splom.thpr.Rd
===================================================================
--- pkg/lme4/man/splom.thpr.Rd (rev 0)
+++ pkg/lme4/man/splom.thpr.Rd 2012-06-16 17:57:36 UTC (rev 1767)
@@ -0,0 +1,28 @@
+\name{splom.thpr}
+\alias{splom.thpr}
+\title{Profile pairs plot}
+\usage{
+ \method{splom}{thpr} (x, data,
+ levels = sqrt(qchisq(pmax.int(0, pmin.int(1, conf)), 2)),
+ conf = c(50, 80, 90, 95, 99)/100, ...)
+}
+\arguments{
+ \item{x}{the result of \code{\link{profile}} (or very
+ similar structure)}
+
+ \item{data}{unused - only for compatibility with generic}
+
+ \item{levels}{the contour levels to be shown; usually
+ derived from \code{conf}}
+
+ \item{conf}{numeric vector of confidence levels to be
+ shown as contours}
+
+ \item{...}{further arguments passed to
+ \code{\link{splom}}}
+}
+\description{
+ Draws profile pairs plots. Contours are for the marginal
+ two-dimensional regions (i.e. using df = 2).
+}
+
Modified: pkg/lme4/tests/nlmer.Rout.save
===================================================================
--- pkg/lme4/tests/nlmer.Rout.save 2012-06-16 14:33:07 UTC (rev 1766)
+++ pkg/lme4/tests/nlmer.Rout.save 2012-06-16 17:57:36 UTC (rev 1767)
@@ -1,5 +1,5 @@
-R version 2.15.0 beta (2012-03-17 r58782)
+R version 2.15.1 beta (2012-06-14 r59562) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -18,6 +18,7 @@
> library(lme4)
Loading required package: lattice
+Loading required package: Matrix
>
> allEQ <- function(x,y, tolerance = 4e-4, ...)
+ all.equal.numeric(x,y, tolerance=tolerance, ...)
@@ -49,7 +50,7 @@
scal 0.362 0.762
> fixef(nm1)
Asym xmid scal
-192.0533 727.9074 348.0737
+192.0534 727.9074 348.0738
>
> ## 'Theoph' Data modeling
> Th.start <- c(lKe = -2.5, lKa = 0.5, lCl = -3)
@@ -58,7 +59,7 @@
+ (lKe+lKa+lCl|Subject),
+ Theoph, start = Th.start, tolPwrss=1e-8))
user system elapsed
- 5.664 0.000 5.682
+ 5.480 0.000 5.497
> print(nm2, corr=FALSE)
Nonlinear mixed model fit by maximum likelihood ['nlmerMod']
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) ~ (lKe + lKa + lCl | Subject)
@@ -77,15 +78,15 @@
Fixed effects:
Estimate Std. Error t value
-lKe -2.44546 0.06257 -39.08
-lKa 0.47635 0.19395 2.46
-lCl -3.21454 0.08100 -39.68
+lKe -2.44551 0.06257 -39.08
+lKa 0.47646 0.19395 2.46
+lCl -3.21462 0.08100 -39.69
>
> system.time(nm3 <- nlmer(conc ~ SSfol(Dose, Time,lKe, lKa, lCl) ~
+ (lKe|Subject) + (lKa|Subject) + (lCl|Subject),
+ Theoph, start = Th.start))
user system elapsed
- 2.644 0.000 2.655
+ 2.728 0.004 2.739
> print(nm3, corr=FALSE)
Nonlinear mixed model fit by maximum likelihood ['nlmerMod']
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) ~ (lKe | Subject) + (lKa | Subject) + (lCl | Subject)
@@ -104,15 +105,15 @@
Fixed effects:
Estimate Std. Error t value
-lKe -2.46572 0.05187 -47.54
-lKa 0.48088 0.20000 2.40
-lCl -3.23074 0.05955 -54.26
+lKe -2.46537 0.05187 -47.53
+lKa 0.48103 0.20000 2.41
+lCl -3.23045 0.05955 -54.25
>
> ## dropping lKe from random effects:
> system.time(nm4 <- nlmer(conc ~ SSfol(Dose, Time,lKe, lKa, lCl) ~ (lKa+lCl|Subject),
+ Theoph, start = Th.start, tolPwrss=1e-8))
user system elapsed
- 1.640 0.000 1.646
+ 1.624 0.000 1.630
> print(nm4, corr=FALSE)
Nonlinear mixed model fit by maximum likelihood ['nlmerMod']
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) ~ (lKa + lCl | Subject)
@@ -131,7 +132,7 @@
Fixed effects:
Estimate Std. Error t value
lKe -2.46592 0.05188 -47.53
-lKa 0.48282 0.20004 2.41
+lKa 0.48283 0.20004 2.41
lCl -3.22897 0.05953 -54.24
>
> system.time(nm5 <- nlmer(conc ~ SSfol(Dose, Time,lKe, lKa, lCl) ~
@@ -139,7 +140,7 @@
+ Theoph,
+ start = Th.start, tolPwrss=1e-8))
user system elapsed
- 0.876 0.000 0.880
+ 0.872 0.000 0.874
> print(nm5, corr=FALSE)
Nonlinear mixed model fit by maximum likelihood ['nlmerMod']
Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) ~ (lKa | Subject) + (lCl | Subject)
@@ -157,9 +158,9 @@
Fixed effects:
Estimate Std. Error t value
-lKe -2.46579 0.05185 -47.56
-lKa 0.48262 0.20064 2.41
-lCl -3.23053 0.05957 -54.23
+lKe -2.46584 0.05185 -47.56
+lKa 0.48229 0.20064 2.40
+lCl -3.23055 0.05957 -54.23
>
> if (require("PKPDmodels")) {
+ oral1cptSdlkalVlCl <-
@@ -187,10 +188,10 @@
Fixed effects:
Estimate Std. Error t value
-lV -0.77187 0.04208 -18.34
-lka 0.47106 0.19864 2.37
-lCl -3.21992 0.07995 -40.28
+lV -0.77190 0.04208 -18.35
+lka 0.47101 0.19864 2.37
+lCl -3.21971 0.07995 -40.27
>
> proc.time()
user system elapsed
- 20.041 0.092 20.238
+ 19.785 0.080 19.919
More information about the Lme4-commits
mailing list