
Thanks very much explaining so clearly why this is complicated! Your help is much appreciated. Given that parametric bootstrapping a binomial glmmPQL model is so tricky, I think I might back out and go down the GEE route instead. GEEs give similar parameter estimates to glmmPQL for my problem, and in light of the comments I've received would appear to be a more suitable choice.

Thanks again.

> 
> From: Ben Bolker <bbolker@gmail.com>
> Subject: Re: [R-sig-ME] Simulate samples from glmmPQL
> Date: 24 October 2012 10:35:52 NZDT
> To: r-sig-mixed-models@r-project.org
> 
> 
> Andrew Digby <andrewdigby@...> writes:
> 
>> I want to use a parametric bootstrap on a GLMM to to determine
>> parameter bias
>> (http://www.bris.ac.uk/cmm/team/hg/bias-correction.pdf) and
>> hopefully to estimate effect confidence intervals
>> (http://glmm.wikidot.com/faq).
> 
>> I'm using glmmPQL{MASS} since I need to include an AR1 correlation
>> structure, but there are no existing parametric bootstrap routines
>> for this. I'm writing my own, but am coming unstuck when trying to
>> simulate a response from the model:
> 
>> m1<-glmmPQL(y ~ trt , random = ~ 1 | ID, correlation=corAR1(form=~week), 
>> family = binomial, data = bacteria)
>> simulate(m1)
>> 
>> 	Error in simulate.lme(m1) : 
>>  		Models with corStruct and/or varFunc objects not allowed.
> 
>> I'm trying to write my own simulation code by modifying that from
>> lme4 [getMethod('simulate',signature='mer')]. I'm not sure if or how
>> I should be including the correlation structure though, and guess
>> there's a very good reason why simulate.lme doesn't work with a
>> corStruct.
> 
>> If anyone can give me an idea of how to simulate samples from a
>> glmmPQL object with a corStruct, or if it's possible, then I'd be
>> very grateful. Or just a way of obtaining a rough estimate of the
>> bias would be great.  Would a non-parametric bootstrap be completely
>> out of the question here
>> (e.g.
>> http://cran.r-project.org/doc/contrib/Fox-Companion/
>>  appendix-bootstrapping.pdf)?
> 
>  This seems tricky to me.
> 
> 1. generating correlated _Gaussian_ variables is no big deal (you can use
> the RandomFields package, which is a bit of a sledgehammer/nutcracker scenario
> for an AR(1) variable, or more simply construct your variance-covariance matrix
> and use MASS::mvrnorm), but ...
> 
> 2. Integrating this in order to write your own version of simulate.lme() that
> did corStruct models might be hairy -- the internals of simulate.lme() are
> a little complicated (but it still might be the best solution).
> 
> 3. As Doug Bates has said here, it's not entirely clear whether the
> model glmmPQL is fitting (essentially, a marginal model of the
> variance-covariance structure) actually corresponds to anything you
> could actually simulate.  A well-formed model of this sort would look
> like:
> 
>   y ~ Binomial(plogis(X %*% beta + Z %*% b + eps))
>   b ~ Sigma_b (from G-side/block structure)
>   eps ~ MVN(AR(1))
> 
> this is a little crude, but the idea is that you actually break the model into 
> well-defined hierarchical bits.  It's not clear (to me at least) whether the
> model glmmPQL fits in this case actually corresponds to this statistical model
> or not.  (INLA appears like it might be the best way forward on this front ...)
> (The advantage of hacking simulate.lme is that you'd be taking advantage
> of whatever simulate.lme is doing ...)
> 
> 4. It turns out that corAR1 is not the only problem: even a 'plain' glmmPQL
> call (without any correlations) uses a varStruct internally, which breaks
> simulate.lme (there is no simulate.glmmPQL, so simulate falls back on
> simulate.lme).  It seems to me that if you could hack simulate.lme you
> could get a result with the right _variance-covariance structure_, but it
> wouldn't be Bernoulli -- it would be normally distributed.  And, just
> to add insult to injury, simulate.lme doesn't work the way simulate
> methods work in the rest of R -- simulate.lme was written earlier,
> and does a bunch of refitting (essentially, parametric bootstrapping!)
> stuff rather than returning the raw simulation values.
> 
> 
> 
> 
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]

