1
$\begingroup$
  • Power analyses allow one to determine the minimum sample size required to detect a certain effect (when there is one)
  • To calculate power, one determines the minimum effect one wants to detect
  • With NHST, I can set this with library(simr), and calculate the lowest sample size needed to find a significant result to reject the null hypothesis that there is no effect whatsoever
  • The problem with point-null hypotheses is that they are (quasi-)always false (Cohen, 1990)
  • A better way is to test "good enough" range hypotheses (Serlin and Lapsley, 1985) or interval or minimum effect tests (Lakens, 2022 Chapter 9, see figure)
    • From Rao and Lovric (2016): Rather than testing problematic and false point-null hypotheses ($H_0: \theta=\theta_0$ and $H_1: \theta≠\theta_0$), we should test negligible null hypotheses ($H_0: |\theta-\theta_0|≤\delta$, effect size is negligible; and $H_1: |\theta-\theta_0|>\delta$, effect size is practically meaningful)
  • By doing so, our predictions will always be either be supported or rejected when the sample is large enough
  • With a minimum effect hypothesis, I have to specify a minimum effect
    • therefore I need more power to find a significant result
    • However, the minimum effect I am interested in detecting is also my power value, which basically cancel each other out, so I would need infinite power
      • How to deal with this?

Lakens 2022 Figure 9.1

Example

  • I have a new fertilizer
  • I set up a lab experiment with 10 seedlings with the old fertilizer ("control") and 10 with the new one ("treatment")
  • I measure seedling dry weight after a pre-defined period
  • I expect dry weight in the control group to be 100 mg (SD 3), and the new fertilizer to be more than 5% more effective (105 mg, same SD)
    • As such, my test hypothesis is $H_{test}: |\theta_{trt}-\theta_{ctrl}|>5$, and alternative hypothesis $H_{alt}: |\theta_{trt}-\theta_{ctrl}|≤5$
  • I want at least 90% power at alpha = .05
treatment <- rep(c("control", "treatment"), each = 10)

weight <- c(rnorm(10, mean = 100, sd = 3),
        rnorm(10, mean = 105, sd = 3))

data_weight <- data.frame(treatment, weight)

library(ggplot2)

ggplot(data_weight, aes(x = treatment, y = weight)) +
  geom_point()

Conventional NHST power

  • With conventional NHST, I want to be able to detect a statistically significant difference either way
  • library(simr) allows me to create a linear model object and run a power analysis
  • To have at least 90% power at alpha = .05, the power analysis demonstrates at least 6 samples are needed in both groups
model_weight1 <- lm(weight ~ treatment, data = data_weight)

# Conventional NHST power
simr::powerSim(model_weight1)
> Power for predictor 'treatment', (95% confidence interval):====================|
>       99.80% (99.28, 99.98)
> 
> Test: t-test
> 
> Based on 1000 simulations, (0 warnings, 0 errors)
> alpha = 0.05, nrow = 20
> 
> Time elapsed: 0 h 0 m 5 s
> 
> nb: result might be an observed power calculation
> Warning message:
> In observedPowerWarning(sim) :
>   This appears to be an "observed power" calculation

simr::powerCurve(model_weight1, within = "treatment")
> Calculating power at 8 sample sizes within treatment
> Power for predictor 'treatment', (95% confidence interval),====================|
> by number of observations within treatment:
>       3: 50.50% (47.35, 53.64) - 6 rows
>       4: 71.30% (68.39, 74.09) - 8 rows
>       5: 84.40% (82.00, 86.60) - 10 rows
>       6: 91.90% (90.03, 93.52) - 12 rows
>       7: 96.20% (94.82, 97.30) - 14 rows
>       8: 98.20% (97.17, 98.93) - 16 rows
>       9: 99.50% (98.84, 99.84) - 18 rows
>      10: 99.90% (99.44, 100.00) - 20 rows
> 
> Time elapsed: 0 h 0 m 48 s
> Warning message:
> In observedPowerWarning(sim) :
>   This appears to be an "observed power" calculation

Power for the non-nil hypothesis

  • I specified I expect the new fertilizer to be at least 5% more effective than the old one
  • Basically, I could subtract the expectation from each observation in the treatment group, and then run the power analysis
  • However, this basically makes the mean difference between the groups equal to 0, and, therefore, power will be extremely low no matter the sample size
data_weight <- data_weight |>
    dplyr::mutate(
        weight_mod = dplyr::case_when(treatment == 'treatment' ~ weight - 5,
TRUE ~ weight)
)

model_weight2 <- lm(weight_mod ~ treatment, data = data_weight)

simr::powerSim(model_weight2)
> Power for predictor 'treatment', (95% confidence interval):====================|
>        8.20% ( 6.57, 10.08)
> 
> Test: t-test
> 
> Based on 1000 simulations, (0 warnings, 0 errors)
> alpha = 0.05, nrow = 20
> 
> Time elapsed: 0 h 0 m 6 s
> 
> nb: result might be an observed power calculation
> Warning message:
> In observedPowerWarning(sim) :
>   This appears to be an "observed power" calculation


