Hello,
I'm fitting group-level data for success, failure of a given exercise using
a binomial-normal model. My data look like this:
> Exer
   Region Group Size Success
1       1     1   51       4
2       1     2  102      14
3       1     3  104      13
4       1     4  103      11
5       2     1  124      25
6       2     2  125      15
7       3     1   45      11
8       3     2   44       8
9       3     3   44      19
10      4     1   89      13
...

The model I choose is:
> m1<-glmer(cbind(Success, Size-Success)~-1+Group+(1|Region),
family="binomial", Exer)

The model converges and gives reasonable parameter estimates.
Now I want to simulate the Success rate for the same regions and groups
than the ones I have in my Exer dataset, while taking into account the
between-Region variability (i.e. not equivalent to fitted(m1)):

require(arm)
predict.lmerBin<-function(object, X, Zt, z){
# @X: model matrix for the fixed-effects
# @Zt: transpose of model matrix for the random effects
# z: random effect (etas)
b<-fixef(object)
#      plogis(as.vector(X%*%b+t(Zt)%*%z)) # As proposed by Thierry Okelinx
invlogit(as.vector(X%*%b+t(Zt)%*%z))         #invlogit() from arm package
}

Question 1: In a previous post (
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q4/004858.html),
Thierry Onkelinx used plogis() instead of invlogit(). Can anyone shade some
light on why plogis()?

omega<-sigma.hat(m1)$sigma$Region        # SD (region) from model fit
X<- m1@X         # Keep the same design matrix as original
Zt<- m1@Zt
netas<-length(unique(Exer$Region))
nSim<-1000

my.predict<-function(indx){
newZ<-rnorm(netas, 0, omega)         # Btw region random effect
return(list(predict.lmerBin(finalmod, X, Zt, newZ)))
}

predblock<-as.vector(t(do.call("rbind", sapply(1:nSim, function(i)
my.predict(i))))) # Run the prediction nSim times
desgblock<-Exer[rep(row.names(Exer), nSim),)] # Design nSim times
preddf<-data.frame(desgblock, predblock)
# Extract the mean, median predicted success rate per Group
require(Hmisc)
require(plyr)
fmpi<-ddply(preddf, .(Group), summarise,
avg=wtd.mean(predblock, weights=Size, normwt=TRUE),
p500=wtd.quantile(predblock, weights=Size, .500, normwt=TRUE))

I'm NOT 100% sure that this code is correct. Any sanity check or insight on
how to make this code more robust is welcome.
Thanks a lot for your help,

Regards

	[[alternative HTML version deleted]]


