# Power analysis for minimum effect tests and good enough range hypotheses

• 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?

## 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