[Lme4-commits] r1531 - in branches/roxygen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jan 29 00:25:35 CET 2012
Author: dmbates
Date: 2012-01-29 00:25:35 +0100 (Sun, 29 Jan 2012)
New Revision: 1531
Modified:
branches/roxygen/R/AllClass.R
branches/roxygen/R/AllGeneric.R
branches/roxygen/R/lmer.R
branches/roxygen/src/respModule.cpp
branches/roxygen/src/respModule.h
Log:
Add a refit generic and a method for merMod, add a setResp method for the lmResp class.
Modified: branches/roxygen/R/AllClass.R
===================================================================
--- branches/roxygen/R/AllClass.R 2012-01-28 18:57:02 UTC (rev 1530)
+++ branches/roxygen/R/AllClass.R 2012-01-28 23:25:35 UTC (rev 1531)
@@ -7,10 +7,7 @@
##' Class \code{"lmList"} is an S4 class with basically a list of objects of
##' class \code{\link{lm}} with a common model.
##' @name lmList-class
-##' @aliases lmList-class lmList.confint-class coef,lmList-method
-##' formula,lmList-method confint,lmList-method plot,lmList-method
-##' show,lmList-method update,lmList-method plot,lmList.confint-method
-##' plot,lmList.confint,ANY-method
+##' @aliases lmList-class show,lmList-method
##' @docType class
##' @section Objects from the Class: Objects can be created by calls of the form
##' \code{new("lmList", ...)} or, more commonly, by a call to
@@ -77,7 +74,7 @@
methods =
list(
### FIXME: probably don't need S as Xwts is required for nlmer
- initialize = function(X, Zt, Lambdat, Lind, theta, S, ...) {
+ initialize = function(X, Zt, Lambdat, Lind, theta, n, ...) {
if (!nargs()) return
X <<- as(X, "matrix")
Zt <<- as(Zt, "dgCMatrix")
@@ -85,15 +82,11 @@
Lind <<- as.integer(Lind)
theta <<- as.numeric(theta)
N <- nrow(X)
- S <- as.integer(S)[1]
- n <- N %/% S
p <- ncol(X)
q <- nrow(Zt)
stopifnot(length(theta) > 0L,
length(Lind) > 0L,
- all(sort(unique(Lind)) == seq_along(theta)),
- S > 0L,
- n * S == N)
+ all(sort(unique(Lind)) == seq_along(theta)))
RZX <<- array(0, c(q, p))
Utr <<- numeric(q)
V <<- array(0, c(n, p))
@@ -105,7 +98,7 @@
delu <<- numeric(q)
uu <- list(...)$u0
u0 <<- if (is.null(uu)) numeric(q) else uu
- Ut <<- if (S == 1) Zt + 0 else
+ Ut <<- if (n == N) Zt + 0 else
Zt %*% sparseMatrix(i=seq_len(N), j=as.integer(gl(n, 1, N)), x=rep.int(1,N))
LamtUt <<- Lambdat %*% Ut
Xw <- list(...)$Xwts
@@ -158,7 +151,7 @@
assign(field, current, envir = vEnv)
}
}
- do.call(new, c(as.list(vEnv), Class=def))
+ do.call(new, c(as.list(vEnv), n=nrow(vEnv$V), Class=def))
},
ldL2 = function() {
'twice the log determinant of the sparse Cholesky factor'
@@ -362,10 +355,18 @@
}
Ptr
},
- setOffset = function(oo) {
+ setOffset = function(oo) {
'change the offset in the model (used in profiling)'
.Call(lm_setOffset, ptr(), as.numeric(oo))
},
+ setResp = function(rr) {
+ 'change the response in the model, usually after a deep copy'
+ .Call(lm_setResp, ptr(), as.numeric(rr))
+ },
+ setWeights = function(ww) {
+ 'change the prior weights in the model'
+ .Call(lm_setWeights, ptr(), as.numeric(ww))
+ },
updateMu = function(gamma) {
'update mu, wtres and wrss from the linear predictor'
.Call(lm_updateMu, ptr(), as.numeric(gamma))
Modified: branches/roxygen/R/AllGeneric.R
===================================================================
--- branches/roxygen/R/AllGeneric.R 2012-01-28 18:57:02 UTC (rev 1530)
+++ branches/roxygen/R/AllGeneric.R 2012-01-28 23:25:35 UTC (rev 1531)
@@ -1,4 +1,3 @@
-
## utilities, these *exported*:
##' @export getL
setGeneric("getL", function(x) standardGeneric("getL"))
@@ -53,6 +52,22 @@
##' @export
refitML <- function(x, ...) UseMethod("refitML")
+##' Refit a model with a different response vector
+##'
+##' Refit a model after modifying the response vector. This could be done using
+##' an \code{\link{update}} method but this approach should be faster because
+##' it bypasses the creation of the model representation and goes directly to
+##' the optimization step.
+##' @title Refit a model by maximum likelihood criterion
+##' @param object a fitted model, usually of class \code{"\linkS4class{lmerMod}"},
+##' to be refit with a new response
+##' @param newresp a numeric vector providing the new response. Must be
+##' 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
+##' @export
+refit <- function(object, newresp, ...) UseMethod("refit")
+
if (FALSE) {
setGeneric("HPDinterval",
function(object, prob = 0.95, ...) standardGeneric("HPDinterval"))
Modified: branches/roxygen/R/lmer.R
===================================================================
--- branches/roxygen/R/lmer.R 2012-01-28 18:57:02 UTC (rev 1530)
+++ branches/roxygen/R/lmer.R 2012-01-28 23:25:35 UTC (rev 1531)
@@ -133,7 +133,7 @@
if ((qrX <- qr(X))$rank < p)
stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
rho <- new.env(parent=parent.env(environment()))
- rho$pp <- do.call(merPredD$new, c(list(X=X, S=1L), reTrms[c("Zt","theta","Lambdat","Lind")]))
+ rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")], n=nrow(X), list(X=X)))
rho$resp <- mkRespMod(fr, if(REML) p else 0L)
devfun <- mkdevfun(rho, 0L)
@@ -178,7 +178,7 @@
cmp <- c(ldL2=rho$pp$ldL2(), ldRX2=rho$pp$ldRX2(), wrss=wrss,
ussq=sqrLenU, pwrss=pwrss,
drsum=NA, dev=if(REML)NA else opt$fval, REML=if(REML)opt$fval else NA,
- sigmaML=sqrt(pwrss/n), sigmaREML=sqrt(pwrss/(n-p)))
+ sigmaML=sqrt(pwrss/n), sigmaREML=sqrt(pwrss/(n-p)), tolPwrss=NA)
new("lmerMod", call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms,
theta=rho$pp$theta, beta=rho$pp$delb, u=rho$pp$delu, lower=reTrms$lower,
@@ -334,7 +334,7 @@
p <- ncol(X)
rho <- as.environment(list(verbose=verbose, tolPwrss=tolPwrss))
parent.env(rho) <- parent.frame()
- rho$pp <- do.call(merPredD$new, c(list(X=X, S=1L), reTrms[c("Zt","theta","Lambdat","Lind")]))
+ rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")], n=nrow(X), list(X=X)))
# response module
rho$resp <- mkRespMod(fr, family=family)
# initial step from working response
@@ -418,7 +418,7 @@
cmp <- c(ldL2=rho$pp$ldL2(), ldRX2=rho$pp$ldRX2(), wrss=wrss,
ussq=sqrLenU, pwrss=pwrss,
drsum=rho$resp$resDev(), dev=opt$fval, REML=NA,
- sigmaML=NA, sigmaREML=NA)
+ sigmaML=NA, sigmaREML=NA, tolPwrss=tolPwrss)
new("glmerMod", call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms,
Gp=reTrms$Gp, theta=rho$pp$theta, beta=rho$pp$beta0, u=rho$pp$u0,
@@ -469,7 +469,7 @@
parent=parent.frame())
rho$pp <- do.call(merPredD$new,
c(vals$reTrms[c("Zt","theta","Lambdat","Lind")],
- list(X=X, S=length(vals$pnames), Xwts=vals$respMod$sqrtXwt,
+ list(X=X, n=length(vals$respMod$mu), Xwts=vals$respMod$sqrtXwt,
beta0=qr.coef(qrX, unlist(lapply(vals$pnames, get,
envir = rho$resp$nlenv))))))
rho$u0 <- rho$pp$u0
@@ -570,7 +570,7 @@
cmp <- c(ldL2=rho$pp$ldL2(), ldRX2=rho$pp$ldRX2(), wrss=wrss,
ussq=sqrLenU, pwrss=pwrss, drsum=NA,
drsum=wrss, dev=opt$fval, REML=NA,
- sigmaML=sqrt(pwrss/n), sigmaREML=NA)
+ sigmaML=sqrt(pwrss/n), sigmaREML=NA, tolPwrss=tolPwrss)
new("nlmerMod", call=mc, frame=vals$fr, flist=vals$reTrms$flist, cnms=vals$reTrms$cnms,
Gp=vals$reTrms$Gp, theta=rho$pp$theta, beta=rho$pp$beta0, u=rho$pp$u0,
@@ -1124,7 +1124,7 @@
ranef.merMod <- function(object, postVar = FALSE, drop = FALSE,
whichel = names(ans), ...)
{
- ans <- as.vector(crossprod(object at pp$Lambdat, object at u))
+ ans <- object at pp$b(1.)
if (!is.null(object at flist)) {
## evaluate the list of matrices
levs <- lapply(fl <- object at flist, levels)
@@ -1178,7 +1178,7 @@
rho$resp <- new(class(rr), y=rr$y, offset=rr$offset, weights=rr$weights, REML=0L)
xpp <- x at pp
rho$pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat,
- Lind=xpp$Lind, theta=xpp$theta, S=1L)
+ Lind=xpp$Lind, theta=xpp$theta, n=nrow(xpp$X))
devfun <- mkdevfun(rho, 0L)
opt <- bobyqa(x at theta, devfun, x at lower)
n <- length(rr$y)
@@ -1474,7 +1474,6 @@
stop(sprintf("Mixed-Effects extraction of '%s' is not available for class \"%s\"",
name, class(object))))
}## {getME}
-
##' @importMethodsFrom Matrix t %*% crossprod diag tcrossprod
##' @importClassesFrom Matrix dgCMatrix dpoMatrix corMatrix
NULL
@@ -1904,3 +1903,48 @@
##' @importFrom stats anova
##' @S3method anova merMod
anova.merMod <- anovaLmer
+
+##' @S3method refit merMod
+refit.merMod <- function(object, newresp, ...) {
+ rr <- object at resp$copy()
+ stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == length(rr$y))
+ rr$setResp(newresp)
+ pp <- object at pp$copy()
+ lower <- object at lower
+ tolPwrss <- object at devcomp$cmp["tolPwrss"]
+ ff <- mkdevfun(list2env(pp=pp, resp=rr, u0=pp$u0, beta0=pp$beta0, verbose=0L,
+ tolPwrss=tolPwrss, dpars=seq_along(lower)),
+ nAGQ=object at devcomp$dims["nAGQ"])
+ xst <- c(rep.int(0.1, length(rho$dpars)), sqrt(diag(pp$unsc())))
+ nM <- NelderMead$new(lower=lower, upper=rep.int(Inf, length(lower)), xst=0.2*xst,
+ x0=c(pp$theta, pp$beta0), xt=xst*0.0001)
+### FIXME: Probably should save the control settings in the merMod object
+ nM$setMaxeval(10000L)
+ nM$setFtolAbs(1e-5)
+ nM$setFtolRel(1e-15)
+ nM$setMinfMax(.Machine$double.xmin)
+ nM$setIprint(0L)
+ while ((nMres <- nM$newf(devfun(nM$xeval()))) == 0L) {}
+ if (nMres < 0L) {
+ if (nMres > -4L)
+ stop("convergence failure, code ", nMres, " in NelderMead")
+ else
+ warning("failure to converge in 1000 evaluations")
+ }
+ opt <- list(fval=nM$value(), pars=nM$xpos(), code=nMres)
+### FIXME: Abstract these operations into another function
+ devcomp <- object at devcomp
+ cmp.old <- devcomp$cmp
+ sqrLenU <- rho$pp$sqrL(0.)
+ wrss <- rho$resp$wrss()
+ pwrss <- wrss + sqrLenU
+ cmp <- c(ldL2=rho$pp$ldL2(), ldRX2=rho$pp$ldRX2(), wrss=wrss, ussq=sqrLenU, pwrss=pwrss,
+ dev=opt$fval, REML=NA, sigmaML=NA, sigmaREML=NA, drsum=opt$fval-sqrLenU,
+ tolPwrss=tolPwrss)
+
+ new(class(object), call=object at call, frame=object at frame, flist=object at flist, cnms=object at cnms,
+ Gp=reTrms$Gp, theta=rho$pp$theta, beta=rho$pp$beta0, u=rho$pp$u0,
+ lower=reTrms$lower, devcomp=list(cmp=cmp, dims=dims), pp=rho$pp,
+ resp=rho$resp)
+}
+
Modified: branches/roxygen/src/respModule.cpp
===================================================================
--- branches/roxygen/src/respModule.cpp 2012-01-28 18:57:02 UTC (rev 1530)
+++ branches/roxygen/src/respModule.cpp 2012-01-28 23:25:35 UTC (rev 1531)
@@ -31,6 +31,12 @@
updateWrss();
}
+ /**
+ * Update the (conditional) mean response and weighted residuals
+ *
+ * @param gamma New value of the linear predictor
+ * @return updated sum of squared, weighted residuals
+ */
double lmResp::updateMu(const VectorXd& gamma) {
if (gamma.size() != d_offset.size())
throw invalid_argument("updateMu: Size mismatch");
@@ -51,12 +57,24 @@
return d_wrss;
}
+/**
+ * Set a new value of the offset.
+ *
+ * The values are copied into the d_offset member because it is mapped.
+ * @param oo New value of the offset
+ */
void lmResp::setOffset(const VectorXd& oo) {
if (oo.size() != d_offset.size())
throw invalid_argument("setOffset: Size mismatch");
- d_offset = oo;
+ d_offset = oo; // this copies the values
}
+ void lmResp::setResp(const VectorXd& yy) {
+ if (yy.size() != d_y.size())
+ throw invalid_argument("setResp: Size mismatch");
+ d_y = yy;
+ }
+
void lmResp::setWeights(const VectorXd& ww) {
if (ww.size() != d_weights.size())
throw invalid_argument("setWeights: Size mismatch");
Modified: branches/roxygen/src/respModule.h
===================================================================
--- branches/roxygen/src/respModule.h 2012-01-28 18:57:02 UTC (rev 1530)
+++ branches/roxygen/src/respModule.h 2012-01-28 23:25:35 UTC (rev 1531)
@@ -23,30 +23,57 @@
class lmResp {
protected:
- double d_wrss;
- const MVec d_y;
- MVec d_weights, d_offset, d_mu, d_sqrtXwt, d_sqrtrwt, d_wtres;
+ double d_wrss; /**< current weighted sum of squared residuals */
+ MVec d_y, /**< response vector */
+ d_weights, /**< prior weights - always present even if unity */
+ d_offset, /**< offset in the model */
+ d_mu, /**< mean response from current linear predictor */
+ d_sqrtXwt, /**< Square roots of the "X weights". For
+ * lmResp and lmerResp these are the same as
+ * the sqrtrwt. For glmResp and nlsResp they
+ * incorporate the gradient of the eta to mu
+ * mapping.*/
+ d_sqrtrwt, /**< Square roots of the residual weights */
+ d_wtres; /**< Current weighted residuals */
public:
lmResp(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
const MVec& sqrtXwt() const {return d_sqrtXwt;}
- const MVec& mu() const {return d_mu;}
+ /**< return a const reference to d_sqrtXwt */
+ const MVec& mu() const {return d_mu;}
+ /**< return a const reference to d_mu */
const MVec& offset() const {return d_offset;}
+ /**< return a const reference to d_offset */
const MVec& sqrtrwt() const {return d_sqrtrwt;}
+ /**< return a const reference to d_sqrtrwt */
const MVec& weights() const {return d_weights;}
+ /**< return a const reference to d_weights */
const MVec& wtres() const {return d_wtres;}
+ /**< return a const reference to d_wtres */
const MVec& y() const {return d_y;}
+ /**< return a const reference to d_y */
double wrss() const {return d_wrss;}
+ /**< return the weighted sum of squared residuals */
double updateMu(const Eigen::VectorXd&);
double updateWts() {return updateWrss();}
- double updateWrss();
+ /**< update the weights. For a
+ * glmResp this done separately from
+ * updating the mean, because of the
+ * iterative reweighting. */
+ double updateWrss(); /**< update the weighted residuals and d_wrss */
void setOffset(const Eigen::VectorXd&);
+ /**< set a new value of the offset */
+ void setResp(const Eigen::VectorXd&);
+ /**< set a new value of the response, y */
void setWeights(const Eigen::VectorXd&);
+ /**< set a new value of the prior weights */
+
};
class lmerResp : public lmResp {
private:
- int d_reml;
+ int d_reml; /**< 0 for evaluating the deviance, p
+ * for evaluating the REML criterion. */
public:
lmerResp(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
More information about the Lme4-commits
mailing list