[Lme4-commits] r1487 - in pkg/lme4Eigen: R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 22 20:51:52 CET 2011
Author: dmbates
Date: 2011-12-22 20:51:52 +0100 (Thu, 22 Dec 2011)
New Revision: 1487
Modified:
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/man/lmer.Rd
pkg/lme4Eigen/src/external.cpp
Log:
Incorporated aGQ code in the glmer function itself. Eliminated the doFit argument which is redundant when devFunOnly is also used.
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-12-21 23:00:19 UTC (rev 1486)
+++ pkg/lme4Eigen/R/lmer.R 2011-12-22 19:51:52 UTC (rev 1487)
@@ -1,9 +1,8 @@
### FIXME: Should there be both a doFit and a devFunOnly argument?
lmer <- function(formula, data, REML = TRUE, sparseX = FALSE,
control = list(), start = NULL,
- verbose = 0L, doFit = TRUE,
- subset, weights, na.action, offset,
- contrasts = NULL, devFunOnly=0L, ...)
+ verbose = 0L, subset, weights, na.action, offset,
+ contrasts = NULL, devFunOnly=FALSE, ...)
{
if (sparseX) warning("sparseX = TRUE has no effect at present")
mf <- mc <- match.call()
@@ -54,14 +53,11 @@
rho$resp <- mkRespMod2(fr, if(REML) p else 0L)
devfun <- mkdevfun(rho, 0L)
+ devfun(reTrms$theta) # one evaluation to ensure all values are set
+
if (devFunOnly) return(devfun)
- # one evaluation to ensure all values are set
- opt <- list(fval=devfun(reTrms$theta))
-
- if (doFit) { # optimize estimates
- if (verbose) control$iprint <- as.integer(verbose)
- opt <- bobyqa(reTrms$theta, devfun, reTrms$lower, control = control)
- }
+ if (verbose) control$iprint <- as.integer(verbose)
+ opt <- bobyqa(reTrms$theta, devfun, reTrms$lower, control = control)
sqrLenU <- rho$pp$sqrL(1.)
wrss <- rho$resp$wrss()
pwrss <- wrss + sqrLenU
@@ -86,14 +82,22 @@
if (is(rho$resp, "lmerResp"))
ff <- function(theta) .Call(lmer_Deviance, pp$ptr(), resp$ptr(), theta)
else if (is(rho$resp, "glmResp")) {
- rho$glmerLaplace <- glmerLaplace
- ff <- switch(nAGQ + 1L,
- function(theta)
- .Call(glmerLaplace, pp$ptr(), resp$ptr(), theta,
- u0, beta0, verbose, FALSE, tolPwrss),
- function(pars)
- .Call(glmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], u0,
- pars[-dpars], verbose, TRUE, tolPwrss))
+ if (nAGQ < 2L) {
+ rho$glmerLaplace <- glmerLaplace
+ ff <- switch(nAGQ + 1L,
+ function(theta)
+ .Call(glmerLaplace, pp$ptr(), resp$ptr(), theta,
+ u0, beta0, verbose, FALSE, tolPwrss),
+ function(pars)
+ .Call(glmerLaplace, pp$ptr(), resp$ptr(), pars[dpars], u0,
+ pars[-dpars], verbose, TRUE, tolPwrss))
+ } else {
+ rho$glmerAGQ <- glmerAGQ
+ rho$GQmat <- GHrule(nAGQ)
+ ff <- function(pars)
+ .Call(glmerAGQ, pp$ptr(), resp$ptr(), fac, GQmat, pars[dpars],
+ u0, pars[-dpars], tolPwrss)
+ }
}
if (is.null(ff)) stop("code not yet written")
environment(ff) <- rho
@@ -102,8 +106,8 @@
glmer <- function(formula, data, family = gaussian, sparseX = FALSE,
control = list(), start = NULL, verbose = 0L, nAGQ = 1L,
- doFit = TRUE, compDev=TRUE, subset, weights, na.action, offset,
- contrasts = NULL, mustart, etastart, devFunOnly = 0L,
+ compDev=TRUE, subset, weights, na.action, offset,
+ contrasts = NULL, mustart, etastart, devFunOnly = FALSE,
tolPwrss = 1e-10, ...)
{
verbose <- as.integer(verbose)
@@ -133,9 +137,12 @@
if(is.function(family)) family <- family()
if (family$family %in% c("quasibinomial", "quasipoisson", "quasi"))
stop('"quasi" families cannot be used in glmer')
- nAGQ <- as.integer(nAGQ)[1]
- if (nAGQ > 1L) warning("nAGQ > 1 has not been implemented, using Laplace")
- stopifnot(length(formula <- as.formula(formula)) == 3)
+
+ if (nAGQ < 0L || nAGQ > 25L) warning("nAGQ must be in the range [0, 25]")
+ stopifnot(length(nAGQ <- as.integer(nAGQ)) == 1L,
+ nAGQ >= 0L,
+ nAGQ <= 25L,
+ length(formula <- as.formula(formula)) == 3)
if (missing(data)) data <- environment(formula)
# evaluate and install the model frame :
m <- match(c("data", "subset", "weights", "na.action", "offset",
@@ -176,27 +183,29 @@
rho$u0 <- rho$pp$u0
rho$beta0 <- rho$pp$beta0
+ rho$lower <- reTrms$lower
- opt <- list(fval=rho$resp$Laplace(rho$pp$ldL2(), rho$pp$ldRX2(), rho$pp$sqrL(0.)))
-
- if (doFit || devFunOnly) { # create the deviance function
- parent.env(rho) <- parent.frame()
- devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
- if (devFunOnly && !nAGQ) return(devfun)
- control$iprint <- min(verbose, 3L)
- opt <- bobyqa(rho$pp$theta, devfun, lower=reTrms$lower, control=control)
- if (nAGQ > 0L) {
- rho$u0 <- rho$pp$u0
- rho$dpars <- seq_along(rho$pp$theta)
- if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
- if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
- rho$control <- control
- devfun <- mkdevfun(rho, nAGQ)
- if (devFunOnly) return(devfun)
- opt <- bobyqa(c(rho$pp$theta, rho$pp$beta0), devfun,
- lower=c(reTrms$lower, rep.int(-Inf, length(rho$pp$beta0))),
- control=control)
+ parent.env(rho) <- parent.frame()
+ devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
+ if (devFunOnly && !nAGQ) return(devfun)
+ control$iprint <- min(verbose, 3L)
+ opt <- bobyqa(rho$pp$theta, devfun, rho$lower, control=control)
+ if (nAGQ > 0L) {
+ rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
+ rho$u0 <- rho$pp$u0
+ rho$beta0 <- rho$pp$beta0
+ rho$dpars <- seq_along(rho$pp$theta)
+ if (nAGQ > 1L) {
+ if (length(reTrms$flist) != 1L || length(reTrms$cnms[[1]]) != 1L)
+ stop("nAGQ > 1 is only available for models with a single, scalar random-effects term")
+ rho$fac <- reTrms$flist[[1]]
}
+ if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
+ if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
+ rho$control <- control
+ devfun <- mkdevfun(rho, nAGQ)
+ if (devFunOnly) return(devfun)
+ opt <- bobyqa(c(rho$pp$theta, rho$beta0), devfun, rho$lower, control=control)
}
sqrLenU <- rho$pp$sqrL(0.)
wrss <- rho$resp$wrss()
Modified: pkg/lme4Eigen/man/lmer.Rd
===================================================================
--- pkg/lme4Eigen/man/lmer.Rd 2011-12-21 23:00:19 UTC (rev 1486)
+++ pkg/lme4Eigen/man/lmer.Rd 2011-12-22 19:51:52 UTC (rev 1487)
@@ -13,16 +13,15 @@
\usage{
lmer(formula, data, REML = TRUE, sparseX = FALSE,
control = list(), start = NULL, verbose = 0L,
- doFit = TRUE, %compDev = TRUE, _support too - FIXME? _
- subset, weights,
- na.action, offset, contrasts = NULL,
- devFunOnly=0L,
+ %doFit = TRUE, %compDev = TRUE, _support too - FIXME? _
+ subset, weights, na.action, offset, contrasts = NULL,
+ devFunOnly = FALSE,
\dots)
glmer(formula, data, family = gaussian, sparseX = FALSE,
control = list(), start = NULL, verbose = 0, nAGQ = 1,
- doFit = TRUE, compDev = TRUE, subset, weights, na.action, offset,
- contrasts = NULL, mustart, etastart, devFunOnly = 0L,
+ compDev = TRUE, subset, weights, na.action, offset,
+ contrasts = NULL, mustart, etastart, devFunOnly = FALSE,
tolPwrss = 1e-06, \dots)
nlmer(formula, data, family = gaussian, start = NULL,
@@ -66,11 +65,6 @@
starting value of \code{theta}. In \code{nlmer} a numeric
\code{start} argument is used as the starting values of the
\code{fixef} slot.}
- \item{doFit}{logical scalar. When \code{doFit = FALSE} the model is
- not fit but instead a structure with the model matrices for the
- random-effects terms is returned, so they can be modified for
- special model forms. When \code{doFit = TRUE}, the default, the
- model is fit immediately.}
\item{compDev}{logical - should compiled code be used for the deviance
evaluation during the optimization of the parameters estimates?}
\item{subset, weights, na.action, offset, contrasts}{further model
@@ -84,7 +78,7 @@
\code{> 1} verbose output is generated during the individual PIRLS
steps.}
\item{devFunOnly}{logical - return only the deviance evaluation
- function - useful mainly for debugging.}% but MM thinks we want to keep.
+ function - useful mainly for debugging.}% but MM thinks we want to keep; DB agrees
\item{tolPwrss}{numeric scalar - the tolerance for declaring
convergence in the penalized weighted residual sum-of-squares
step. Defaults to 1e-06.}
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2011-12-21 23:00:19 UTC (rev 1486)
+++ pkg/lme4Eigen/src/external.cpp 2011-12-22 19:51:52 UTC (rev 1487)
@@ -222,7 +222,8 @@
static double sqrt2pi = std::sqrt(2. * PI);
- SEXP glmerAGQ(SEXP pp_, SEXP rp_, SEXP fac_, SEXP GQmat_, SEXP tol_) {
+ SEXP glmerAGQ(SEXP pp_, SEXP rp_, SEXP fac_, SEXP GQmat_, SEXP theta_,
+ SEXP u0_, SEXP beta0_, SEXP tol_) {
BEGIN_RCPP;
XPtr<glmResp> rp(rp_);
@@ -231,6 +232,9 @@
if (fac.size() != rp->mu().size())
throw std::invalid_argument("size of fac must match dimension of response vector");
+ pp->setTheta(as<MVec>(theta_));
+ pp->setU0(as<MVec>(u0_));
+ pp->setBeta0(as<MVec>(beta0_));
pwrssUpdate(rp, pp, 0, true, ::Rf_asReal(tol_)); // should be a no-op
const ArrayXd u0(pp->u0());
const ArrayXd devc0(devcCol(fac, u0, rp->devResid()));
@@ -690,7 +694,7 @@
CALLDEF(glmFamily_muEta, 2),
CALLDEF(glmFamily_variance, 2),
- CALLDEF(glmerAGQ, 5),
+ CALLDEF(glmerAGQ, 8),
CALLDEF(glmerPwrssUpdate, 5),
CALLDEF(glmerWrkIter, 2),
CALLDEF(glmerLaplace, 8),
More information about the Lme4-commits
mailing list