Authors: Satoshi Aoki
Categories: Article, Cohen's d, Effect size, Hedges' d, Heteroscedasticity, Mathematics, Psychology, Standardized mean difference
Source: Heliyon
Effect sizes of the difference, or standardized mean differences, are widely used for meta-analysis or power-analysis. However, common effect sizes of the difference such as Cohen's d or Hedges' d assume variance equality that is fragile and is often violated in practical applications. Based on Welch's t tests, we defined a new effect size of the difference between means, which did not assume variance equality, thereby providing a more accurate value for data with unequal variance. In addition, we presented the unbiased estimator of an effect size of the difference between a mean and a known constant. An R package is also provided to compute these effect sizes with their variance and confidence interval.
Keywords: Psychology, Mathematics, Effect size, Standardized mean difference, Cohen's d, Hedges' d, Heteroscedasticity
Effect sizes of the difference or, more precisely, standardized mean differences between two groups, are widely used to estimate the magnitude of effect independent of the sample size [1], to conduct meta-analysis [2], or to conduct power-analysis [3]. The American Educational Research Association (AERA) or the American Psychological Association (APA) strongly recommend effect sizes are reported in the corresponding fields [4], [5]. Furthermore, the misuse and misunderstanding of p-value have become public [6], and use of effect sizes is spreading beyond pedagogy and psychology, where effect sizes have developed, into areas such as in biology [1]. In spite of such importance, the classical effect sizes of the difference assume variance equality (homoscedasticity), which is hard to assume practically or is even expected to be violated a priori in clinical data [7]. While Bonett [8] defined a confidence interval of an effect size estimator which did not assume homoscedasticity, its parameter was not defined. This problem of variance inequality (heteroscedasticity) has been long debated [9], [10]. In addition, the unbiased estimator of an effect size of the difference between a mean and a constant was undefined. To solve these problems, based on Welch's t test [11], [12], we defined an effect size of the difference between means that does not assume homoscedasticity and calculated the unbiased estimator of an effect size of the difference between a mean and a constant.
Effect size of the difference was developed by Cohen [13], who studied in the field of psychology. Cohen [3], [13] defined the effect size as a parameter for two independently and normally distributed populations, Y1∼N(μ1,σ2) and Y2∼N(μ2,σ2):
which is expressed as d in the original articles [3], [13]. Note that both populations share the common variance σ2. The estimator of this parameter was represented as ds in [3]. However, we refer to this estimating statistic as g to distinguish it from the other d we introduce later. The statistic g is defined as
where
and, for i=1,2,
Here, Y¯1, Yj1, and n1 are the mean of the sample, the sample (random variable), and the sample size of group 1, respectively, while Y¯2, Yj2, and n2 are those of group 2. For the denominator, this effect size uses the pooled standard deviation, which suggests the most precise population variance under the assumption of equal variance [14].
In the field of pedagogy, Glass [2] suggested another effect size of the difference, independently of Cohen's works. He defined it as “the mean difference on the outcome variable between treated and untreated subjects divided by the within group standard deviation,” where “the within group standard deviation” corresponds to the standard deviation of the untreated group. He clearly distinguished the treated (experimental) group from the untreated (control) group, and there was no assumption regarding the two groups. His effect size was subsequently formulated and named Glass' Δ by Hedges [14], which is
where YE¯ is the mean of the variable in the experimental group, YC¯ is that in the control group, and SC is the unbiased standard deviation of the control group.
Hedges [14] also defined the δ (1) and the g (2) independently of Cohen. Furthermore, Hedges [14] indicated that g (2) is biased from δ (1), making it unsuitable for analyses that do not treat the entire population. The unbiased estimator of δ (1) is defined as gU in [14] and d in [15]. In this study, we call it d, which is
Using the gamma function, the correction coefficient J is defined as
The effect sizes g (2) and d (5) are used in various fields of science, but they assume homoscedasticity just like Student's t-test [16], [17]. When this assumption of homoscedasticity is violated, Grissom [9] recommended the use of Glass's Δ (4) instead of d (5). However, Glass's Δ (4) and d (5) have different meaning because of the difference in denominator. Therefore, Glass's Δ (4) cannot substitute for d (5) in a strict sense. Behavior of g (2), Δ (4), and d (5) under heteroscedasticity was studied in [10], although the justification for using effect size parameter α, that they defined, to measure the statistic bias under heteroscedasticity was not shown.
Bonett [8] in psychology proposed a confidence interval (CI) of effect size which does not assume homoscedasticity. First, he defined a general effect size estimator
where ∑j=1kcj=0, Y¯j is a sample mean, and s=k−1∑j=1ksj2. Concerning effect size of the difference between two means, substituting k=2, c1=1 and c2=−1 gives
Then, he assumed its corresponding parameter and its CI. The CI was calculated using approximation of CI [18] and variance of the estimator which was approximately calculated without assuming homoscedasticity. The parameters estimated by δˆ (7) or (8) were not formulated. Namely, he defined the CI for heteroscedasticity without defining a parameter, and this can be a problem. When the estimator does not always correspond to a single parameter, the CI of an undefined parameter loses its consistency in what to estimate, and heteroscedasticity or difference of sample sizes can change the correspondence between an estimator and a parameter (see section 5.2). Although his CI was effective relative to the other CIs in his simulation experiment where the parameter was given a value, what the value meant could change depending on the variance and sample size, and the change could not be expected since the parameter was not formulated.
It should be noted, Cohen [3] also defined a parameter of an effect size of the difference between a mean and a constant for a normally distributed population N1(μ,σ12) and a known constant C as
Cohen [3] originally referred to this as d3′, but we refer to this as γ (9) to clearly distinguish it from d (5). Cohen [3] also defined a biased estimator of an effect size for a normally distributed population with the sample value Yi1 (i=1,...,n1), the sample mean Y1¯, and a known constant C as
The s1 is the square root of (3). Cohen [3] originally referred to this as ds′, but we refer this to cbiased for the reason described above. To the best of my knowledge, the unbiased estimator of γ (9) has not been shown.
There are other effect sizes of the difference that do not assume normality or independence. Since their assumption is different from that of effect size we focus on, we do not treat them in detail and briefly introduce them. Dunlap et al. [19] invented effect size of the difference between two correlated paired groups. Algina et al. [20] proposed robust effect size of the difference, which is based on g (2) using 20% trimmed mean and 20% Winsorized variance assuming that samples are taken from an observing population and another contaminating population.
First, we define the parameter of an effect size of the difference between means for two independently and normally distributed populations N1(μ1,σ12) and N2(μ2,σ22) as
where r is a non-negative real number. This parameter is not generalization of δ (1) and is different from it. Then, suppose two independently and normally distributed populations with the samples Yi1 (i=1,...,n1) and Yi2 (i=1,...,n2), and the sample mean Y¯1 and Y¯2. Based on the statistic tw, the so-called Welch's t [11], [12], a biased estimator of ϵr (11) is defined as
where
si2 is the same as (3), and
Finally, e, the unbiased estimator of ϵr (11), is
Therefore,
Here, r corresponds to the ratio n1/n2. J is the correction coefficient that is defined in equation (6). The degree of freedom f is approximately calculated using the Welch-Satterthwaite equation [11], [21] as
The variance of e (15) is
Although this effect size is derived from the difference, we refer to it as e not d. This is because Cohen's d (2) and Hedges' d (5) already exist, and more d would cause further confusion. The proof of the bias correction and variance derivation does not assume homoscedasticity (see the Appendix). In addition, e (15) is a consistent estimator of ϵr (11) at the same time. See the Appendix for the proof of the consistency.
Using cbiased (10), the unbiased estimator of the effect size parameter γ (9) is defined for a normally distributed population with the sample value Yi1 (i=1,...,n1), the sample mean Y1¯, and a known constant C as
Therefore,
The correction coefficient J (6) is the same as the one used above. The variance of c is
See the Appendix for proofs of the bias correction and the derivation of the variance. In addition, c (17) is a consistent estimator of γ (9) (see the Appendix for the proof). When interested in constants rather than variables, c′ defined as
can be used instead of c.
In terms of the effect sizes of the difference, the CI based on a noncentral t variate is not directly given by a formula [22]. The CI is derived from that of noncentral parameters of noncentral t-distribution, which is in turn obtained by some searching method. The CI based on the biased effect sizes are given
and
where ncpL is the noncentral parameter that gives the upper limit of cumulative probability (e.g., 0.975 cumulative probability for 95% CI) for noncentral t-distribution with the corresponding t value (see the discussion section) and the degree of freedom, and ncpH is that which gives the lower limit (e.g., 0.025 cumulative probability for 95% CI), and n˜ and n1 are the same as (14) and (10). The CIs based on the unbiased estimator of the effect sizes are given by multiplying the corresponding correction coefficient J (6) of the corresponding degree of freedom to the above intervals.
The CI by Bonett [8] is calculated using variance of the estimator which is approximately calculated without assuming homoscedasticity and approximate assumption of CI [18]. Therefore, it is not necessary to apply Bonett's CI to e (15) or c (17), because the derivation of their CIs does not assume homoscedasticity, and their exact CIs can be calculated without approximation.
I developed a new package es.dif for R [23]. It enables the statistics d (5), e (15), c (17), their biased statistics, variance, and CI based on the two samples or their mean, variance, and sample size to be computed. In this package, approximation of J (6) [14] is not employed unless its degree of freedom exceeds 342, when the gamma function returns values that are too large to be treated in R. The CI is obtained by binary search.
The remainder of this section presents some examples of the package. First, the following script calculates d (5), e (15), their variances and 95% CIs for data 1 (0,1,2,3,4) and data 2 (0,0,1,2,2).
library(es.dif) > data1<-c(0,1,2,3,4) > data2<-c(0,0,1,2,2) > es.d(data1,data2) [,1] [,2] [1,] "Hedges' " "0.682379579593354" [2,] "variance:" "0.484026380702367" [3,] "CI:" "[ -0.503527216375147 , 1.82938058482178 ]" > es.e(data1,data2) [,1] [,2] [1,] "Unbiased " "0.668264936033828" [2,] "variance:" "0.506830833214916" [3,] "CI:" "[ -0.50334965496395 , 1.7965317007171 ]"
Using options of the function, you can change the type I error rate for the CI, calculate biased effect sizes, and output results in the vector style. For example, cbiased (10) with 99% CI in the vector style is calculated by this script.
library(es.dif) > data1<-c(0,0,1,2,2) > data2<-c(2) > es.c(data1,data2,alpha=0.01, unbiased=FALSE,vector_out=TRUE) [1] -1.0000000 0.9292037 -2.5390625 0.5778885
In the vector-style output, the four values in the vector show the effect size, its variance, and lower and higher limits of the CI. In addition, this package includes functions that can output effect sizes from the (estimated) parameters and the sample sizes. The following scripts compute d (5) and e (15) for two populations, N(1,2) and N(0,1) with the sample size 5 and 10, respectively.
library(es.dif) > mean1<-1 > mean2<-0 > var1<-2 > var2<-1 > n1<-5 > n2<-10 > es.para.d(mean1,mean2,var1,var2,n1,n2) [,1] [,2] [1,] "Hedges' " "0.82286529714397" [2,] "variance:" "0.349443397657368" [3,] "CI:" "[ -0.248827687382689 , 1.86616833367494 ]" > es.para.e(mean1,mean2,var1,var2,n1,n2) [,1] [,2] [1,] "Unbiased " "0.674259756444758" [2,] "variance:" "0.41613476136966" [3,] "CI:" "[ -0.354146439977423 , 1.65626025590509 ]"
These types of functions also have the options for the type I error rate, the biased effect size, and the vector-style output.
While the situation to use c (17) is clearly different, the e (15) and d (5) have a similar application range in practice. Therefore, we prepared an example of the applications in which the sample variances are not equal. Table 1 shows well-known data of three Iris species by Fisher [24], which can also be checked in R [23] using a command “iris”. Note that only the petal width of I. setosa has fewer significant digits. For this data, we calculated d (5), e (15), the ratio of d (5) to e (15), and the ratio of the standard deviations of the two comparing data. Theoretically, e (15) is a more precise estimator of its own parameter than d (5) under this heteroscedasticity.
The calculated result is shown in Table 2. When considering their significant digits, the comparing pair of the sepal length of I. setosa and I. virginica showed the different effect size of d (5) and e (15) (in bold in Table 2). Even though most pairs showed identical values of d (5) and e (15), the result revealed that violation of the assumption of homoscedasticity in d (5) can affect the result even in two significant digits.
Fig. 1 shows the ratio of d (5) to e (15) plotted against the ratio of standard deviations of the comparing data. This figure shows that the similar two standard deviations give similar d (5) and e (15). In other words, the more different two standard deviations encourage the use of e (15) over d (5) more.
Figure 1 Plotted graph of Table 2.
To examine the nature of d (5) and e (15), we also conducted a simulation study. In addition to d (5) and e (15), Bonett's statistic δˆ (8) was also included as a reference, although its accuracy cannot be discussed because of the lack of the parameter definition. The above effect sizes and their width of 95% CI were calculated for 100,000 Monte Calro replications from N(1,σ12) and N(0,σ22) for each condition, and they were represented by their average values. The population means were fixed to 1 and 0. The sample sizes were changed from 10 to 30 by 10. The population standard deviation σ1 was fixed to 1 and σ2 was changed 1 to 10 by 1. However, some redundant data were omitted from the result. The calculation was conducted using es.dif R package shown above and metafor R package [25]. The R source code used for the simulation was shown in the Appendix.
Table 3 shows the result of the simulation. When the sample size ratio was conserved under σ1≠σ2, e (15) gave more similar and concordant values than d (5). For example, e (15) for n1=n2=10,20,30 under σ2=10 were 0.142, 0.140 and 0.141, whereas the corresponding d (5) were 0.148, 0.143 and 0.143. This is the nature and advantage of e (15) which is designed to estimate the same parameter under heteroscedasticity and the same sample size ratio. The width of CI was narrowest for d (5) under σ1=σ2, and e (15) had the second narrowest. Under σ1≠σ2, e (15), e (15) and δˆ (8) had the narrowest CI under n1=n2, n1>n2 and n1<n2, respectively. The narrowest CIs of e (15) were followed by d (5), whereas what followed the narrowest CIs of δˆ (8) was not fixed. It was shown that e (15) had wider situation under which it had the narrowest or second narrowest CI than d (5) or δˆ (8). Bonett's statistic δˆ (8) equaled to d (5) under n1=n2 as their definition. Under n1≠n2 and σ1≠σ2, e (15) was closer to δˆ (8) than d (5). This might imply relative accuracy of δˆ (8) over d (5) under heteroscedasticity.
Comparison of t tests and the effect sizes of the difference except δˆ (8) shows the clear correspondence between them (Table 4). Statistic d (5) corresponds to the unpaired two-sample t test [16], [17], whose statistic is the basis of g (2). Statistic ebiased (12) uses the statistic (13) of Welch's t test [12], which aims to test two means with unequal variances, and cbiased (10) uses the same statistic as the one-sample t test [17]. Considering this, it is natural that power analyses should be conducted, using the corresponding pair of the effect size and t test. In other words, power analyses of Student's one-sample t test, Student's unpaired two-sample t test, and Welch's t test should be conducted based on the c statistic (17), d (5), and the e statistic (15), respectively. Co-use of non-corresponding t-test and effect size causes inconsistence of the assumption about the population(s).
In this subsection, the relationship between the effect sizes of the difference and sample sizes is described. The value of g (2), a biased estimator of the effect size of the difference under homoscedasticity, is independent of the sample sizes when the assumption of homoscedasticity (s1=s2) is fulfilled. When s1≠s2, it depends on the ratio q=(n1−1)/(n2−1) as implied in [9]. This is because g (2) is no longer an estimator of δ (1) under s1≠s2, and it will be a biased estimator of the other parameter δq′, which is
Note that even d (5) cannot be the unbiased estimator of δq′ when s1≠s2, because g (2) is not distributed as non-central t variate in this situation. Even if n1 and n2 vary, g (2) roughly estimates the same parameter, given the ratio q is fixed.
Next, the ebiased (12) is a biased estimator of ϵr (11), but ϵr (11) equals to the other parameters in the particular situation. When s1=s2, ϵr=δ, and ebiased (12) equals to g (2), and is independent of the sample sizes. When s1≠s2 and n1=n2, ϵr=δq′. In this case, ebiased (12) equals to g (2) and is also independent of the sample sizes. While d (5) is a biased estimator of δq′, e (15) is its unbiased estimator. Therefore, usage of e (15) is always preferable to d (5) in this situation. When s1≠s2 and n1≠n2, ebiased (12) depends on the rate r=n1/n2. Therefore, strictly speaking, multiple ebiaseds can be comparable only when the sample size ratio r is identical.
The effect size estimator δˆ (8) did not have a defined parameter, but when n1=n2 and s1=s2, δˆ (8) equals to g (2) and ebiased (12), and is independent of sample size. Under n1=n2 and s1≠s2, δˆ (8) also equals to g (2) and suffers from the same problem as it. Under n1≠n2 and s1≠s2, the value of δˆ (8) is no longer the same as g (2), and precise discussion on its behavior is hinderred by the lack of its parameter definition. When trying to consider δˆ (8) as a noncentral t-variate like the other effect sizes, its degree of freedom should be about n1+n2−2, and n1 and n2 should affect the degree of freedom under s1≠s2.
Unlike g (2) or ebiased (12), cbiased (10) is always independent of the sample size.
The behavior of the unbiased estimator of the effect sizes (d (5), e (15), and c (17)) are almost identical to those that are biased, but they slightly increase as the sample sizes become large. This is because of the correction coefficient J (6), and its behavior is illustrated in detail in [14].
In summary, in terms of the effect size of the difference between two means, usage of e (15) is preferable to d (5) or δˆ (8), and e (15) can be the remedy for application of effect size of the difference under heteroscedasticity. However, when the ratio of the two sample sizes cannot be set as uniform under heteroscedasticity, neither d (5) nor e (15) can be precisely compared. This is a form of the Behrens-Fisher problem, which cannot be solved strictly.
The effect size e (15) has a vast applicable range covering all kinds of natural and social sciences. This is because e (15) corresponds to Welch's t test, whose use is nowadays encouraged over Student's t test (e.g., [26]). The effect size e (15) is the best option, especially when the ratio of the sample sizes of two groups can be fixed. The effect size c (17) has a relatively narrower range regarding the application. In comparison of paired two groups (the difference in pairs vs. 0) and in some simulation studies (result of simulation vs. the optimal value) or physics (result of experiment vs. physical constant), an effect size of the constant may be needed.
Satoshi Aoki: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.
This study was partly supported by National Bioresource Project from AMED, grant number 16km0210053j0005.
The authors declare no conflict of interest.
Supplementary content related to this article has been published online at https://CRAN.R-project.org/package=es.dif.
I would like to thank Prof. Motomi Ito and Dr. Fumio Tajima for their contribution to the manuscript and Dr. Masakazu Shimada for highlighting a miscalculation in the manuscript.
In summary, this proof is an application of the proof in [14] to the statistic v in [11]. Suppose two independently and normally distributed populations N1(μ1,σ12) and N2(μ2,σ22). Their sample means are Y¯1 and Y¯2, and their samples are Yi1 (i=1,...,n1) and Yi2 (i=1,...,n2). The statistic ebiased (12) between them can be converted into
where
Here, since N1 and N2 are independently and normally distributed, the numerator of (18) has the normal distribution of N(θ,1), where
and the si2 is the same as (3). In the denominator, wf is approximately distributed as χ2(f) [11]. Therefore, n˜ebiased is distributed as a non-central t variate with the non-centrality parameter θ and approximate degree of freedom f (16). From the nature of the non-central t distribution [27], the expected value of ebiased (12) is
Now, assuming r=n1/n2, then θ/n˜=ϵr. In this case, the expected value of e (15) is
Thus, e (15) is an unbiased estimator of ϵr (11). The variation of ebised (12) is
Therefore, the variation of e (15) is
The bias correction and derivation of the variance can be proved in the same way as that of d (5). The statistic cbiased (10) can be converted into
and this (19) is distributed as a non-central t variate with non-centrality parameter
and degree of freedom n1−1. Therefore, the expected value cbiased (10) is
Because c=cbiasedJ(n1−1), the expected value of c (17) is
Thus, c is an unbiased estimator of the effect size parameter γ (9). The variation of cbiased (10) is
Therefore, the variation of c (17) is
First, we treat the proof of c, which is simpler than that of e. For the proof, we introduce a lemma.
Now, we move on to the proof of c (17). When n1→∞, the variance of c will be
Thus, limn1→∞var(c)→0, and c is an unbiased estimator of γ. Therefore, based on Lemma 1, c (17) is a consistent estimator of γ (9). □
On the other hand, e (15) consists of two populations. Therefore, a variation of the previous lemma is necessary.
This lemma can be proved in the same way as Lemma 1.
Now, consider n1=rϕ and n2=ϕ to assume ϕ→∞, which equals to (n1,n2)→(∞,∞). Note that r>0 and θ>0, since n1≥1 and n2≥1. Using r and ϕ, f (6) and n˜ (14) can be expressed as
and
Therefore, when ϕ→∞, the variance of e (15) will be
The limit does not contain r, meaning lim(n1,n2)→(∞,∞)var(e) always gives an identical value 0. Also, e is an unbiased estimator of ϵr (11). Therefore, based on Lemma 2, e (15) is a consistent estimator of ϵr (11). □
The simulation study in this article was conducted using the following code in R: