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

[Phillip Alday]

On 06/22/2017 01:43 PM, Diana Michl wrote:
::snip::

> My intercept and effect size is already on the logit scale - the
> intercept is 7.2. But whenever I enter it directly into the qlogis
> function, I'm warned that NaNs are produced. This happens whenever I
> enter anything over 1.0. Does this mean I actually need to enter a
> percentage rather than the absolute value...? Or is the problem an
> entirely different one?

I'm not sure this makes sense. Unfortunately, the terminology in
statistics for things like "link" and "response" scale is rather
confusing and unintuitive at times and it's easy to get mixed up which
scale we're talking about at the moment (0s and 1s in the dependent
variable? the scale of the coefficients, whatever that is called?). Any
of the q*() functions will return NAs if you give it something bigger
than 1.0 because 1.0 is per definition the maximum probability. So
you're asking for the value where the cumulative probability exceeds the
maximum probability, which is of course "huh, what?" or in R terms, NA.

If you know your coefficient is supposed to be 7.2 based on a previous
model, then it's already on the correct scale and you don't need to
convert it to logit scale. If you're starting off from percentages, then
you need to convert it.

(I had a mistake in my previous message -- the q*() functions don't give
the density, it gives the place where the cumulative probability. So
qlogis(0.5) give 0 because the probability of having a value between
-Inf and 0 on a logistic distribution is 0.5)

>
::snip::
>> 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).
> So if I'm seeing this correctly, this way of calculating power is not a
> direct one, but you're rather meant to draw your own conclusions from
> looking at the graphs of the simulations and get a general feel of how
> far off or on your effect size is?

More or less. You can also use the simulation results to directly give
you a power estimate without a graphical presentation. But yes, it's a
simulation based approach and not an analytic (formula-based) one.


>>> 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
> So you can't properly attempt this very way with 4 effects? The author
> implies that you can, but I'm getting the impression it'd be quite
> complicated... I'm not sure running each effect separately would be
> correct.

The interaction of main effects is just another effect from the model's
perspective, so you can model these separately. Each simulation includes
random draws from all of the model terms, so in that sense, you're
modelling everything at once. The CRAN version of the simr package
(which makes simulation-based power analyses a lot easier) unfortunately
requires a separate set of simulation runs for each model term. I had
been working on simultaneous testing, i.e. re-using simulation runs for
all the terms, but for various reasons never got around to making it
stable enough for the CRAN version. You can install it with devtools:

> devtools::install_github("palday/simr",ref="simultaneous")
> library(simr)
> ?powerSimSimultaneous
> ?powerSimMultiple


Best,
Phillip

>> 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/
> I'll try those, thanks very much!
> 
> Best
> Diana