simr::powerCurve(model_weight2, within = "treatment")
> Calculating power at 8 sample sizes within treatment
> Power for predictor 'treatment', (95% confidence interval),====================|
> by number of observations within treatment:
>       3:  5.80% ( 4.43,  7.43) - 6 rows
>       4:  4.80% ( 3.56,  6.31) - 8 rows
>       5:  5.50% ( 4.17,  7.10) - 10 rows
>       6:  4.40% ( 3.21,  5.86) - 12 rows
>       7:  5.30% ( 3.99,  6.88) - 14 rows
>       8:  7.10% ( 5.59,  8.87) - 16 rows
>       9:  7.30% ( 5.77,  9.09) - 18 rows
>      10:  7.80% ( 6.21,  9.64) - 20 rows
> 
> Time elapsed: 0 h 0 m 53 s
> Warning message:
> In observedPowerWarning(sim) :
>   This appears to be an "observed power" calculation

Example results

Three scenarios. See simulations below.

  1. Effect size is 0. P-value of the null hypothesis is nonsignificant, of the hypothesis that the effect is smaller than 5 is very small, and of the hypothesis that the effect is larger than five very large. Conclusion: support that there is no effect.
  2. Effect size is 5. P-value of the null hypothesis is significant, of the hypothesis that the effect is smaller than 5 is non-significant, and of the hypothesis that the effect is larger than five is also non-significant. Conclusion: no support that there is or is not an effect.
  3. Effect size is 10. P-value of the null hypothesis is significant, of the hypothesis that the effect is smaller than 5 is non-significant, and of the hypothesis that the effect is larger than five is significant. Conclusion: support that there is a real effect.
set.seed(102)

# create data frame
treatment <- rep(c("control", "treatment"), each = 10)

weight_100 <- c(rnorm(10, mean = 100, sd = 3),
                rnorm(10, mean = 100, sd = 3))

weight_105 <- c(rnorm(10, mean = 100, sd = 3),
            rnorm(10, mean = 105, sd = 3))

weight_110 <- c(rnorm(10, mean = 100, sd = 3),
              rnorm(10, mean = 110, sd = 3))

data_weight_results <- data.frame(treatment, weight_95, weight_105, weight_110)


# run analyses with emmeans
## weight 105
model105 <- lm(weight_105 ~ treatment, data = data_weight_results)

EMM_trt_105 <- emmeans::emmeans(model105, specs = ~ treatment)

EMM_trt_105

PRS_trt_105 <- pairs(EMM_trt_105)

emmeans::test(PRS_trt_105)
>     contrast            estimate   SE df t.ratio p.value
>     control - treatment    -7.04 1.56 18  -4.521  0.0003
emmeans::test(PRS_trt_105, null = -5, side = ">")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment    -7.04 1.56 18   -5  -1.311  0.8968
emmeans::test(PRS_trt_105, null = -5, side = "<")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment    -7.04 1.56 18   -5  -1.311  0.1032
     
## weight 100
model100 <- lm(weight_100 ~ treatment, data = data_weight_results)

EMM_trt_100 <- emmeans::emmeans(model100, specs = ~ treatment)

EMM_trt_100

PRS_trt_100 <- pairs(EMM_trt_100)

emmeans::test(PRS_trt_100)
>     contrast            estimate   SE df t.ratio p.value
>     control - treatment     2.83 1.45 18   1.956  0.0661
emmeans::test(PRS_trt_100, null = -5, side = ">")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment     2.83 1.45 18   -5   5.415  <.0001
emmeans::test(PRS_trt_100, null = -5, side = "<")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment     2.83 1.45 18   -5   5.415  1.0000

## weight 110
model110 <- lm(weight_110 ~ treatment, data = data_weight_results)

EMM_trt_110 <- emmeans::emmeans(model110, specs = ~ treatment)

EMM_trt_110

PRS_trt_110 <- pairs(EMM_trt_110)

emmeans::test(PRS_trt_110)
>     contrast            estimate   SE df t.ratio p.value
>     control - treatment    -10.4 1.41 18  -7.410  <.0001
emmeans::test(PRS_trt_110, null = -5, side = ">")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment    -10.4 1.41 18   -5  -3.854  0.9994
emmeans::test(PRS_trt_110, null = -5, side = "<")
>     contrast            estimate   SE df null t.ratio p.value
>     control - treatment    -10.4 1.41 18   -5  -3.854  0.0006
ES null < >
0 .066 <.001 1.000
5 .003 .103 .897
10 <.001 .999 .001

Table 1. P-values of simulated effect sizes

References

$\endgroup$

1 Answer 1

0
$\begingroup$

I'll try to answer my own question, for posterity. I think I can thanks to discussions with Daniël Lakens on the Fediverse.

The primary issue here is that it is impossible to verify null hypotheses (i.e. the effect size is exactly zero). If my hypothesis were ES_trt = ES_ctl, no amount of data would be able to prove this. Hence, I would need infinite power. With my hypothesis ES_trt - ESS_ctl > 5, I basically shifted the null from 0 to 5, so I'm testing a non-nil null hypothesis.

The way forward is to calculate power for an effect size of 5 (which I did and for 90% power I need 6 samples in each group), and test a conventional null hypothesis for this sample size together with an equivalence test with the same effect size. By doing so, I am (relatively) certain I can show support for or against the hypothesis, while not inflating my Type I error rate.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.