[R-sig-ME] binomial fixed-effect p-values by simulation
Daniel Ezra Johnson
danielezrajohnson at gmail.com
Mon Aug 25 00:30:38 CEST 2008
A quick update on my progress implementing the simulation significance
test. If anyone is interested I can try to provide reproducible
examples illustrating these points.
1) I found that occasionally, refit(model.mer,response.vector)
produced a different result than glmer(response.vector ~ ...). My
solution is to use glmer() rather than refit() for each simulation,
and it doesn't appear to make it that much slower.
2) I found that the LRT approach does not seem to work as well as
simulating the z-statistic. In the LRT approach, occasionally either
the null model or the alternative model would come out with a "crazy"
log-likelihood -- e.g. one LL would be -400, the other -21000 -- and
this happened often enough in the right direction to throw off the
empirical p-values. So I'm going ahead with the z-statistic approach
(the same datasets don't give a "crazy" z-statistic).
3) There's still the problem of what to do when the simulated model
does not converge or has a false convergence - whether or not it's
biasing the empirical p-value to throw away those runs. It probably
is, in that whether or not the false convergence occurs is probably
not unrelated to the size of the fixed effect. But it doesn't appear
to be a direct relationship, so I'm comfortable tossing those runs for
Probably all of these things are happening more often than they
usually would because I'm testing this on a rather extreme type of
data: binomial, 20 subjects with a large subject standard deviation
(5) and only 100 observations per subject, so there end up being a lot
of subjects with invariant responses, 0% or 100%.
And since my "observed" data has a Wald fixed effect p-value of 0.0005
(the p-value from the z-score), it's not that easy to judge whether
the simulation is slightly more conservative in a reasonable amount of
Well, 1000 runs just finished and I got 2 where abs(z) >= z.obs, so
empirical p is 0.002.
I can provide more details if anyone is interested. Thanks again for
your help with this.
More information about the R-sig-mixed-models