For the paired design, which is traditionally used to obtain data for the paired t-test, we can calculate a standardized mean difference, Cohen’s d, using the average of the standard deviations of the two conditions. Cohen’s d for paired designs can be calculated as follows.
d_{av} =\frac{(M_1 - M_2)}{s_{av}}, \tag{1}
where s_{av} equals
s_{av}= \sqrt{\frac{1}{2}(S^2_1+S^2_2)}. \tag{2}
Now, (1) is of course an estimate of (3), the population value of Cohen’s d for paired designs, but we do not only need a point estimate, but also a confidence interval.
The following R-code can be used to obtain a 95% confidence interval for the estimate of \delta_{av}, the population mean difference standardized by using the average of the two standard deviations:
\delta_{av} = \frac{\mu_1 - \mu_2}{\sqrt{\frac{1}{2}(\sigma_1^2 + \sigma_2^2)}} , \tag{3}
This procedure uses the approximate procedure by Algina & Keselman (2003), which is also used by ESCI (Cumming, 2012; Cumming and Calin-Jageman, 2017), as Kline (2013) explains. The following steps are performed to obtain the 95% confidence interval.
- Use the obtained t-value of the paired t-test to estimate the non-centrality parameter \lambda. Steps 2 and 3 are for calculating a 95% confidence interval for the non-centrality parameter.
- Use an iterative procedure to find the non-centrality parameter of the t-distribution for which the observed t-value is the .025 quantile. This is the upper limit of the confidence interval for the non-centrality parameter.
- Use an iterative procedure to find the non-centrality parameter of the t-distribution for which the observed t-value is the .975 quantile. This is the lower limit of the confidence interval for the non-centrality parameter.
- To obtain a CI for \delta_{av} multiply the limits of the confidence interval for the non-centrality parameter by the value \sqrt{\frac{2S_D^2}{n(S_1^2+S_2^2)}} where n is the sample size, S_D^2 the variance of the difference scores , and S_1^2 and S_2^2 the variances of the two variables.
The following R-function does all the work. Note that with large potential values for the noncentrality parameter R issues warnings that “full precision has not been achieved in ‘pnt{final}'”. These warnings can be ignored ( I checked many examples against ESCI’s output), but in order to prevent them, I let the optimize function search from the observed value of t to maximally five times its value, and I have included the option to suppress warnings or not (just set warn = TRUE to get the warnings).
ci.d.av <- function(t, n, s.1, s.2, s.diff, warn = FALSE) {
if (t < 1) {
lims <- c(-5, 5)
} else {
lims <- c(-5*abs(t), 5*abs(t))
}
df = n - 1
multiplier = sqrt((2*s.diff^2) / (n*(s.1^2 + s.2^2)) )
loss <- function(x, prob) (pt(t, df, x) - prob)^2
if (warn == FALSE) {
ul <- suppressWarnings(optimize(loss, lims, prob=.025))$minimum
ll <- suppressWarnings(optimize(loss, lims, prob=.975))$minimum
} else {
ul <- optimize(loss, lims, prob=.025)$minimum
ll <- optimize(loss, lims, prob=.975)$minimum
}
return(round(c(ll, ul), 4)*multiplier)
}
The arguments of the function are t, the t-value of the t-test testing the hypothesis of equal population means, the sample size (n), the standard deviations (s.1 and s.2) of the two variables and the standard deviation of the difference scores (s.diff).
Calculating Cohen’s d for paired designs: an example
Here is a quick example.
library(MASS)
set.seed(1234)
# generate random multivariate normal data with sample size n = 20
theData <- mvrnorm(20, c(.5, .0), matrix(c(1, .8, .8, 1), 2, 2))
# calculate the standard deviations
sds <- apply(theData, 2, sd)
# calculate the standard deviation of difference scores
sDiff <- sd(theData[,1] - theData[,2])
# get t.value and a value for d.av
# here I use the output of the t-test in R to obtain t and the mean
# difference score (needed for calculating d.av)
theTest <- t.test(theData[,1], theData[,2], paired=TRUE)
t = theTest$statistic
d.av = theTest$estimate / mean(sds)
ci.d.av(t = t, n = 20, s.1 = sds[1], s.2 = sds[2], sDiff)
The results are that the estimate equals d_{av} = 0.87, 95\% \text{CI} [0.51, 1.22].
Alternatively, we can make use of the conf.limits.nct function of the MBESS package (Kelley, 2007a, 2007b), and proceed as follows (using the data generated above).
library(MBESS)
ci.d.av.2 <- function(t, n, s.1, s.2, s.diff) {
df = n - 1
multiplier = sqrt((2*s.diff^2) / (n*(s.1^2 + s.2^2)) )
unlist(conf.limits.nct(t, df)[c(1,3)])*multiplier
}
ci.d.av.2(t = t, n = 20, s.1 = sds[1], s.2 = sds[2], s.diff = sDiff).
References
Algina, J. & Keselman, H. J. (2003). Approximate confidence intervals for effect sizes. Educational and Psychological Measurement, 63, 721-734.
Cumming, G. (2012). Understanding the New Statistics. Effect Sizes, Confidence Interval, and Meta-Analysis. New York: Routledge.
Cumming, G. & Calin-Jageman, R. (2017). Introduction fo the New Statistics. Estimation, Open Science, and Beyond. New York: Routledge.
Kelley, K. (2007a). Methods for the Behavioral, Educational, and Social Sciences: An R Package. Behavior Research Methods, 39, 979–984.
Kline, R.b. (2013). Beyond Significance Testing. Statistics Reform in the Behavioral Sciences. (Second Edition). Washington: APA.