[Lme4-commits] r1603 - in pkg/lme4Eigen: . R inst/tests tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 21 03:16:23 CET 2012
Author: bbolker
Date: 2012-02-21 03:16:21 +0100 (Tue, 21 Feb 2012)
New Revision: 1603
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/inst/tests/test-lmer.R
pkg/lme4Eigen/tests/simulate.R
Log:
work on refit
* add nAGQ to environment for refit
* don't add entries to lower-bound vector unnecessarily
simulate:
* cleanup: allow prob/weight formulation for binomial,
bernoulli response specified as factor
add weights method
extend simulate.R testing
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-02-15 21:54:04 UTC (rev 1602)
+++ pkg/lme4Eigen/DESCRIPTION 2012-02-21 02:16:21 UTC (rev 1603)
@@ -1,11 +1,11 @@
Package: lme4Eigen
-Version: 0.9996875-9
+Version: 0.9996875-10
Date: $Date$
Revision: $Revision$
Title: Linear mixed-effects models using Eigen and S4
Author: Douglas Bates <bates at stat.wisc.edu>,
Martin Maechler <maechler at R-project.org> and
- Ben Bolker <bbolker at gmail.com>
+ Ben Bolker <bolker at mcmaster.ca>
Maintainer: <lme4-authors at R-forge.wu-wien.ac.at>
Description: Fit linear and generalized linear mixed-effects models.
The models and their components are represented using S4 classes and
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-02-15 21:54:04 UTC (rev 1602)
+++ pkg/lme4Eigen/NAMESPACE 2012-02-21 02:16:21 UTC (rev 1603)
@@ -80,12 +80,12 @@
importFrom(stats,terms)
importFrom(stats,update)
importFrom(stats,vcov)
+importFrom(stats,weights)
importMethodsFrom(Matrix,"%*%")
importMethodsFrom(Matrix,coerce)
importMethodsFrom(Matrix,crossprod)
importMethodsFrom(Matrix,diag)
importMethodsFrom(Matrix,t)
-importMethodsFrom(Matrix,tcrossprod)
S3method(anova,merMod)
S3method(as.function,merMod)
S3method(coef,lmList)
@@ -136,5 +136,6 @@
S3method(VarCorr,merMod)
S3method(vcov,merMod)
S3method(vcov,summary.mer)
+S3method(weights,merMod)
S3method(xyplot,thpr)
useDynLib(lme4Eigen,.registration=TRUE)
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-15 21:54:04 UTC (rev 1602)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-21 02:16:21 UTC (rev 1603)
@@ -482,6 +482,7 @@
##' minqa::bobyqa(1, dd, 0)
##'
mkdevfun <- function(rho, nAGQ=1L) {
+ ## FIXME: should nAGQ be automatically embedded in rho?
stopifnot(is.environment(rho), is(rho$resp, "lmResp"))
ff <- NULL
if (is(rho$resp, "lmerResp")) {
@@ -1166,6 +1167,13 @@
rr <- object at resp$copy()
## FIXME: want this to work for binomial (two-column) data?
## or tell people they have to work with prop/weights formulation?
+ if (is.matrix(newresp) && ncol(newresp)==2 && rr$family$family=="binomial") {
+ ntot <- rowSums(newresp)
+ ## FIXME: maybe unnecessary?
+ if (any(ntot==0)) stop("totals must be >0 ")
+ newresp <- newresp[,1]/ntot
+ rr$setWeights(ntot)
+ }
stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == length(rr$y))
rr$setResp(newresp)
pp <- object at pp$copy()
@@ -1173,7 +1181,8 @@
nAGQ <- dc$dims["nAGQ"]
nth <- dc$dims["nth"]
ff <- mkdevfun(list2env(list(pp=pp, resp=rr, u0=pp$u0, beta0=pp$beta0, verbose=0L,
- tolPwrss=dc$cmp["tolPwrss"], dpars=seq_len(nth))),
+ tolPwrss=dc$cmp["tolPwrss"], dpars=seq_len(nth),
+ nAGQ=unname(nAGQ))),
nAGQ=nAGQ)
xst <- rep.int(0.1, nth)
x0 <- pp$theta
@@ -1181,7 +1190,7 @@
if (!is.na(nAGQ) && nAGQ > 0L) {
xst <- c(xst, sqrt(diag(pp$unsc())))
x0 <- c(x0, pp$beta0)
- lower <- c(lower, rep(-Inf,length(pp$beta0)))
+ lower <- c(lower, rep(-Inf,length(x0)-length(lower)))
}
control <- list(...)$control
if (is.null(control)) control <- list()
@@ -1338,27 +1347,28 @@
val <- switch(family$family,
poisson=rpois(ntot,lambda=musim),
binomial={
- resp <- model.response(object at frame)
- bernoulli <- !is.matrix(resp)
- if (bernoulli) {
- rbinom(ntot,prob=musim,size=1)
- } else {
- nresp <- nrow(resp)
- ## FIXME: should "N-size" (column 2) be named?
- ## copying structures from stats/R/family.R
- sizes <- rowSums(resp)
- Y <- rbinom(ntot, size = sizes, prob = musim)
- YY <- cbind(Y, sizes - Y)
- yy <- lapply(split(YY,gl(nsim,nresp,2*nsim*nresp)),
- matrix, ncol=2,
- dimnames=list(NULL,colnames(resp)))
- ## colnames() <- colnames(resp)
- ## yy <- split(as.data.frame(YY),
- ## rep(1:nsim,each=length(sizes)))
- names(yy) <- paste("sim",seq_along(yy),sep="_")
- yy
- }
- },
+ w <- weights(object)
+ Y <- rbinom(ntot,prob=musim,size=w)
+ resp <- model.response(object at frame)
+ if (!is.matrix(resp)) { ## bernoulli, or weights specified
+ if (is.factor(resp)) {
+ if (any(weights(object)!=1)) stop("non-uniform weights with factor response??")
+ factor(levels(resp)[Y+1],levels=levels(resp))
+ } else {
+ Y/w
+ }
+ } else {
+ ## FIXME: should "N-size" (column 2) be named?
+ ## copying structures from stats/R/family.R
+ nresp <- nrow(resp)
+ YY <- cbind(Y, w - Y)
+ yy <- lapply(split(YY,gl(nsim,nresp,2*nsim*nresp)),
+ matrix, ncol=2,
+ dimnames=list(NULL,colnames(resp)))
+ names(yy) <- paste("sim",seq_along(yy),sep="_")
+ yy
+ }
+ },
stop("simulation not implemented for family",
family$family))
} else {
@@ -1557,7 +1567,7 @@
##' @keywords utilities
##' @examples
##'
-##' ## shows many methds you should consider *before* getME():
+##' ## shows many methods you should consider *before* getME():
##' methods(class = "merMod")
##'
##' (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
@@ -1966,3 +1976,8 @@
mc[[1]] <- as.name("dotplot.ranef.mer")
eval(mc)
}
+
+##' @S3method weights merMod
+weights.merMod <- function(object, ...) {
+ object at resp$weights
+}
Modified: pkg/lme4Eigen/inst/tests/test-lmer.R
===================================================================
--- pkg/lme4Eigen/inst/tests/test-lmer.R 2012-02-15 21:54:04 UTC (rev 1602)
+++ pkg/lme4Eigen/inst/tests/test-lmer.R 2012-02-21 02:16:21 UTC (rev 1603)
@@ -21,6 +21,7 @@
expect_that(extractAIC(fm1ML), equals(c(3, 333.327059881135)))
expect_that(vcov(fm1)[1,1], equals(375.720278729861))
expect_that(vcov(fm1ML)[1,1], equals(313.097224695739))
+ ## FIXME: recent version gets 313.09721874266512032 instead?
expect_that(fm2 <- refit(fm1, Dyestuff2$Yield), is_a("lmerMod"))
expect_that(fixef(fm2), is_equivalent_to(5.6656))
expect_that(VarCorr(fm2)[[1]][1,1], equals(0))
Modified: pkg/lme4Eigen/tests/simulate.R
===================================================================
--- pkg/lme4Eigen/tests/simulate.R 2012-02-15 21:54:04 UTC (rev 1602)
+++ pkg/lme4Eigen/tests/simulate.R 2012-02-21 02:16:21 UTC (rev 1603)
@@ -1,10 +1,46 @@
library(lme4Eigen)
-## debug(simulate.merMod)
+
+## binomial (2-column and prob/weights)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
-fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+gm2 <- glmer(incidence/size ~ period + (1 | herd), weights=size,
+ data = cbpp, family = binomial)
-simulate(gm1)
-simulate(fm1)
+s1 <- simulate(gm1,seed=101)[[1]]
+s2 <- simulate(gm2,seed=101)[[1]]
+stopifnot(all.equal(s1[,1]/rowSums(s1),s2))
-## FIXME: would like better tests ...
+## Bernoulli
+## works, but too slow
+if (FALSE) {
+ data(guImmun,package="mlmRev")
+ g1 <- glmer(immun~kid2p+mom25p+ord+ethn+momEd+husEd+momWork+rural+pcInd81+
+ (1|comm/mom),family="binomial",data=guImmun)
+ s2 <- simulate(g1)
+}
+
+d <- data.frame(f=rep(LETTERS[1:10],each=10))
+d$x <- runif(nrow(d))
+u <- rnorm(10)
+d$eta <- with(d,1+2*x+u[f])
+d$y <- rbinom(nrow(d),plogis(d$eta),size=1)
+
+g1 <- glmer(y~x+(1|f),data=d,family="binomial")
+s1 <- simulate(g1,seed=102)[[1]]
+
+d$y <- factor(c("N","Y")[d$y+1])
+g1B <- glmer(y~x+(1|f),data=d,family="binomial")
+s1B <- simulate(g1B,seed=102)[[1]]
+stopifnot(all.equal(s1,as.numeric(s1B)-1))
+
+## another Bernoulli
+data(Contraception,package="mlmRev")
+fm1 <- glmer(use ~ urban+age+livch+(1|district), Contraception, binomial)
+s3 <- simulate(fm1)
+
+d$y <- rpois(nrow(d),exp(d$eta))
+g2 <- glmer(y~x+(1|f),data=d,family="poisson")
+s4 <- simulate(g2)
+
+fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+s5 <- simulate(fm1)
More information about the Lme4-commits
mailing list