As an exercise I wanted to perform a paired t-test manually in R to refresh a lecture I had in the past. Everything went well, but then I thought about calculating the power of this paired t-test and that's where the trouble began.
I know that the power is the area under the alternative distribution minus the area of the type II error ($\beta$), which is delimited by the $\alpha$ significance level. So basically, in this example I need to find $P(X ≤ \alpha)$ of the alternative distribution that is centered around the observed mean difference I calculated, but to be frank I'm not sure how to construct that distribution. I tried to use the same procedure as for the t-statistic under the null, but that doesn't make sense, since the expected mean and the observed mean would be the same, thus the whole term would just equal 0 (1-pt((expMean - obsMean)*stdError, df
). And as far as I know, t-distributions are only used under the assumption that the null hypothesis is true. From here on I'm just getting more confused and I think that I'm missing something obvious.
I used the pwr.t.test function from the pwr package to compare my result.
It would be very helpful if somebody could help me to do such tests manually, as most solutions I found elsewhere skip the part I'm trying to do manually and simply use some sort of power calculator.
The code I used:
# data
aP <- c(0.5331039, 0.4578532, 0.3129205, 0.5144858, 0.8149759, 0.4136268)
aM <- c(0.2750040, 0.5056830, 0.4828734, 0.4439654, 0.2738658, 0.3081768)
# difference between P and M
Diff <- aM - aP
# INIT t test
obsMean <- mean(Diff)
expMean <- 0
stdError <- (sqrt(length(Diff))/sd(Diff))
n <- length(aP)
df <- n - 1
alpha = 0.05
# T-statistic
T_stat <- (obsMean-expMean)*stdError; T_stat
# critical value
crit_values <- qt(c(0.025,0.975),df) # lower bound = -2.570582
p_value <- 2*(pt(T_stat, df)); p_value
p_value < alpha
# comparison
t.test(aM, aP, paired = TRUE, alternative = "two.sided")
# INIT power
obsMean <- mean(Diff)
expMean <- mean(Diff)
# power???
power <- 1-pt((expMean - obsMean)*stdError, df); power
# comparison
cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))
pwr.t.test(n = 6,d = cohensD, type = "paired", alternative = "two.sided")
# power = 0.4210006
```