[Lme4-commits] r1725 - in pkg/lme4.0: . src tests

[noreply at r-forge.r-project.org]

Author: bbolker
Date: 2012-05-09 17:23:45 +0200 (Wed, 09 May 2012)
New Revision: 1725

Added:
   pkg/lme4.0/tests/agq.R
Modified:
   pkg/lme4.0/DESCRIPTION
   pkg/lme4.0/src/lmer.c
   pkg/lme4.0/tests/getME.R
Log:

  Wayne Zhang's fixes for AGQ ("backport")
  note that there are still potential outstanding issues with the multiple-RE case (!)



Modified: pkg/lme4.0/DESCRIPTION
===================================================================
--- pkg/lme4.0/DESCRIPTION	2012-05-08 16:49:10 UTC (rev 1724)
+++ pkg/lme4.0/DESCRIPTION	2012-05-09 15:23:45 UTC (rev 1725)
@@ -1,6 +1,6 @@
 Package: lme4.0
-Version: 0.9999-1
-Date: 2012-03-22
+Version: 0.9999-2
+Date: 2012-05-08
 Title: Linear mixed-effects models using S4 classes
 Description: Fit linear and generalized linear mixed-effects models.
   This is the implementation of lme4 available on CRAN and developed up to 2011.

Modified: pkg/lme4.0/src/lmer.c
===================================================================
--- pkg/lme4.0/src/lmer.c	2012-05-08 16:49:10 UTC (rev 1724)
+++ pkg/lme4.0/src/lmer.c	2012-05-09 15:23:45 UTC (rev 1725)
@@ -1376,13 +1376,13 @@
 	    double *ans = Calloc(nl, double);    /* current penalized residuals in different levels */
 
 				/* update abscissas and weights */
+	    /* fixes from Wayne Zhang */
 	    for(int i = 0; i < nre; ++i){
 		for(int j = 0; j < nl; ++j){
 		    z[i + j * nre] = ghx[pointer[i]];
 		}
 		w_pro *= ghw[pointer[i]];
-		if(!MUETA_SLOT(x))
-		  z_sum += z[pointer[i]] * z[pointer[i]];
+		z_sum += ghx[pointer[i]] * ghx[pointer[i]];
 	    }
 
 	    CHM_DN cz = N_AS_CHM_DN(z, q, 1), sol;
@@ -1391,7 +1391,10 @@
 	    Memcpy(z, (double *)sol->x, q);
 	    M_cholmod_free_dense(&sol, &c);
 
-	    for(int i = 0; i < q; ++i) u[i] = uold[i] + sigma * z[i];
+	    /* fix from Wayne Zhang */
+	    for(int i = 0; i < q; ++i) 
+		u[i] = uold[i] + sqrt(2) * sigma * z[i];
+
 	    update_mu(x);
 
 	    AZERO(ans, nl);
@@ -1404,7 +1407,7 @@
 
 	    for(int i = 0; i < nl; ++i)
 		tmp[i] += exp( factor * ans[i] + z_sum) * w_pro / sqrt(PI);
-				/* move pointer to next combination of weights and abbsicas */
+	    /* move pointer to next combination of weights and abscissas */
 	    int count = 0;
 	    pointer[count]++;
 	    while(pointer[count] == nAGQ && count < nre - 1){

Added: pkg/lme4.0/tests/agq.R
===================================================================
--- pkg/lme4.0/tests/agq.R	                        (rev 0)
+++ pkg/lme4.0/tests/agq.R	2012-05-09 15:23:45 UTC (rev 1725)
@@ -0,0 +1,27 @@
+library(lme4.0)
+
+## test AGQ=10 against glmmML results
+g1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+                         family = binomial, data = cbpp, nAGQ=10)
+getinfo.merMod <- function(m) c(fixef(m),deviance=deviance(m))
+
+i1 <- getinfo.merMod(g1)
+
+if (FALSE) {
+    ## comment to avoid R CMD check complaints
+    ## library(glmmML)
+    getinfo.glmmML <- function(m) c(coef(m),deviance=deviance(m))
+    i2 <- getinfo.glmmML(glmmML(cbind(incidence, size - incidence) ~ period,
+                                family = binomial,
+                                cluster=herd,
+                                method="ghq",
+                                n.points=10,
+                                data = cbpp))
+}
+i2 <- structure(c(-1.39924138019006, -0.991381753243723, -1.12782730029348, 
+                  -1.57948092000465, 100.010029977086), .Names = c("(Intercept)", 
+                                                        "period2", "period3", "period4", "deviance"))
+
+stopifnot(all.equal(i1,i2,tol=1e-6))
+
+

Modified: pkg/lme4.0/tests/getME.R
===================================================================
--- pkg/lme4.0/tests/getME.R	2012-05-08 16:49:10 UTC (rev 1724)
+++ pkg/lme4.0/tests/getME.R	2012-05-09 15:23:45 UTC (rev 1725)
@@ -46,3 +46,10 @@
 chkMEs(fm4, nmME)
 
 isREML(fm1)
+
+L <- as(fm1 at L,"Matrix")
+Z <- getME(fm1,"Z")
+A <- fm1 at A
+dim(L)
+dim(Z)
+dim(A)