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 
  • $\begingroup$ There are three "approaches" to this: (1) Use 'power and sample size' procedure in statistical software (or if you trust the site, an online calculator). (2) Simulation, which you attempt in your Question. (3) Use of non-central t distribution, where the non-centrality parameter depends on the size of difference you want to detect. // Some of the Q*As in the right margin may be relevant. Or search this site for something like power of t test, non-centrality// Not clear which of the three approaches you want to learn and which are only for verification. $\endgroup$
    – BruceET
    Nov 13 '20 at 21:02

I'm cheating here... I just looked up the code for pwr.t.test and I extracted the relevant parts to produce the power of a paired, two-sided t-test:

Your input:

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)

cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))

pwr.t.test(n = length(aP), d = cohensD, type = "paired", alternative = "two.sided", sig.level= 0.05)
# power = 0.4210006

To manually reproduce:

n <- length(aP)
tsample <- 1 # 1 because type is paired
tside <- 2
sig.level <- 0.05
d <- cohensD

nu <- (n - 1) * tsample
qu <- qt(sig.level/tside, nu, lower = FALSE)
pt(qu, nu, ncp = sqrt(n/tsample) * d, lower = FALSE) +
    pt(-qu, nu, ncp = sqrt(n/tsample) * d, lower = TRUE)

# [1] 0.4210006

EDIT Here's an annotated version of the code above:

We want to calculate the power of a paired t-test with given type 1 error $\alpha = 0.05$ and effect size (as Cohen's d) determined by the sample pairs aP, aM. So the input is:

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)

sig.level <- 0.05
cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))

First, we need to find the critical value of the t-statistics that incorrectly accepts the null hypothesis in 5% of the cases. Since the test is two-sided, this means finding the values of $x$ that define the two shaded tails in the probability density function pictured below, each shaded area being 2.5% of the total area:

enter image description here

For this we can use the quantile function qt with $n - 1$ degrees of freedom:

df <- (length(aP) - 1)
qu <- qt(sig.level/2, df, lower = FALSE)

# Code for plot
x <- seq(-6, 6, length.out= 100)
y <- dt(x, df= df)
plot(x, y, type= 'l', lwd= 1.5, xlab= 'Value of T', ylab= 'Density')
polygon(c(x[x > qu], qu), c(y[x > qu], 0), col= "grey", border= 'black')
polygon(c(x[x < -qu], -qu), c(y[x < -qu], 0), col= "grey", border= 'black')

We can verify that the critical value qu (and -qu) defines 2.5% of the area by integrating the PDF between -Inf and -qu and between qu and Inf:

integrate(dt, -Inf, -qu, df= df) # -> 0.025 with absolute error < 6.1e-05
integrate(dt, qu, Inf, df= df) # -> 0.025 with absolute error < 6.1e-05

Now we assume the null hypothesis is false and difference between means is not zero but has the desired Cohen's d. So we are looking at the t-distribution with non-centrality parameter that makes it skewed in the direction of the effect size. This is how the R documetatiion describes the NCP:

The most used applications are power calculations for t-tests: Let T= (mX - m0) / (S/sqrt(n)) where mX is the ‘mean’ and S the sample standard deviation (‘sd’) of X_1, X_2, ..., X_n which are i.i.d. N(mu, sigma^2) Then T is distributed as non-central t with ‘df’= n - 1 degrees of freedom and non-centrality parameter ‘ncp’ = (mu - m0) * sqrt(n)/sigma.

So we have:

ncp <- sqrt(length(aP)) * cohensD

We want to know the percentage area of the t-distribution with this NCP and degrees of freedom that falls outside the critical values -qu and qu from above. I.e. we want the shaded areas below (the area on the right tail is practically invisible):

enter image description here

right <- pt(qu, df, ncp = ncp, lower = FALSE)
left <- pt(-qu, df, ncp = ncp, lower = TRUE)
right + left
[1] 0.42 # As per pwr.t.test()

# Code for plot
x <- seq(-12, 5, length.out= 200)
y <- dt(x, df= df, ncp= ncp)
plot(x, y, type= 'l', lwd= 1.5, xlab= 'Value of T', ylab= 'Density')
polygon(c(x[x > qu], qu), c(y[x > qu], 0), col= "grey", border= 'black')
polygon(c(x[x < -qu], -qu), c(y[x < -qu], 0), col= "grey", border= 'black')
abline(v= c(-qu, qu), lty= 'dashed', col= 'blue')

Again we can verify by integrating the PDF:

integrate(dt, -Inf, -qu, df= df, ncp= ncp) # -> 0.42 with absolute error < 1.3e-05
integrate(dt, qu, Inf, df= df, ncp= ncp) # -> 6.9e-05 with absolute error < 2.8e-08

Hope this helps (and check it's correct)!

  • 2
    $\begingroup$ It would presumably help the OP if you would unpack what the code is doing & what that means. Can you use that to clear up the OP's confusion & any mistaken assumptions they're making? $\endgroup$ Nov 13 '20 at 14:14
  • 1
    $\begingroup$ @gung-ReinstateMonica I agree, it's not a great answer. But since the OP is specifically asking for the R code to manually reproduce power, I thought this is better than nothing (I'll try to find some time to explain). $\endgroup$
    – dariober
    Nov 13 '20 at 14:25
  • 2
    $\begingroup$ I didn't say it's a bad answer. I think it's the start of a good answer. However, although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics (machine learning, etc.) not code. It can be good to provide code as well, but you should try to elaborate your substantive answer in text for people who don't read R well enough to recognize & extract the answer from the code. $\endgroup$ Nov 13 '20 at 14:30
  • 1
    $\begingroup$ @gung-ReinstateMonica - I edited my answer with some explanation $\endgroup$
    – dariober
    Nov 16 '20 at 10:31

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

Difference    Size     Power
         2      25  0.944343

enter image description here

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.$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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