[R-sig-ME] Fwd: Re: Power calculations for lmer mixed designs

[Phillip Alday]

Hi Diane,

On 06/19/2017 03:06 PM, Diana Michl wrote:
> Dear List,
> 
> I want to do a power calculation on a mixed model as described here by 
> Ben Bolker: https://rpubs.com/bbolker/11703
> 
> However, there are three things I can't figure out:
> 
> - The qlogis function, I don't get what it does exactly and what to 
> enter instead of 0.7 as in the example
> 
> |beta <- c(qlogis(0.7), -0.2)|

See the paragraph above that code: "baseline range (fraction of grid
cells sampled) is 70% (logit(0.7)=0.8473)". qlogis() is the quantile
function for the logistic distribution, so qlogis returns the point in
that distribution where the density is equal to the value mentioned.
This is used here to compute logit(0.7):

> qlogis(0.7)
[1] 0.8472979

Which is the value in the paragraph given for logit(0.7). This is just
converting effect size to the link (logit) scale (as mentioned in the
introduction: ' "beta" is the fixed-effects parameters, in this case
(intercept,treat) – also all on the logit scale.').

Note that this is the intercept term ("baseline range").

> - What exactly the resulting graphs mean, both here and in my own 
> example. What does the read bar at -0.2 tell me? 

Skimming the text, -0.2 is the coefficient/coefficient for the treatment
(already on the link scale -- note that it is not transformed, so it
represents log odds). In the explanatory text, "decrease log-odds of
range by 0.4", i.e. -0.4 is used -- maybe Ben Bolker can explain why the
number changed / what obvious transformation step I'm missing.

So the redline represents the "ground truth", i.e. the real parameter
value (because you chose it) and the other plots tell your simulated
estimates -- how close they were (see esp. the CI plot).

> Is there a way to give power as a number between 0 and 1?
> - How to do this for more than 1 treatment (my attempts to do it for 4 
> all failed, I seem to miss places in the exampe I need to change. For 
> example, I tried adding values to the beta-qlogis line above and others, 
> but it gave errors)

You calculate your power based on the ability of models you fit to the
simulated datasets to detect the effect. The general work flow is:

1. define "ground truth" model based on some theoretical assumptions
2. simulate lots of datasets of a certain size
3. fit models to those datasets.
4. number of models able to detect effect / number of models = power
5. if power to low, repeat with bigger simulated datasets

You can do this all by hand, but there are several tools for automating
this. I've tinkered with the simr package in the past and found it nice
(even nice enough to contribute a few things because I plan on using it
again).

There are also "analytic" tools that try to compute power without
simulation, see for exmaple Jake Westfall's Shiny apps (and read his
papers on power issues in mixed effects models!):

https://jakewestfall.shinyapps.io/pangea/
https://jakewestfall.shinyapps.io/crossedpower/

Best,
Phillip

> Thanks very much for any input!
> 
> 
> Diana
>