# umx!

## Power

Power is the probability of detecting an effect when it is present. Knowing the power of a study is important: Weak studies spew out 5% false positives (i.e., at the rate alpha), but without power, they fail to detect true results at rate beta. Power is 1-beta.

This post covers:

1. Computing power to detect a difference between 2 models.
2. Computing power given a certain `n`.

`umx` supports power calculations for the `umxACE` model (via `power.ACE.test`) and a generic `umxPower` function, which can handle a range of power-related questions given a trueModel, and a parameter which you wish to drop.

## Comparing two models with umxPower

`R` has great support for power. For instance, testing power to detect a correlation of .3 in a sample of 36 subjects, R quickly tells us we’re going to need some more subjects…

``````pwr::pwr.r.test(n = 36, r = .3)
n = 36
r = 0.3
sig.level = 0.05
power = 0.4361314
alternative = two.sided
``````

Computing power for SEM is rather more complex, as the models are arbitrary. Many more questions can also be asked than, say, `pwr.t.test` could answer.

`umxPower` tackles this by using an interface very similar to `umxModify`. You offer up a true model, and paths that you want test.

``````umxPower(trueModel = , update = "path_to_test", n= , power = , sig.level = .05, method = c("ncp", "empirical"))
``````

To emulate the `pwr.r.test` example, you would build a `trueModel` of two manifests (“X”, and “Y”), with a covariance between them equal to r= .3. You can then input this to `umxPower` test, setting update to “X_with_Y” (value= 0 is the default for updated parameters).

### Demonstration: Power to detect a correlation of .3 between two variables: X and Y

First: Make and run a model of the true correlation (.3)

Make data for the model (this can be tricky for more complex models. the easiest solution is often to use the model with mxGenerateData to generate data corresponding to the model).

``````tmp = umx_make_raw_from_cov(qm(1, .3| .3, 1), n=200, varNames= c("X", "Y"), empirical= TRUE)
``````

Build the model: in this case a correlation

``````m1 = umxRAM("corXY", data = tmp,
umxPath("X", with = "Y"),
umxPath(var = c("X", "Y"))
)

``````

Next: Using `umxPower`, test the path you want to drop,

``````umxPower(m1, "X_with_Y", n= 50, method="empirical")
``````

After a few moments, this yields:

``````method = empirical
n = 50  sig.level = 0.05
power = 0.6266667
probes = 300  statistic = LRT
``````

We can also take advantage of the noncentrality parameter, to run power. This is near-instant and accurate if the model is valid.

``````umxPower(m1, "X_with_Y", n= 50, method="ncp")
method = ncp
n = 50
sig.level = 0.05
power = 0.5837944
statistic = LRT
``````

Now, let’s compare the results using a `cor.test` doing the same thing?

``````pwr::pwr.r.test(n = 50, r = .3)
n = 50
r = 0.3
sig.level = 0.05
power = 0.5715558
alternative = two.sided
``````

### Another example

If someone publishes this simple model claiming that larger engines lower miles/gallon in vehicles, what power do we have to detect `engine_litres_to_mpg` is > 0 if we replicate with the same `n`, and assume that the obtained effect is the true one?

First we build the claimed model:

``````df = mtcars
df\$engine_litres = .016 * df\$disp

m1 = umxRAM(data = df, "#mpg_model
mpg ~ engine_litres + wt
wt ~~ engine_litres"
)
``````

And test power:

``````umxPower(m1, "engine_litres_to_mpg", n= 30, method = "empirical")

method = empirical
n = 30
power = 0.46
probes = 300
sig.level = 0.05
statistic = LRT
``````

Now let’s solve for N to give us 80% power.

``````# solve for N
umxPower(m1, update = "engine_litres_to_mpg", power = .8, method = "ncp")

####################
# Estimating n #
####################

method = ncp
n = 68
power = 0.8
sig.level = 0.05
statistic = LRT

``````

Finally, we can make a table showing how power changes across n:

``````
tmp = umx_make_raw_from_cov(qm(1, .3| .3, 1), n= 200, varNames= c("X", "Y"), empirical= TRUE)

m1 = umxRAM("corXY", data = tmp,
umxPath("X", with = "Y"),
umxPath(var = c("X", "Y"))
)

umxPower(m1, update = "X_with_Y", explore = TRUE)

N power
1   17.68  0.25
2   25.36  0.34
3   33.03  0.42
4   40.70  0.50
5   48.38  0.57
6   56.05  0.63
7   63.72  0.69
8   71.40  0.74
9   79.07  0.78
10  86.74  0.82
11  94.42  0.85
12 102.09  0.87
13 109.76  0.90
14 117.44  0.91
15 125.11  0.93
16 132.78  0.94
17 140.46  0.95
18 148.13  0.96
19 155.80  0.97
20 163.48  0.98

``````