-
Notifications
You must be signed in to change notification settings - Fork 80
/
Copy path11-power_tests.Rmd
217 lines (148 loc) · 21.9 KB
/
11-power_tests.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# Power Analysis to Estimate Required Sample Sizes for Modeling {#power-tests}
`r if (knitr::is_latex_output()) '\\index{power analysis|(}'`
In the vast majority of situations in people analytics, researchers and analysts have limited control over the size of their samples. The most common situation is, of course, that analyses are run with whatever data can be gleaned and cleaned in the time available. At the same time, as we have seen in all of our previous work, even if a certain difference might exist in real life in the populations being studied, it is by no means certain that a specific analysis on samples from these populations will elucidate that difference. Whether that difference is visible depends on the statistical properties of the samples used. Therefore, researchers and analysts are living in the reality that when they conduct inferential analysis, the usefulness of their work depends to a very large degree on the samples they have available.
This suggests that a conscientious analyst would be well advised to do some up-front work to determine if their samples have a chance of yielding results that are of some inferential value. In a practical context, however, this is only partly true (and that is an important reason why this chapter has been left towards the end of this book). Estimating required sample sizes is an imprecise science. Although the mathematics suggest that in theory it should be precise, in reality we are guessing most of the inputs to the mathematics. In many cases we are so clueless about those inputs that we move into the realms of pure speculation and produce ranges of required sample sizes that are so wide as to be fairly meaningless in practice.
That said, there are situations where conducting *power analysis*---that is, analysis of the required statistical properties of samples in order to have a certain minimum probability of observing a true difference---makes sense. Power analysis is an important element of experimental design. Experiments in people analytics usually take one of two forms:
`r if (knitr::is_latex_output()) '\\index{experiments}'`
1. *Prospective experiments* involve running some sort of test or pilot on populations to determine if a certain measure has a hypothesized effect. For example, introducing a certain new employee benefit for a specific subset of the company for a limited period of time, and determining if there was a difference in the impact on employee satisfaction compared to those who did not receive the benefit.
2. *Retrospective experiments* involve the use of historical data to test if a certain measure has a hypothesized effect. This usually occurs opportunistically when it is apparent that a certain measure has occurred in the past and for a limited time, and data can be drawn to test whether or not that measure resulted in the hypothesized effect.
Both prospective and retrospective experiments can involve a lot of work---either in setting up experiments or in extracting data from history. There is a natural question as to whether the chances of success justify the required resources and effort. Before proceeding in these cases, it is sensible to get a point of view on the likely power of the experiment and what level of sample size might be needed in order to establish a meaningful inference. For this reason, power analysis is a common component of research proposals in the medical or social sciences.
Power analysis is a relatively blunt instrument whose primary value is to make sure that substantial effort is not being wasted on foolhardy research. If the analyst already has reasonably available data and wants to test for the effect of a certain phenomenon, the most direct approach is to just go and run the appropriate model assuming that it is relatively straightforward to do so. Power analysis should only be considered if there is clearly some substantial labor involved in the proposed modeling work.
## Errors, effect sizes and statistical power
`r if (knitr::is_latex_output()) '\\index{power analysis!logic and principles}'`Before looking at practical ways to conduct power tests on proposed experiments, let's review an example of the logical and mathematical principles behind power testing, so that we understand what the results of power tests mean. Recall from Section \@ref(hyp-tests) the logical mechanisms for testing hypotheses of statistical difference. Given data on samples of two groups in a population, the *null hypothesis* $H_0$ is the hypothesis that a difference does not exist between the groups in the overall population. If the null hypothesis is rejected, we accept the alternative hypothesis $H_1$ that a difference does exist between the groups in the population.
Recall also that we use the statistical properties of the samples to make inferences about the null and alternative hypotheses based on statistical likelihood. This means that four possible situations can occur when we run hypothesis tests:
1. We fail to reject $H_0$, and in fact $H_1$ is false. This is a good outcome.
2. We reject $H_0$, but in fact $H_1$ is false. This is known as a *Type I error*.
3. We fail to reject $H_0$, but in fact $H_1$ is true. This is known as a *Type II error*.
4. We reject $H_0$, and in fact $H_1$ is true. This is a good outcome and one which is most often the motivation for the hypothesis test in the first place.
`r if (knitr::is_latex_output()) '\\index{statistical power}'`*Statistical power* refers to the fourth situation and is the *probability that $H_0$ is rejected and $H_1$ is true*. Statistical power depends *at a minimum* on three criteria:
* The significance level $\alpha$ at which the analysis wishes to reject $H_0$ (see Section \@ref(hyp-tests)). Usually $\alpha = 0.05$.
* The size $n$ of the sample being used.
* The size of the difference observed in the sample, known as the *effect size*. There are numerous definitions of the effect size that depend on the specific type of power test being conducted.
As an example to illustrate the mathematical relationship between these criteria, let's assume that we run an experiment on a group of employees of size $n$ where we introduce a new benefit and then test their satisfaction levels before and after its introduction. As a statistic of a random variable, we can expect the mean difference in satisfaction to have a normal distribution. Let $\mu_0$ be the mean of the population under the null hypothesis and let $\mu_1$ be the mean of the population under the alternative hypothesis. Now let's assume that in our sample we observe a mean satisfaction of $\mu^*$ after the experiment. Recall from Chapter \@ref(found-stats) that to meet a statistical significance standard of $\alpha$, we will need $\mu^*$ to be greater than a certain multiple of the standard error $\frac{\sigma}{\sqrt{n}}$ above $\mu_0$ based on the normal distribution. Let's call that multiple $z_{\alpha}$. Therefore, we can say that the statistical power of our hypothesis test is:
$$
\begin{aligned}
\mathrm{Power} &= P(\mu^* > \mu_0 + z_{\alpha}\frac{\sigma}{\sqrt{n}}\vert{\mu = \mu_1}) \\
&= P(\frac{\mu^* - \mu_1}{\frac{\sigma}{\sqrt{n}}} > -\frac{\mu_1 - \mu_0}{\frac{\sigma}{\sqrt{n}}} + z_{\alpha}\vert{\mu = \mu_1}) \\
&= 1 - \Phi(-\frac{\mu_1 - \mu_0}{\frac{\sigma}{\sqrt{n}}} + z_{\alpha}) \\
&= 1 - \Phi(-\frac{\mu_1 - \mu_0}{\sigma}\sqrt{n} + z_{\alpha}) \\
&= 1 - \Phi(-d\sqrt{n} + z_{\alpha})
\end{aligned}
$$
where $\Phi$ is the cumulative normal probability distribution function, `r if (knitr::is_latex_output()) '\\index{Cohen\'s effect size}'`and $d = \frac{\mu_1 - \mu_0}{\sigma}$ is known as *Cohen's effect size*. Therefore, we can see that power depends on a measure of the observed effect size between our two samples (defined as Cohen's $d$) the significance level $\alpha$ and the sample size $n$^[We will also need to know the expected distribution of the statistics that we are analyzing in order to determine the power probability.].
The reader may immediately observe that many of these measures are not known at the typical point at which we would wish to do a power analysis. We can assert a minimum level of statistical power that we would wish for---usually this is somewhere between 0.8 and 0.9. We can also assert our $\alpha$. But at a point of experimental design, we usually do not know the sample size and we do not know what difference would be observed in that sample (the effect size). This implies that we are dealing with a single equation with more than one unknown, and this means that there is no unique solution^[In reality there are more unknowns that this math would imply, due to the imperfection of what we are trying to measure. For example measurement error and reliability will often be an unmeasurable unknown. For this reason you will often need a larger sample size than that indicated by power tests.]. Practically speaking, looking at ranges of values will be common in power analysis.
## Power analysis for simple hypothesis tests {#simple-stats}
Usually we will run power analyses to get a sense of required sample sizes. Given the observations on unknowns in the previous section, we will have to assert certain possible statistical results in order to estimate required sample sizes. Most often, we will need to suggest the observed effect size in order to obtain the minimum sample size for that effect size to return a statistically significant result at a desired level of statistical power.
`r if (knitr::is_latex_output()) '\\index{R!packages!WebPower@\\texttt{WebPower}}'`Using our example from the previous section, let's assume that we would see a 'medium' effect size on our samples. *Cohen's Rule of Thumb* for $d$ states that $d = 0.2$ is a small effect size, $d = 0.5$ a medium effect size and $d = 0.8$ a large effect size. We can use the `wp.t()` function from the `WebPower` package in R to do a power analysis on a paired two-sample $t$-test and return a minimum required sample size. We can assume $d = 0.5$ and that we require a power of 0.8---that is, we want an 80% probability that the test will return an accurate rejection of the null hypothesis.
```{r}
library(WebPower)
# get minimum n for power of 0.8
(n_test <- WebPower::wp.t(d = 0.5, p = 0.8, type = "paired"))
```
This tells us that we need an absolute minimum of `r ceiling(n_test$n)` individuals in our sample for an effect size of 0.5 to return a significant difference at an alpha of 0.05 with 80% probability. Alternatively we can test the power of a specific proposed sample size.
```{r}
# get power for n of 40
(p_test <- WebPower::wp.t(n1 = 40, d = 0.5, type = "paired"))
```
This tells us that a minimum sample size of 40 would result in a power of `r round(p_test$p, 2)`. A similar process can be used to plot the dependence between power and sample size under various conditions as in Figure \@ref(fig:powersampleplot). This is known as a *power curve*.
```{r powersampleplot, fig.cap="Plot of power against sample size for a paired t-test", fig.align="center", out.width = if(knitr::is_latex_output()) {"90%"}}
# test a range of sample sizes
sample_sizes <- 20:100
power <- WebPower::wp.t(n1 = sample_sizes, d = 0.5, type = "paired")
plot(power)
```
We can see a 'sweet spot' of approximately 40--60 minimum required participants, and a diminishing return on statistical power over and above this. Similarly we can plot a proposed minimum sample size against a range of effect sizes as in Figure \@ref(fig:powerdplot).
```{r powerdplot, fig.pos = "!h", fig.cap="Plot of power against effect size for a paired t-test", fig.align="center", out.width = if(knitr::is_latex_output()) {"90%"}}
# test a range of effect sizes
effect_sizes <- 2:8/10
samples <- WebPower::wp.t(n1 = rep(40, 7),
d = effect_sizes,
type = "paired")
plot(samples$d, samples$power, type = "b",
xlab = "Effect size", ylab = "Power")
```
Similar power test variants exist for other common simple hypothesis tests. Let's assume that we want to institute a screening test in a recruiting process, and we want to validate this test by running it on a random set of employees with the aim of proving that the test score has a significant non-zero correlation with job performance. If we assume that we will see a moderate correlation of $r = 0.3$ in our sample^[Cohen's rule of thumb for correlation coefficients is Weak: 0.1, Moderate: 0.3 and Strong: 0.5.], we can use the `wp.correlation()` function in `WebPower` to do a power analysis, resulting in Figure \@ref(fig:correlpower).
```{r correlpower, fig.cap="Plot of power against sample size for a correlation test", fig.align="center", out.width = if(knitr::is_latex_output()) {"90%"}}
sample_sizes <- 50:150
correl_powers <- WebPower::wp.correlation(n = sample_sizes, r = 0.3)
plot(correl_powers)
```
Figure \@ref(fig:correlpower) informs us that we will likely want to be hitting at least 100 employees in our study to have any reasonable chance of establishing possible validity for our screening test.
## Power analysis for linear regression models
In power tests of linear regression models, the effect size is a statistic of the difference in model fit between the two models being compared. Most commonly this will be a comparison of a 'full' fitted model involving specific input variables compared to a 'reduced' model with fewer input variables (often a random variance model with no input variables).
`r if (knitr::is_latex_output()) '\\index{power analysis!f-squared}'`The $f^2$ statistic is defined as follows:
$$
f^2 = \frac{R_{\mathrm{full}}^2 - R_{\mathrm{reduced}}^2}{1 - R_{\mathrm{full}}^2}
$$
where the formula refers to the $R^2$ fit statistics for the two models being compared. As an example, imagine we already know that GPA in college has a significant relationship with job performance, and we wish to determine if our proposed screening test had incremental validity on top of knowing college GPA. We might run two linear regression models, one relating job performance to GPA, and another relating job performance to *both* GPA and screening test score. Assuming we would observe a relatively small effect size for our screening test, we assume $f^2 = 0.05$^[Cohen's rule of thumb for $f^2$ effect sizes is Small: 0.02, Medium: 0.15, Large: 0.35.], we can plot sample size against power in determining whether the two models are significantly different. We will also need to define the number of predictors in the full model (`p1 = 2`) and the reduced model (`p2 = 1`). The plot is shown in Figure \@ref(fig:regressionpower).
```{r regressionpower, fig.cap="Plot of power against sample size for a small effect of a second input variable in a linear regression model", fig.align="center", out.width = if(knitr::is_latex_output()) {"90%"}}
sample_sizes <- 100:300
f_sq_power <- WebPower::wp.regression(n = sample_sizes,
p1 = 2, p2 = 1, f2 = 0.05)
plot(f_sq_power)
```
`r if (knitr::is_latex_output()) '\\newpage'`
## Power analysis for log-likelihood regression models
In Chapter \@ref(bin-log-reg), we reviewed how measures of fit for log-likelihood models are still the subject of some debate. Given this, it is unsurprising that measures of effect size for log-likelihood models are not well established. The most well-developed current method appeared in @demidenko, and works when we want to do a power test on a single input variable $x$ using the Wald test on the significance of model coefficients (see Section \@ref(wald) for a reminder of the Wald test).
In this method, the statistical power of a significance test on the input variable $x$ is determined using multiple inputs as follows:
1. The likelihood of a positive outcome when $x = 0$ is used to determine the intercept (`p0` in the code below).
2. The likelihood of a positive outcome when $x = 1$ is then used to determine the regression coefficient for $x$ (`p1` in the code below).
3. A distribution for $x$ is inputted (`family` below) and the parameters of that distribution are also entered (`parameter` below). For example, if the distribution is assumed to be normal then the mean and standard deviation would be entered as parameters.
4. This information is fed into the Wald test, and the power for specific sample sizes is calculated.
For example, let's assume that we wanted to determine if our new screening test had a significant effect on promotion likelihood by running an experiment on employees who were being considered for promotion. We assume that our screening test is scored on a percentile scale and has a mean of 53 and a standard deviation of 21. We know that approximately 50% of those being considered for promotion will be promoted, and we believe that the screening test may have a small effect whereby those who score zero would still have a 40% chance of promotion and every additional point scored would increase this chance by 0.2 percentage points. We run the `wp.logistic()` function in `WebPower` to plot a power curve for various sample sizes as in Figure \@ref(fig:logisticpower).
```{r logisticpower, fig.cap="Plot of power against sample size for a single input variable in logistic regression", fig.align="center", out.width = if(knitr::is_latex_output()) {"90%"}}
sample_sizes <- 50:2000
logistic_power <- WebPower::wp.logistic(n = sample_sizes,
p0 = 0.4, p1 = 0.402,
family = "normal",
parameter = c(53, 21))
plot(logistic_power)
```
This test suggests that we would need over 1000 individuals in our experiment in order to have at least an 80% chance of establishing the statistical significance of a true relationship between screening test score and promotion likelihood.
## Power analysis for hierarchical regression models
Power tests for explicit hierarchical models usually originate from the context of the design of clinical trials, which not only concern themselves with the entire sample size of a study but also need to determine the split of that sample between treatment and control. It is rare that power analysis would need to be conducted for hierarchical models in people analytics but the technology is available in the `WebPower` package to explore this.
*Cluster randomized trials* are trials where it is not possible to allocate individuals randomly to treatment or control groups and where entire clusters have been allocated at random instead. This creates substantial additional complexity in understanding statistical power and required sample sizes. The `wp.crt2arm()` function in `WebPower` supports power analysis on 2-arm trials (treatment and control), and the `wp.crt3arm()` function supports power analysis on 3-arm trials (Two different treatments and a control).
*Multisite randomized trials* are trials where individuals are assigned to treatment or control groups at random, but where these individuals also belong to different clusters which are important in modeling---for example, they may be members of clinical groups based on pre-existing conditions, or they may be being treated in different hospitals or outpatient facilities. Again, this makes for a substantially more complex calculation of statistical power. The `wp.mrt2arm()` and `wp.mrt3arm()` functions offer support for this.
Power tests are also available for structural equation models. This involves comparing a more 'complete' structural model to a 'subset' model where some of the coefficients from the more 'complete' model are set to zero. Such power tests can be valuable when structural models have been applied previously on responses to survey instruments and there is an intention to test alternative models in the future. They can provide information on required future survey participation and response rates in order to establish whether the improved fit can be established for the alternative models.
There are two approaches to power tests for structural equation models, using a chi square test and a root mean squared error (RMSEA) approach. Both of these methods take a substantial number of input parameters, consistent with the complexity of structural equation model parameters and the various alternatives for measuring fit of these models. The chi square test approach is implemented by the `wp.sem.chisq()` function, and the RMSEA approach is implemented by the `wp.sem.rmsea()` function in `WebPower`.
## Power analysis using Python
A limited set of resources for doing power analysis is available in the `stats.power` module of the `statsmodels` package. As an example, here is how we would conduct the power analysis for a paired $t$-test as in Section \@ref(simple-stats) above.
```{python}
import math
from statsmodels.stats.power import TTestPower
power = TTestPower()
n_test = power.solve_power(effect_size = 0.5,
power = 0.8,
alpha = 0.05)
print(math.ceil(n_test))
```
And a power curve can be constructed as in Figure \@ref(fig:pythonpowerplot).
```{python, eval = FALSE}
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
fig = TTestPower().plot_power(dep_var = 'nobs',
nobs = np.arange(20, 100),
effect_size = np.array([0.5]),
alpha = 0.05)
plt.show()
```
```{python, echo = FALSE, fig.show='hide'}
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
# matplotlib.rcParams['text.usetex'] = True
plt.cla()
fig = plt.figure()
fig = TTestPower().plot_power(dep_var = 'nobs',
nobs = np.arange(20, 100),
effect_size = np.array([0.5]),
alpha = 0.05)
plt.savefig("www/11/pythonpowerplot-1.png")
```
```{r pythonpowerplot, echo = FALSE, fig.cap="Plot of power against sample size for a paired t-test", fig.align="center", out.width = "90%"}
knitr::include_graphics("www/11/pythonpowerplot-1.png")
```
`r if (knitr::is_latex_output()) '\\index{power analysis|)}'`