[Lme4-commits] r1413 - pkg/lme4Eigen/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 1 01:09:25 CEST 2011
Author: dmbates
Date: 2011-10-01 01:09:24 +0200 (Sat, 01 Oct 2011)
New Revision: 1413
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
Log:
Use inheritance on reference classes, being careful about package installation. Create a cleaner version of mkRespMod2.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-09-30 23:05:24 UTC (rev 1412)
+++ pkg/lme4Eigen/R/AllClass.R 2011-09-30 23:09:24 UTC (rev 1413)
@@ -1,4 +1,4 @@
-## Class definitions for the package
+### Class definitions for the package
setClass("lmList",
representation(call = "call",
@@ -9,8 +9,8 @@
setClass("lmList.confint", contains = "array")
### FIXME
-## shouldn't we have "merPred" with two *sub* classes "merPredD" and "merPredS"
-## for the dense and sparse X cases ?
+### shouldn't we have "merPred" with two *sub* classes "merPredD" and "merPredS"
+### for the dense and sparse X cases ?
merPredD <-
setRefClass("merPredD", # Predictor class for mixed-effects models with dense X
@@ -165,21 +165,34 @@
fields =
list(y = "numeric",
Ptr = "externalptr",
- ptr = function (value) {
+ ptr = function (value) { # generates external pointer if necessary
if (!missing(value)) stop("ptr field cannot be assigned")
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(y)) stop("field y must have positive length")
- Ptr <<- .Call(lm_Create, y)
- }
+ try(
+ {
+ if (.Call(isNullExtPtr, Ptr)) {
+ Ptr <<- .Call(lm_Create, y)
+ if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
+ if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
+ }
+ }, silent=TRUE)
Ptr
},
- offset = function(value)
- if (missing(value)) .Call(lm_offset, ptr) else
- .Call(lm_setOffset, ptr, as.numeric(value)),
-
- weights = function(value)
- if (missing(value)) .Call(lm_weights, ptr) else
- .Call(lm_setWeights, ptr, as.numeric(value))),
+ Offset = "numeric",
+ offset = function(value) {
+ if (missing(value))
+ return(tryCatch(.Call(lm_offset, ptr), error=function(e)Offset))
+ ## lm_setOffset checks value's length, etc.
+ .Call(lm_setOffset, ptr, value <- as.numeric(value))
+ Offset <<- value
+ },
+ Weights = "numeric",
+ weights = function(value) {
+ if (missing(value))
+ return(tryCatch(.Call(lm_weights, ptr), error=function(e)Weights))
+ ## lm_setWeights checks value's length, non-negativity, etc.
+ .Call(lm_setWeights, ptr, value <- as.numeric(value))
+ Weights <<- value
+ }),
methods =
list(
fitted = function() {
@@ -195,7 +208,7 @@
.Call(lm_sqrtrwt, ptr)
},
updateMu = function(gamma) {
- 'from the linear predictor, gamma, update the\nconditional mean response, mu, residuals and weights'
+ 'update mu, wtres and wrss from the linear predictor'
.Call(lm_updateMu, ptr, as.numeric(gamma))
},
wrss = function() {
@@ -213,20 +226,35 @@
lmerResp <-
setRefClass("lmerResp",
fields=
- list(REML=function(value)
- if (missing(value)) .Call(lmer_REML, ptr) else .Call(lmer_setREML, ptr, value),
+ list(
ptr = function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(y)) stop("field y must have positive length")
- Ptr <<- .Call(lmer_Create, y)
- }
+ try(
+ if (.Call(isNullExtPtr, Ptr)) {
+ if (!length(y)) stop("field y must have positive length")
+ Ptr <<- .Call(lmer_Create, y)
+ if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
+ if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
+ }, silent=TRUE)
Ptr
- }
- ),
+ },
+ REML="integer",
+ reml=function(value) {
+ if (missing(value))
+ return(tryCatch(.Call(lmer_REML, ptr), error=function(e)0L))
+ stopifnot(length(value <- as.integer(value)) > 0L,
+ (value <- value[1]) >= 0L)
+### FIXME: Change the cls_set<whatever> C functions to return the
+### value that was set so these can be collapsed to one-liners like
+### REML <<- .Call(lmer_setREML, ptr, value) ?
+ .Call(lmer_setREML, ptr, value)
+ REML <<- value
+ }),
contains="lmResp",
methods=
list(
+### FIXME: Should this method be called "deviance" or "objective" or ...?
+### The name Laplace only makes sense for glmer or nlmer.
Laplace = function(ldL2, ldRX2, sqrL) {
'returns the profiled deviance or REML criterion'
.Call(lmer_Laplace, ptr, ldL2, ldRX2, sqrL)
@@ -235,25 +263,32 @@
setOldClass("family")
-
glmResp <-
setRefClass("glmResp",
fields=
- list(family="family",
+ list(
+ family="family",
ptr = function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(y)) stop("field y must have positive length")
- if (!length(family$family)) stop("field family must be initialized")
- Ptr <<- .Call(glm_Create, family, y)
- }
+ try (
+ if (.Call(isNullExtPtr, Ptr)) {
+ if (!length(y))
+ stop("field y must have positive length")
+ if (!length(family$family))
+ stop("field family must be initialized")
+ Ptr <<- .Call(glm_Create, family, y)
+ if (length(Offset)) .Call(lm_setOffset, Ptr, Offset)
+ if (length(Weights)) .Call(lm_setWeights, Ptr, Weights)
+ }, silent=TRUE)
Ptr
- }
- ##,
- ## n=function(value) # only used in aic function for binomial family
- ## if (missing(value)) .Call(glm_n, ptr) else
- ## .Call(glm_setN, ptr, as.numeric(value))
- ),
+ },
+ N="integer",
+ n=function(value) { # only used in aic function for binomial family
+ if (missing(value))
+ return(tryCatch(.Call(glm_n, ptr),
+ error=function(e)rep.int(1, length(y))))
+ .Call(glm_setN, ptr, as.numeric(value))
+ }),
contains="lmResp",
methods=
list(
@@ -298,7 +333,7 @@
.Call(glm_sqrtWrkWt, ptr)
},
updateMu = function(gamma) {
- 'from the linear predictor, gamma, update the\nconditional mean response, mu, residuals and weights'
+ 'update mu, residuals, weights, etc. from the linear predictor'
.Call(glm_updateMu, ptr, as.numeric(gamma))
},
updateWts = function() {
@@ -331,10 +366,9 @@
pnames="character",
ptr = function (value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(y)) stop("field y must have positive length")
- Ptr <<- .Call(nls_Create, nlmod, nlenv, y)
- }
+ try (if (.Call(isNullExtPtr, Ptr))
+ Ptr <<- .Call(nls_Create, nlmod, nlenv, y),
+ silent=TRUE)
Ptr
}),
contains="lmResp",
@@ -345,7 +379,7 @@
.Call(nls_Laplace, ptr, ldL2, ldRX2, sqrL)
},
updateMu=function(gamma) {
- 'from the linear predictor, gamma, update the\nconditional mean response, mu, residuals and weights'
+ 'update mu, residuals, gradient, etc. given the linear predictor matrix'
.Call(nls_updateMu, ptr, as.numeric(gamma))
})
)
@@ -358,10 +392,11 @@
family="family",
ptr=function(value) {
if (!missing(value)) stop("ptr field cannot be assigned")
- if (.Call(isNullExtPtr, Ptr)) {
- if (!length(family$family)) stop("field family must be initialized")
- Ptr <<- .Call(glmFamily_Create, family)
- }
+ try(
+ if (.Call(isNullExtPtr, Ptr)) {
+ if (!length(family$family)) stop("field family must be initialized")
+ Ptr <<- .Call(glmFamily_Create, family)
+ }, silent=TRUE)
Ptr
}),
methods=
@@ -395,7 +430,8 @@
)
setClass("merMod",
- representation(call = "call",
+ representation(Gp = "integer",
+ call = "call",
frame = "data.frame", # "model.frame" is not S4-ized yet
flist = "list",
cnms = "list",
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2011-09-30 23:05:24 UTC (rev 1412)
+++ pkg/lme4Eigen/R/lmer.R 2011-09-30 23:09:24 UTC (rev 1413)
@@ -242,7 +242,6 @@
##' @param nlmod the nonlinear model function (nlsRespMod only)
##' @return a lmerResp (or glmerResp) instance
mkRespMod2 <- function(fr, family = NULL, nlenv = NULL, nlmod = NULL) {
- n <- nrow(fr)
y <- model.response(fr)
if(length(dim(y)) == 1) {
## avoid problems with 1D arrays, but keep names
@@ -250,50 +249,47 @@
dim(y) <- NULL
if(!is.null(nm)) names(y) <- nm
}
- if (is.null(nlmod)) {
- if (is.null(family)) {
- ans <- new("lmerResp", y)
- } else {
- stopifnot(inherits(family, "family"))
- rho <- new.env()
- rho$etastart <- model.extract(fr, "etastart")
- rho$mustart <- model.extract(fr, "mustart")
- rho$nobs <- n
- rho$y <- y
- eval(family$initialize, rho)
- family$initialize <- NULL # remove clutter from str output
- ans <- new("glmerResp", family, rho$y)
- ans$updateMu(family$linkfun(unname(rho$mustart)))
- }
- if (!is.null(offset <- model.offset(fr))) {
- if (length(offset) == 1L) offset <- rep.int(offset, n)
- stopifnot(length(offset) == n)
- ans$offset <- unname(offset)
- }
- if (!is.null(weights <- model.weights(fr))) {
- stopifnot(length(weights) == n, all(weights >= 0))
- ans$weights <- unname(weights)
- }
- return(ans)
+ rho <- new.env()
+ rho$y <- y
+ n <- N <- nrow(fr)
+ if (!is.null(nlenv)) {
+ stopifnot(is.numeric(val <- eval(nlmod, nlenv)),
+ length(val) == N,
+ is.matrix(gr <- attr(val, "gradient")),
+ mode(gr) == "numeric",
+ nrow(gr) == N,
+ !is.null(pnames <- colnames(gr)))
+ N <- length(gr)
}
- stopifnot(is.language(nlmod), is.environment(nlenv))
- val <- eval(nlmod, env=nlenv)
- stopifnot(is.numeric(val),
- length(val) == length(y),
- is.matrix(gr <- attr(val, "gradient")),
- mode(gr) == "numeric",
- nrow(gr) == length(y))
- ans <- new("nlsResp", nlenv, nlmod, colnames(gr), y)
if (!is.null(offset <- model.offset(fr))) {
- if (length(offset) == 1L) offset <- rep.int(offset, length(gr))
- stopifnot(length(offset) == length(gr))
- ans$offset <- unname(offset)
+ if (length(offset) == 1L) offset <- rep.int(offset, N)
+ stopifnot(length(offset) == N)
+ rho$offset <- unname(offset)
}
if (!is.null(weights <- model.weights(fr))) {
stopifnot(length(weights) == n, all(weights >= 0))
- ans$weights <- unname(weights)
+ rho$weights <- unname(weights)
}
- ans
+ if (is.null(family) && is.null(nlenv))
+ return(do.call(new, c(list(Class="lmerResp"), as.list(rho))))
+ if (!is.null(family)) {
+ stopifnot(inherits(family, "family"))
+ # need weights for initialize evaluation
+ if (!exists("weights", rho)) rho$weights <- rep.int(1, n)
+ rho$nobs <- n
+ eval(family$initialize, rho)
+ family$initialize <- NULL # remove clutter from str output
+ ll <- as.list(rho)
+ ans <- do.call(new, c(list(Class="glmResp", family=family),
+ ll[setdiff(names(ll),
+ c("m", "nobs", "mustart"))]))
+ ans$updateMu(if (!is.null(es <- model.extract(fr, "etastart"))) es else
+ family$linkfun(get("mustart", rho)))
+ return(ans)
+ }
+ stopifnot(is.language(nlmod), is.environment(nlenv))
+ do.call(new, c(list(Class="nlsResp", nlenv=nlenv, nlmod=nlmod,
+ pnames=pnames), as.list(rho)))
}
###' Create Z, Lambda, Lind, etc.
More information about the Lme4-commits
mailing list