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