Calculate the power of a paired t-test in R manually 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 
```

 A: Remember that a paired t test is a one-sample test on
differences $D_i = X_i-Y_i,$ for $i=1,2, \dots, n$ and
$D_i$ are independently $\mathsf{Norm}(\mu_D, \sigma_D).$
Consider a test of $H_0:\mu=0$ vs. $H_a:\mu > 0$ at the 5% level with $n = 25.$ You seek the power of the test against the specific alternative $\mu = \mu_a = 2 > 0.$
In order to find the power, you need to have an educated guess of the value of $\sigma.$ With $\alpha = 0.05, n = 25, \sigma = 3,$ it is possible to find $P(\mathrm{Rej\;} H_0\,|\, \mu=\mu_a).$ [Of course, if you
knew the exact value of $\sigma,$ then you would be doing a z-test instead of a t-test.]
Minitab software: Here is relevant output from a recent release of Minitab. [R and other statistical software programs have similar procedures. @dariober's Answer (+1) gives a brief mention of that--for a two-tailed test.]
The power for the specified parameters is $\pi = 0.944.$ [The probability of Type II error is $\beta = 1 - \pi = 0.065.]$
Power and Sample Size 

1-Sample t Test

Testing mean = null (versus > null)
Calculating power for mean = null + difference
α = 0.05  Assumed standard deviation = 3


            Sample
Difference    Size     Power
         2      25  0.944343


Simulation. With 100,000 iterations, we can
anticipate about two place accuracy. The approximate
result from the following simulation in R is $\pi = 0.945.$
set.seed(2020)
pv = replicate(10^5, t.test(
         rnorm(25, 2, 3), alt="g")$p.val)
mean(pv <= 0.05)
[1] 0.9449

Using non-central t distribution.
The critical value for a (one-sided) test of $H_0: \mu = 0$ vs. $H_a:\mu > 0$ at the 5% level with $n = 25$ is
$c = 1.7109.$ That is, we reject $H_0$ if
$T_0 = \frac{\bar D - 0}{S_D.\sqrt{n}} \ge c.$
c = qt(.95, 24);  c
[1] 1.710882

We seek $P\left(T_a=\frac{\bar D - \mu_a}{S_D/\sqrt{n}}\ge c\right) = 0.9443,$ where $T_a$ has a noncentral t distribution
with degrees of freedom $\nu = n-1 = 24$ and noncentrality
parameter $\delta = \sqrt{n}(2)/3 = 10/3.$ [Notice that the third paramter of the R CDF function df is the non-centrality parameter.]
del = 5(2)/3
1 - pt(c, 24, del)
[1] 0.9443429

