Thanks again for spotting that problem, Vito. Because of the need for
generality, the C code in the mer_postVar function ends up being
fairly involved. Once I traced through the example, however, I could
readily spot that I was running a loop over the levels of each of the
grouping factors but I wasn't using that loop index to change the
input values. Therefore, the result was exactly what you observed,
those values all ended up corresponding to the results for the first
level of the grouping factor.
Do the results in the enclosed file look like what you expect? I
rewrote the code that you sent because I found some aspects of that
code to be rather frightening. (Is that code directly from the book?)
I will commit the changes to the SVN archive on R-forge by this
evening (North American time). I have been tearing apart my code
again and right now I have lmer and glmer running but nlmer is having
problems. I'll close out the bug report after I commit the changes.
On Jan 17, 2008 7:42 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
Thanks for the thorough bug report, Vito. I have added it to the bug
report list at http://r-forge.r-project.org/projects/lme4 (under the
Tracker tab). Those with R-forge logins can add bug reports there
directly if they wish. Those who do not have an R-forge login should
consider getting one.
I appreciate the example. It will make it a lot easier to isolate the problem.
On Jan 17, 2008 5:39 AM, vito muggeo <vmuggeo at dssm.unipa.it> wrote:
Dear all,
I'm running data from "Gelman&Hill Data Analysis Using regression..
pag259-261" (you can get data by running the code below). In particular
I'm interested in getting variances of the predictions of the random
intercepts. I get "right" (as compared with those from the book) results
using lme 0.99875-9 (downloaded from CRAN).
However using the latest lme4_0.999375-0 obtained via
install.packages("lme4", repos = "http://r-forge.r-project.org")
the variances of the predictions are wrong..(they are constant
regardless of different cluster sizes. You can check this running the
code below..
Many thanks,
vito
srrs2 <-
read.table("http://www.stat.columbia.edu/~gelman/arm/examples/radon/srrs2.dat",
header=T, sep=",")
mn <- srrs2$state=="MN"
radon <- srrs2$activity[mn]
log.radon <- log (ifelse (radon==0, .1, radon))
floor <- srrs2$floor[mn] # 0 for basement, 1 for first floor
n <- length(radon)
y <- log.radon
x <- floor
# get county index variable
county.name <- as.vector(srrs2$county[mn])
uniq <- unique(county.name)
J <- length(uniq)
county <- rep (NA, J)
for (i in 1:J){
county[county.name==uniq[i]] <- i
}
library(arm)
M1<-lmer(y~1+x+(1|county))
ranef(M1)[[1]][1:10,] #OK
se.ranef(M1)[[1]][1:10,] #wrong results: they are constant!!!??!!!
R version 2.6.1 (2007-11-26)
i386-pc-mingw32
locale:
LC_COLLATE=Italian_Italy.1252;LC_CTYPE=Italian_Italy.1252;LC_MONETARY=Italian_Italy.1252;LC_NUMERIC=C;LC_TIME=Italian_Italy.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] foreign_0.8-23 car_1.2-7 arm_1.1-1
R2WinBUGS_2.1-6 coda_0.13-1 lme4_0.999375-0 Matrix_0.999375-4
lattice_0.17-2
[9] MASS_7.2-38
loaded via a namespace (and not attached):
[1] grid_2.6.1
--
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Universit? di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612