[Lme4-commits] r1584 - in pkg/lme4Eigen: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 9 19:39:11 CET 2012
Author: dmbates
Date: 2012-02-09 19:39:11 +0100 (Thu, 09 Feb 2012)
New Revision: 1584
Modified:
pkg/lme4Eigen/R/lme4Eigen-package.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/predict.R
pkg/lme4Eigen/R/utilities.R
pkg/lme4Eigen/man/Nelder_Mead.Rd
pkg/lme4Eigen/man/VarCorr.Rd
pkg/lme4Eigen/man/cbpp.Rd
pkg/lme4Eigen/man/findbars.Rd
pkg/lme4Eigen/man/getME.Rd
Log:
Update documentation and exports
Modified: pkg/lme4Eigen/R/lme4Eigen-package.R
===================================================================
--- pkg/lme4Eigen/R/lme4Eigen-package.R 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/R/lme4Eigen-package.R 2012-02-09 18:39:11 UTC (rev 1584)
@@ -87,18 +87,14 @@
##' ## response as a matrix
##' (m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
##' cbpp, binomial, nAGQ=25L))
-##' dput(unname(fixef(m1)))
-##' dput(unname(ranef(m1, drop=TRUE)[[1]]))
##' ## response as a vector of probabilities and usage of argument "weights"
##' m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
##' cbpp, binomial, nAGQ=25L)
-##' dput(unname(fixef(m1p)))
-##' dput(unname(ranef(m1p, drop=TRUE)[[1]]))
##' ## Confirm that these are equivalent:
##' stopifnot(all.equal(fixef(m1), fixef(m1p), tol = 1e-5),
##' all.equal(ranef(m1), ranef(m1p), tol = 1e-5),
##' TRUE)
-##'
+##' ## Can this section be moved to a test file? I don't think it belongs in an example. DB
##' for(m in c(m1, m1p)) {
##' cat("-------\n\nCall: ",
##' paste(format(getCall(m)), collapse="\n"), "\n")
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-09 18:39:11 UTC (rev 1584)
@@ -624,31 +624,30 @@
##' \code{x}.
##'
##' @param x fitted \code{*lmer()} model, see \code{\link{lmer}},
-##' \code{\link{glmer}}, etc.
+##' \code{\link{glmer}}, etc.
##' @param FUN a \code{\link{function}(x)}, computating the \emph{statistic} of
-##' interest, which must be a numeric vector, possibly named.
+##' interest, which must be a numeric vector, possibly named.
##' @param nsim number of simulations, positive integer; the bootstrap \eqn{B}
-##' (or \eqn{R}).
+##' (or \eqn{R}).
##' @param seed optional argument to \code{\link{set.seed}}.
##' @param use.u logical, indicating, if the spherized random effects should be
-##' simulated / bootstrapped as well. If \code{FALSE}, they are not changed,
-##' and all inference is conditional on these.
+##' simulated / bootstrapped as well. If \code{FALSE}, they are not changed,
+##' and all inference is conditional on these.
##' @param verbose logical indicating if progress should print output
##' @param control an optional \code{\link{list}}, to be passed to the minimizer
-##' (of the log-likelihood, or RE likelihood), which is currently set to
-##' \code{\link[minqa]{bobyqa}} in package \pkg{minqa}.
+##' (of the log-likelihood, or RE likelihood), which is currently set to
+##' \code{\link[minqa]{bobyqa}} in package \pkg{minqa}.
##' @return an object of S3 \code{\link{class}} \code{"boot"}, compatible with
-##' \pkg{boot} package's \code{boot()} result.
+##' \pkg{boot} package's \code{boot()} result.
##' @seealso For inference, including confidence intervals,
-##' \code{\link{profile-methods}}.
+##' \code{\link{profile-methods}}.
##'
-##' \code{\link[boot]{boot}()}, and then \code{\link[boot]{boot.ci}} from
-##' package \pkg{boot}.
+##' \code{\link[boot]{boot}()}, and then \code{\link[boot]{boot.ci}} from
+##' package \pkg{boot}.
##' @references Davison, A.C. and Hinkley, D.V. (1997) \emph{Bootstrap Methods
-##' and Their Application}. Cambridge University Press.
+##' and Their Application}. Cambridge University Press.
##' @keywords models htest
##' @examples
-##'
##' fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
##' ## see ?"profile-methods"
##' mySumm <- function(.) { s <- sigma(.)
@@ -678,8 +677,7 @@
##' }
##' @export
bootMer <- function(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
- verbose = FALSE, control = list())
-{
+ verbose = FALSE, control = list()) {
stopifnot((nsim <- as.integer(nsim[1])) > 0,
is(x, "merMod"), is(x at resp, "lmerResp"))
FUN <- match.fun(FUN)
@@ -708,9 +706,9 @@
sigm.x <- sigma(x)
## Here, and below ("optimize"/"bobyqa") using the "logic" of lmer() itself:
-## lmer..Update <- if(is(x, "lmerSp")) lmerSpUpdate else lmerDeUpdate
-# devfun <- mkdevfun(x)
-## oneD <- length(x at re@theta) < 2
+ ## lmer..Update <- if(is(x, "lmerSp")) lmerSpUpdate else lmerDeUpdate
+ ## devfun <- mkdevfun(x)
+ ## oneD <- length(x at re@theta) < 2
theta0 <- getME(x,"theta")
## just for the "boot" result -- TODOmaybe drop
mle <- list(beta = beta, theta = theta0, sigma = sigm.x)
@@ -734,16 +732,16 @@
## devfun(0) # -> theta := 0 and update the rest
## } else {
opt <- bobyqa(theta0, mkdevfun(resp, x at pp), x at lower, control = control)
-## xx <- updateMod(x, opt$par, opt$fval)
- ## FIXME: also here, prefer \hat\sigma^2 == 0 (exactly)
-## }
+ ## xx <- updateMod(x, opt$par, opt$fval)
+ ## FIXME: also here, prefer \hat\sigma^2 == 0 (exactly)
+ ## }
foo <- tryCatch(FUN(xx), error = function(e)e)
if(verbose) { cat(sprintf("%5d :",i)); str(foo) }
t.star[,i] <- if (inherits(foo, "error")) NA else foo
}
rownames(t.star) <- names(t0)
-## boot() ends with the equivalent of
+ ## boot() ends with the equivalent of
## structure(list(t0 = t0, t = t.star, R = R, data = data, seed = seed,
## statistic = statistic, sim = sim, call = call,
## ran.gen = ran.gen, mle = mle),
@@ -1261,9 +1259,17 @@
}
##' @importFrom stats simulate
-##' @S3method simulate merMod
-##' @aliases simulate.merMod
-##' @param use.u (logical) generate new random-effects values (FALSE) or generate a simulation condition on the current random-effects estimates (TRUE)?
+NULL
+##' Simulate responses from the model represented by a fitted model object
+##'
+##' @title Simulate responses from a \code{\linkS4class{merMod}} object
+##' @param object a fitted model object
+##' @param nsim positive integer scalar - the number of responses to simulate
+##' @param seed an optional seed to be used in \code{set.seed} immediately
+##' before the simulation so as to to generate a reproducible sample.
+##' @param use.u (logical) generate new random-effects values (FALSE) or
+##' generate a simulation condition on the current random-effects estimates (TRUE)?
+##' @param ... optional additional arguments, none are used at present
##' @examples
##' ## test whether fitted models are consistent with the
##' ## observed number of zeros in CBPP data set:
@@ -1273,7 +1279,8 @@
##' zeros <- sapply(gg,function(x) sum(x[,"incidence"]==0))
##' plot(table(zeros))
##' abline(v=sum(cbpp$incidence==0),col=2)
-
+##' @method simulate merMod
+##' @export
simulate.merMod <- function(object, nsim = 1, seed = NULL, use.u = FALSE, ...) {
stopifnot((nsim <- as.integer(nsim[1])) > 0,
is(object, "merMod"))
@@ -1718,6 +1725,7 @@
}
## FIXME: should ... go to formatVC or to print ... ?
+##' @S3method print VarCorr.merMod
print.VarCorr.merMod <- function(x,digits = max(3, getOption("digits") - 2), ...) {
print(formatVC(x, digits = digits, useScale = attr(x,"useSc"), ...),quote=FALSE)
}
Modified: pkg/lme4Eigen/R/predict.R
===================================================================
--- pkg/lme4Eigen/R/predict.R 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/R/predict.R 2012-02-09 18:39:11 UTC (rev 1584)
@@ -1,77 +1,97 @@
+##' \code{\link{predict}} method for \code{\linkS4class{merMod}} objects
+##'
+##' @title Predictions from a model at new data values
+##' @param object a fitted model object
##' @param newdata data frame for which to evaluate predictions
-##' @param REform formula for random effects to include. If NULL, include all random effects; if NA, include no random effects
-##' @param allow.new.levels (logical) if FALSE, then any new levels detected in \code{newdata} will trigger an error; if TRUE, then the prediction will use the unconditional (population-level) values for data with previously unobserved levels
+##' @param REform formula for random effects to include. If NULL,
+##' include all random effects; if NA, include no random effects
+##' @param terms a \code{\link{terms}} object - not used at present
+##' @param type character string - either \code{"link"}, the default,
+##' or \code{"response"} indicating the type of prediction object returned.
+##' @param allow.new.levels (logical) if FALSE, then any new levels
+##' detected in \code{newdata} will trigger an error; if TRUE, then
+##' the prediction will use the unconditional (population-level)
+##' values for data with previously unobserved levels
+##' @param ... optional additional parameters. None are used at present.
+##' @return a numeric vector of predicted values
##' @note explain why there is no option for computing standard errors of predictions!
##' @note offsets not yet handled
-
-## FIXME: take some of the stuff from tests/predict.R and incorporate it as examples
-
+##' @examples
+##' (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 |herd), cbpp, binomial))
+##' str(p0 <- predict(gm1)) # fitted values
+##' str(p1 <- predict(gm1,REform=NA)) # fitted values, unconditional (level-0)
+##' newdata <- with(cbpp, expand.grid(period=unique(period), herd=unique(herd)))
+##' str(p2 <- predict(gm1,newdata)) # new data, all RE
+##' str(p3 <- predict(gm1,newdata,REform=NA)) # new data, level-0
+##' str(p4 <- predict(gm1,newdata,REform=~(1|herd))) # explicitly specify RE
+##' @method predict merMod
+##' @export
predict.merMod <- function(object, newdata=NULL, REform=NULL,
terms=NULL, type=c("link","response"),
allow.new.levels=FALSE, ...) {
- ## FIXME: appropriate names for result vector?
- if (any(getME(object,"offset")!=0)) stop("offsets not handled yet") ## FIXME
- type <- match.arg(type)
- if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
- X_orig <- getME(object, "X")
- ## FIXME/WARNING: how do we do this in an eval-safe way???
- form_orig <- eval(object at call$formula,parent.frame())
- if (is.null(newdata) && is.null(REform)) {
- if (is(object at resp,"lmerResp")) return(fitted(object))
- return(switch(type,response=object at resp$mu, ## fitted(object),
- link=object at resp$eta)) ## fixme: getME() ?
- } else {
- if (is.null(newdata)) {
- X <- X_orig
+ ## FIXME: appropriate names for result vector?
+ if (any(getME(object,"offset")!=0)) stop("offsets not handled yet") ## FIXME
+ type <- match.arg(type)
+ if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
+ X_orig <- getME(object, "X")
+ ## FIXME/WARNING: how do we do this in an eval-safe way???
+ form_orig <- eval(object at call$formula,parent.frame())
+ if (is.null(newdata) && is.null(REform)) {
+ if (is(object at resp,"lmerResp")) return(fitted(object))
+ return(switch(type,response=object at resp$mu, ## fitted(object),
+ link=object at resp$eta)) ## fixme: getME() ?
} else {
- form <- form_orig
- form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
- RHS <- form[-2]
- X <- model.matrix(RHS, newdata, contrasts.arg=attr(X_orig,"contrasts"))
+ if (is.null(newdata)) {
+ X <- X_orig
+ } else {
+ form <- form_orig
+ form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
+ RHS <- form[-2]
+ X <- model.matrix(RHS, newdata, contrasts.arg=attr(X_orig,"contrasts"))
+ }
+ pred <- drop(X %*% fixef(object))
+ if (is.null(REform)) {
+ REform <- form_orig[-2]
+ }
+ ## FIXME: ??? can't apply is.na() to a 'language' object?
+ ## what's the appropriate test?
+ if (is.language(REform)) {
+ re <- ranef(object)
+ ##
+ ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
+ new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],levels)
+ re_x <- mapply(function(x,n) {
+ if (any(!new_levels[[n]] %in% rownames(x))) {
+ if (!allow.new.levels) stop("new levels detected in newdata")
+ ## create an all-zero data frame corresponding to the new set of levels ...
+ newx <- as.data.frame(matrix(0,nrow=length(new_levels[[n]]),ncol=ncol(x),
+ dimnames=list(new_levels[[n]],names(x))))
+ ## then paste in the matching RE values from the original fit/set of levels
+ newx[rownames(x),] <- x
+ x <- newx
+ }
+ x
+ },
+ re,names(re),SIMPLIFY=FALSE)
+ ## separate random effects from orig model into individual columns
+ re_List <- do.call(c,lapply(re_x,as.list))
+ re_names <- names(re_List)
+ z_names <- mapply(paste,names(ReTrms$cnms),ReTrms$cnms,MoreArgs=list(sep="."))
+ ## pick out random effects values that correspond to
+ ## random effects incorporated in REform ...
+ ## FIXME: more tests for possible things going wrong here?
+ m <- match(z_names,re_names)
+ if (any(is.na(m)))
+ stop("random effects specified in REform that were not present in original model")
+ re_new <- unlist(re_List[m])
+ pred <- pred + drop(as.matrix(re_new %*% ReTrms$Zt))
+ } ## REform provided
+ } ## predictions with new data or new REform
+ ## FIXME: would like to have an isGLMM() accessor for merMod objects?
+ if (is(object at resp,"glmResp") && type=="response") {
+ pred <- object at resp$family$linkinv(pred)
}
- pred <- drop(X %*% fixef(object))
- if (is.null(REform)) {
- REform <- form_orig[-2]
- }
- ## FIXME: ??? can't apply is.na() to a 'language' object?
- ## what's the appropriate test?
- if (is.language(REform)) {
- re <- ranef(object)
- ##
- ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
- new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],levels)
- re_x <- mapply(function(x,n) {
- if (any(!new_levels[[n]] %in% rownames(x))) {
- if (!allow.new.levels) stop("new levels detected in newdata")
- ## create an all-zero data frame corresponding to the new set of levels ...
- newx <- as.data.frame(matrix(0,nrow=length(new_levels[[n]]),ncol=ncol(x),
- dimnames=list(new_levels[[n]],names(x))))
- ## then paste in the matching RE values from the original fit/set of levels
- newx[rownames(x),] <- x
- x <- newx
- }
- x
- },
- re,names(re),SIMPLIFY=FALSE)
- ## separate random effects from orig model into individual columns
- re_List <- do.call(c,lapply(re_x,as.list))
- re_names <- names(re_List)
- z_names <- mapply(paste,names(ReTrms$cnms),ReTrms$cnms,MoreArgs=list(sep="."))
- ## pick out random effects values that correspond to
- ## random effects incorporated in REform ...
- ## FIXME: more tests for possible things going wrong here?
- m <- match(z_names,re_names)
- if (any(is.na(m)))
- stop("random effects specified in REform that were not present in original model")
- re_new <- unlist(re_List[m])
- pred <- pred + drop(as.matrix(re_new %*% ReTrms$Zt))
- } ## REform provided
- } ## predictions with new data or new REform
- ## FIXME: would like to have an isGLMM() accessor for merMod objects?
- if (is(object at resp,"glmResp") && type=="response") {
- pred <- object at resp$family$linkinv(pred)
- }
- return(pred)
+ return(pred)
}
Modified: pkg/lme4Eigen/R/utilities.R
===================================================================
--- pkg/lme4Eigen/R/utilities.R 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/R/utilities.R 2012-02-09 18:39:11 UTC (rev 1584)
@@ -189,7 +189,8 @@
##' From the right hand side of a formula for a mixed-effects model,
##' determine the pairs of expressions that are separated by the
-##' vertical bar operator.
+##' vertical bar operator. Also expand the slash operator in grouping
+##' factor expressions.
##'
##' @title Determine random-effects expressions from a formula
##' @seealso \code{\link{formula}}, \code{\link{model.frame}}, \code{\link{model.matrix}}.
@@ -203,6 +204,8 @@
##' ## => list( Days | Subject )
##' findbars(y ~ Days + (1|Subject) + (0+Days|Subject))
##' ## => list of length 2: list ( 1 | Subject , 0+Days|Subject)
+##' findbars(~ 1 + (1|batch/cask))
+##' ## => list of length 2: list ( 1 | cask:batch , 1 | batch)
##' \dontshow{
##' stopifnot(identical(findbars(f1),
##' list(expression(Days | Subject)[[1]])))
Modified: pkg/lme4Eigen/man/Nelder_Mead.Rd
===================================================================
--- pkg/lme4Eigen/man/Nelder_Mead.Rd 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/man/Nelder_Mead.Rd 2012-02-09 18:39:11 UTC (rev 1584)
@@ -39,9 +39,9 @@
}
\value{
a list with 4 components \item{fval}{numeric scalar - the
- minimum function value achieved} \item{pars}{numeric
+ minimum function value achieved} \item{par}{numeric
vector - the value of \code{x} providing the minimum}
- \item{code}{integer scalar - convergence code}
+ \item{ierr}{integer scalar - error code}
\item{control}{list - the list of control settings after
substituting for defaults}
}
Modified: pkg/lme4Eigen/man/VarCorr.Rd
===================================================================
--- pkg/lme4Eigen/man/VarCorr.Rd 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/man/VarCorr.Rd 2012-02-09 18:39:11 UTC (rev 1584)
@@ -18,17 +18,14 @@
Default is \code{3}.}
}
\value{
- a matrix with the estimated variances, standard
- deviations, and correlations for the random effects. The
- first two columns, named \code{Variance} and
- \code{StdDev}, give, respectively, the variance and the
- standard deviations. If there are correlation components
- in the random effects model, the third column, named
- \code{Corr}, and the remaining unnamed columns give the
- estimated correlations among random effects within the
- same level of grouping. The within-group error variance
- and standard deviation are included as the last row in
- the matrix.
+ a list of matrices, one for each random effects grouping
+ term. For each grouping term, the standard deviations and
+ correlation matrices for each grouping term are stored as
+ attributes \code{"stddev"} and \code{"correlation"},
+ respectively, of the variance-covariance matrix, and the
+ residual standard deviation is stored as attribute
+ \code{"sc"} (for \code{glmer} fits, this attribute stores
+ the scale parameter of the model).
}
\description{
This function calculates the estimated variances,
Modified: pkg/lme4Eigen/man/cbpp.Rd
===================================================================
--- pkg/lme4Eigen/man/cbpp.Rd 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/man/cbpp.Rd 2012-02-09 18:39:11 UTC (rev 1584)
@@ -40,18 +40,14 @@
## response as a matrix
(m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
cbpp, binomial, nAGQ=25L))
-dput(unname(fixef(m1)))
-dput(unname(ranef(m1, drop=TRUE)[[1]]))
## response as a vector of probabilities and usage of argument "weights"
m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
cbpp, binomial, nAGQ=25L)
-dput(unname(fixef(m1p)))
-dput(unname(ranef(m1p, drop=TRUE)[[1]]))
## Confirm that these are equivalent:
stopifnot(all.equal(fixef(m1), fixef(m1p), tol = 1e-5),
all.equal(ranef(m1), ranef(m1p), tol = 1e-5),
TRUE)
-
+## Can this section be moved to a test file? I don't think it belongs in an example. DB
for(m in c(m1, m1p)) {
cat("-------\\n\\nCall: ",
paste(format(getCall(m)), collapse="\\n"), "\\n")
Modified: pkg/lme4Eigen/man/findbars.Rd
===================================================================
--- pkg/lme4Eigen/man/findbars.Rd 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/man/findbars.Rd 2012-02-09 18:39:11 UTC (rev 1584)
@@ -13,7 +13,8 @@
\description{
From the right hand side of a formula for a mixed-effects
model, determine the pairs of expressions that are
- separated by the vertical bar operator.
+ separated by the vertical bar operator. Also expand the
+ slash operator in grouping factor expressions.
}
\section{Note}{
This function is called recursively on individual terms
@@ -26,6 +27,8 @@
## => list( Days | Subject )
findbars(y ~ Days + (1|Subject) + (0+Days|Subject))
## => list of length 2: list ( 1 | Subject , 0+Days|Subject)
+findbars(~ 1 + (1|batch/cask))
+## => list of length 2: list ( 1 | cask:batch , 1 | batch)
\dontshow{
stopifnot(identical(findbars(f1),
list(expression(Days | Subject)[[1]])))
Modified: pkg/lme4Eigen/man/getME.Rd
===================================================================
--- pkg/lme4Eigen/man/getME.Rd 2012-02-08 22:40:20 UTC (rev 1583)
+++ pkg/lme4Eigen/man/getME.Rd 2012-02-09 18:39:11 UTC (rev 1584)
@@ -5,7 +5,7 @@
\title{Extract or Get Generalize Components from a Fitted Mixed Effects Model}
\usage{
getME(object,
- name = c("X", "Z", "Zt", "u", "Gp", "L", "Lambda", "Lambdat", "Lind", "RX", "RZX", "beta", "theta", "REML", "n_rtrms", "is_REML", "devcomp"))
+ name = c("X", "Z", "Zt", "u", "Gp", "L", "Lambda", "Lambdat", "Lind", "RX", "RZX", "beta", "theta", "REML", "n_rtrms", "is_REML", "devcomp", "offset"))
}
\arguments{
\item{object}{a fitted mixed-effects model of class
@@ -36,7 +36,8 @@
\item{is_REML}{same as the result of
\code{\link{isREML}}} \item{devcomp}{a list consisting of
a named numeric vector, \dQuote{cmp}, and a named integer
- vector, \dQuote{dims}, describing the fitted model} }}
+ vector, \dQuote{dims}, describing the fitted model}
+ \item{offset}{model offset} }}
}
\value{
Unspecified, as very much depending on the
More information about the Lme4-commits
mailing list