# Example:
<- .7 # observed proportion
p <- .5 # proportion of interest
p0 <- 50 # sample size
n
<- (p-p0) / sqrt(p*(1-p)/n)
z
<- 2*(1-pnorm(z)) # two tailed test
pval
data.frame(z,pval)
z pval
1 3.086067 0.002028231
For dichotomous relationships that involve proportions, there are many variations of effect sizes that one can use. Commonly used effect size measures for statistical procedures on categorical data include: phi coefficient (\(\phi\)), Cramer’s \(V\), Cohen’s \(h\), Cohen’s \(\omega\), odds ratio (\(OR\)), risk difference (\(RD\)), and relative risk (\(RR\)).
Type | Description | Section |
---|---|---|
\(\phi\) - phi coefficient | Pearson correlation between two binary variables (i.e., 2x2 contingency tables). | Section 9.2.1 |
\(V\) - Cramer’s V | Measures the association between categorical variables. Similar to a \(\phi\) coefficient, but meant for contingency tables larger than 2x2. | Section 9.2.2 |
\(h\) - Cohen’s h | Pearson correlation between two binary variables. Difficult to interpret. | Section 9.2.3 |
\(w\) - Cohen’s w | Association between two categorical variables and it is computed identically to the \(\phi\) coefficient. If computed on a 2x2 contingency table, it will have an identical value to \(\phi\). | Section 9.2.4 |
\(\boldsymbol{פ}\) - Ben-Shachar’s Fei | A correction to Cohen’s \(w\) for one dimensional count tables. | Section 9.2.5 |
\(OR\) - Odds Ratio | Ratio of odds of an event occurring between treatment and control groups | Section 9.2.6 |
\(RD\) - Risk Difference | Difference between proportions in treatment and control groups. | Section 9.2.7 |
\(RR\) - Relative Risk | Ratio of proportions in the treatment and control groups. | Section 9.2.8 |
If we have a single sample and we want to assess the difference between a proportion and some proportion of interest. We can first calculate the test statistic by comparing the observed proportion (\(p\)) vs the proportion of interest (\(p_0\)):
\[ z = \frac{p-p_0}{\sqrt{\frac{p(1-p)}{n}}}, \tag{9.1}\] where \(n\) is the sample size. Note that this is only valid if the proportion of interest is chance (\(p_0=.50\)) because the sampling distribution with a proportion of .50 is normal. However if the proportion of interest is not .50, then we should instead compute Cohen’s \(h\) (see Section 9.2.3 for details), which transforms the scale so that the distributions are normal regardless of the proportion. The test-statistic with Cohen’s \(h\),
\[ z = h\sqrt{n} \tag{9.2}\]
Let’s try testing the proportion against chance (\(p_0=.50\)) in R. We can then calculate the p-value in base R by using the pnorm()
function:
# Example:
<- .7 # observed proportion
p <- .5 # proportion of interest
p0 <- 50 # sample size
n
<- (p-p0) / sqrt(p*(1-p)/n)
z
<- 2*(1-pnorm(z)) # two tailed test
pval
data.frame(z,pval)
z pval
1 3.086067 0.002028231
Results show a significant difference from chance with \(z\) = 3.09 and p-val = .002
Phi coefficient (\(\phi\)) is a measure of association between two binary variables (therefore, it ONLY applies to 2 by 2 contingency tables, i.e., each variable has only two levels). It is a special case of the Pearson correlation coefficient and an \(r\) for two binary variables is equal to phi. Note that unlike \(r\) that ranges from -1 to 1, phi ranges from 0 to 1. Also, the sign of \(r\) indicates the direction of association, whereas to get the direction of an association given a 2 by 2 contingency table, we need to look at the table itself; phi only provides a measure of strength. The 2 by 2 contingency table is illustrated by Table 9.1.
\(X=0\) | \(X=1\) | |
---|---|---|
\(Y=0\) | \(n_{00}\) | \(n_{10}\) |
\(Y=1\) | \(n_{01}\) | \(n_{11}\) |
The sample sizes within each cell provide us with the necessary information to estimate the relationship between the two variables. A large phi coefficient would be expected to have relatively large sample sizes in the diagonal cells (\(n_{00}\) and \(n_{11}\)) and relatively low sample sizes in the off-diagonal cells (\(n_{01}\) and \(n_{10}\)). To calculate phi, it can be calculated from the cells of the contingency table directly (adapted from equation 1, Guilford 1965),
\[ \phi = \frac{n_{11}n_{00} -n_{10}n_{01}}{\sqrt{(n_{00} + n_{01})(n_{10} + n_{11})(n_{00} + n_{10})(n_{01} + n_{11})}} \tag{9.3}\]
or more conveniently, from the \(\chi^2\)-statistic (equation 7.2.5, Cohen 1988),
\[ \phi = \sqrt{\frac{\chi^2}{n}} \tag{9.4}\]
Where \(n\) is the total sample size (i.e., the sum of all the cells). Using the effectsize
package in R, we can calculate the the phi coefficient using the phi
function directly from the contingency table:
# Example contingency table:
# 40 17
# 11 45
library(effectsize)
<- matrix(c(40, 11,
contingency_table 17, 45),ncol = 2)
<- phi(contingency_table, alternative = "two.sided")
phi_coefficient
phi_coefficient
Phi (adj.) | 95% CI
-------------------------
0.50 | [0.31, 0.69]
In our example we obtained a phi coefficient of \(\phi\) = .50 [0.31, 0.69].
Cramer’s V, sometimes also referred to as Cramer’s phi (\(\phi\)), is a generalized effect size measure of the association between two nominal variables. It applies to contingency tables of any size (\(2\times 2\), \(3\times 3\), \(3\times 4\), \(5\times 3\), etc.). Cramer’s \(V\) on a \(2\times 2\) contingency table is equivalent to the phi coefficient. For an illustration of a higher order contingency table, Table 9.2 represents a \(3\times 4\) contingency table of two variables.
\(X=0\) | \(X=1\) | \(X=2\) | \(X=3\) | |
---|---|---|---|---|
\(Y=0\) | \(n_{00}\) | \(n_{10}\) | \(n_{21}\) | \(n_{31}\) |
\(Y=1\) | \(n_{01}\) | \(n_{11}\) | \(n_{21}\) | \(n_{31}\) |
\(Y=2\) | \(n_{02}\) | \(n_{12}\) | \(n_{22}\) | \(n_{32}\) |
Similarly to the phi coefficient, the value of Cramer’s \(V\) ranges from 0 to 1 and can interpreted in a similar way to a phi coefficient. Again we can use the \(\chi^2\) statistic to compute the value, however, since there can be more than 2 levels to each variable, we also need to take into account the number of levels, \(k\), of the variable with the least number of levels (e.g., a \(3 \times 4\) contingency table, \(k\) would be equal to 3). Cramer’s \(V\) is defined as (equation 7.2.6, Cohen 1988),
\[ V = \sqrt{\frac{\chi^2}{n(k-1)}} \tag{9.5}\]
The standard error of a Cramer’s \(V\) is similar to that of a Pearson correlation and a \(\phi\) coefficient.
\[ SE_V = \sqrt{\frac{\left(1-V^2\right)^2}{n-1}} \tag{9.6}\]
Where \(n\) is the total sample size (i.e., the sum of all cells). Like the pearson correlation, we can not calculate the confidence interval directly from the standard error, instead, we must convert \(V\) to a Fisher’s Z statistic, \(Z_V = \text{arctanh}(V)\). We can then calculate the 95% confidence interval for \(V\) by back-transforming the confidence interval for \(Z_V\):
\[ SE_{Z_V} = \frac{1}{\sqrt{n-3}} \tag{9.7}\]
\[ CI_{V} = \tanh(Z_V \pm 1.96\times SE_{Z_V}) \tag{9.8}\]
Using the ufs
package (Peters and Gruijters 2023), we can calculate Cramer’s \(V\) and it’s 95% confidence interval using the Fisher’s Z method described above. For the example, we can example data from a 3 \(\times\) 3 contingency table.
# Example contingency table:
# 40 14 12
# 11 27 9
# 5 10 34
library(ufs)
<- matrix(c(40, 11, 5,
contingency_table 14, 27, 10,
12, 9, 34),ncol = 3)
<- cramersV(contingency_table)
V <- confIntV(contingency_table)
CI
# print pearson correlation and confidence intervals
data.frame(V = MOTE::apa(V$output$cramersV),
Vlow = MOTE::apa(CI$output$confIntV.fisher[1]),
Vhigh = MOTE::apa(CI$output$confIntV.fisher[2]))
V Vlow Vhigh
1 0.442 0.309 0.558
In our example we obtained a Cramer’s \(V\) of \(V\) = .44 [.31, .56].
Cohen’s \(h\) is a measure of distance between two proportions or probabilities. It is sometimes also referred to as the “difference between arcsines”. For a given proportion \(p\), its arcsine transformation is given by (equation 6.2.1, Cohen 1988):
\[ \psi = 2\cdot \text{arcsin}(\sqrt{p}). \tag{9.9}\]
Cohen’s \(h\) is the difference between the arcsine transformations of two proportions (equation 6.2.2, Cohen 1988):
\[ h = \psi_1 - \psi_2 \tag{9.10}\]
Cohen’s \(h\) is commonly used for the power analysis of proportion tests. In fact, it is the required effect size measure in the program G Power (Faul et al. 2009). We can calculate the standard error of Cohen’s \(h\),
\[ SE_h = \sqrt{\frac{1}{n_1} + \frac{1}{n_2}} \tag{9.11}\]
Since the sampling distribution of \(h\) is symmetric, we can calculate the confidence intervals from the standard error,
\[ CI_h = h \pm1.96\times SE_h \tag{9.12}\]
To calculate Cohen’s \(h\), we can use the cohens_h
function in the effectsize
package in R.
# install package if not done so already
# install.packages('effectsize')
# Example proportions: p1 = .45, p2 = .30
library(effectsize)
<- matrix(c(40, 11,
contingency_table 14, 27),ncol = 2)
cohens_h(contingency_table)
Cohen's h | 95% CI
------------------------
0.93 | [0.52, 1.34]
From the example, the R code outputted a Cohen’s \(h\) value of \(h\) = .93 [0.52, 1.34].
Cohen’s \(w\) is a measurem of association analogous to the phi coefficient but on tables that are larger than 2x2. Although Cohen’s \(w\) is useful for power analyses, it is not so useful as a stand-alone effect size. As Cohen (1988) states (pp. 221):
As a measure of association, [Cohen’s \(w\)] lacks familiarity and convenience
Cohen’s \(w\) has the exact same formula as the phi coefficient with the only difference being that the \(\chi^2\) statistic comes from a contingency table of any size (equation 7.2.5, Cohen 1988),
\[ w = \sqrt{\frac{\chi^2}{n}} \tag{9.13}\]
And can also be calculated directly from Cramer’s \(V\) (equation 7.2.7, Cohen 1988),
\[ w = V \times \sqrt{k-1} \tag{9.14}\]
Where \(k\) is the number of categories in the variable with the least number of categories. We can use the cohens_w()
function in the effectsize
package (Ben-Shachar, Lüdecke, and Makowski 2020).
# Example contingency table
# 40 14
# 11 27
<- matrix(c(40, 11,
contingency_table 14, 27),ncol = 2)
cohens_w(contingency_table,
alternative = "two.sided")
Cohen's w | 95% CI
------------------------
0.45 | [0.24, 0.65]
From the example code, the cohens_w
function returned Cohen’s \(w\) value of \(w\) = .45 [0.24, 0.65].
Ben-Shachar et al. (2023) introduced a new effect size for one-dimensional tables of counts/proportions that they label with the Hebrew letter, פ. Ben-Shachar’s פ is a correction to Cohen’s \(w\) that adjusts for the expected value and consequently bounds the value between 0 and 1. The equation for פ is defined as,
\[ \mathbf{פ }= \sqrt{\frac{\chi^2}{n \left(\frac{1}{\min\left(P_E\right)} -1\right)}} \tag{9.15}\]
Where \(\min(P_E)\) is the smallest expected probability. The formula for Ben-Schachar’s פ can be also be expressed in terms of Cohen’s \(\omega\),
\[ \mathbf{פ }= \frac{\omega}{\sqrt{\left(\frac{1}{\max(P_E)} -1\right)}} \tag{9.16}\]
In R, we can calculate Ben-Shachar’s פ using the fei()
function in the effectsize
package (Ben-Shachar, Lüdecke, and Makowski 2020).
# Example:
# Observed counts: 20, 50, 100 (observed proportions: .12, .29, .59)
# Expected proportions: .5, .2, .3
<- c(20,50,100)
observed_counts <- c(.5,.2,.3)
expected_probabilities
fei(observed_counts,
p = expected_probabilities,
alternative = "two.sided")
Fei | 95% CI
-------------------
0.39 | [0.31, 0.47]
- Adjusted for uniform expected probabilities.
From the example code, the fei
function returned Ben-Shachar’s פ value of .39 [0.31, 0.47].
Odds ratio measures the effect size between two binary variables. It is commonly used in medical and behavioral intervention research, and notably, in meta-analysis.
Let’s imagine a study conducted to investigate the association between smoking and the development of major depressive disorder (MDD). The study includes a sample of 251 individuals, categorizing them into two groups: 125 smokers and 126 non-smokers. The researchers are interested in understanding the odds of having major depressive disorder (MDD) among smokers compared to non-smokers. Say we find that 25 smokers were diagnosed with MDD while 100 were not, but in the non-smoker group, 12 individuals were diagnosed with MDD while 120 were not. The odds ratio would then be:
\[ OR = \frac{25/100}{12/120}= \frac{.25}{.10} = 2.50 \tag{9.17}\]
In general, we can can compute the odds-ratio from a contingency table between binary variables \(X\) (i.e., the treatment) and \(Y\) (i.e., the outcome; see Table 9.3).
\(X=T\) | \(X=C\) | |
---|---|---|
\(Y=0\) | \(n_{T0}\) | \(n_{C0}\) |
\(Y=1\) | \(n_{T1}\) | \(n_{C1}\) |
Ultimately, we want to compare the outcome between the treatment group (\(X=T\)) and the control group (\(X=C\)). Therefore we can compute the odds ratio as,
\[ OR = \frac{n_{T1}/n_{T0}}{n_{C1}/n_{C0}} \tag{9.18}\]
The standard distribution of the odds-ratio is asymmetric. To calculate confidence intervals, we can first convert the odds ratio to a log odds ratio (\(LOR= \log(OR)\)). Then we can calculate the standard error of the log odds ratio,
\[ SE_{LOR} = \sqrt{\frac{1}{n_{T0}} + \frac{1}{n_{T1}} + \frac{1}{n_{C0}} + \frac{1}{n_{C1}}} \tag{9.19}\]
With the standard error of the log odds ratio we can then calculate the confidence interval of the odds ratio by back-transforming using the exponential function,
\[ CI_{OR} = \exp(LOR \pm 1.96\times SE_{LOR}) \tag{9.20}\]
In R, we can use the effectsize
package to calculate the odds ratio and it’s confidence interval:
# Example:
# Treatment Group: 10 diseased, 43 healthy
# Control Group: 24 diseased, 41 healthy
<- matrix(c(10, 24,
contingency_table 43, 41),ncol = 2)
oddsratio(contingency_table,
alternative = "two.sided")
Odds ratio | 95% CI
-------------------------
0.40 | [0.17, 0.93]
The code output for this example shows an odds ratio of \(OR\) = 0.40 [0.17, 0.93]
Risk difference can be used to interpret the difference between two proportions. If we use the contingency table from Table 9.3, and calculate a risk difference between the treatment group and the control group. We can first calculate the proportion of cases where the outcome is \(Y=1\) within the control group and the treatment group:
\[ p_C=\frac{n_{C1}}{n_{C0}+n_{C1}} \tag{9.21}\]
\[ p_T=\frac{n_{T1}}{n_{T0}+n_{T1}} \tag{9.22}\]
Then using these proportions we can calculate the risk difference (\(RD\)),
\[ RD = p_T - p_C. \tag{9.23}\]
The corresponding standard error is,
\[ SE_{RD} = \sqrt{\frac{p_C(1-p_C)}{n_C} + \frac{p_T(1-p_T)}{n_T} } \tag{9.24}\]
Where \(n_C\) and \(n_T\) are the total sample sizes within the control and treatment group, respectively. The standard error can then be used to compute the 95% confidence intervals,
\[ CI_{RD} = RD \pm 1.96 \times SE_{RD} \tag{9.25}\]
The risk difference formula is fairly simple, so we can compute it using base R.
# Example:
# Treatment group: proportion of cases = .5, sample size = 40
# Control group: proportion of cases = .3, sample size = 45
<- .50
pT <- .30
pC <- 40
nT <- 45
nC
<- pT - pC
RD
<- sqrt( pC*(1-pC)/nC + pT*(1-pT)/nT )
SE
# compute 95% CIs
<- RD - 1.96*SE
RDlow <- RD + 1.96*SE
RDhigh
data.frame(
RD = MOTE::apa(RD),
RDlow = MOTE::apa(RDlow),
RDhigh = MOTE::apa(RDhigh)
)
RD RDlow RDhigh
1 0.200 -0.005 0.405
The relative risk, often referred to as the “risk ratio,” calculates the ratio between the proportion of cases in the treatment group and the proportion of cases in the control group. It provides a straightforward interpretation: “individuals receiving the treatment have a \(RR\) times higher odds of experiencing the outcome compared to controls.” To calculate relative riskm, first we need to calculate the proportion of outcome cases in the treatment and control group
\[ p_C=\frac{n_{C1}}{n_{C0}+n_{C1}} \]
\[ p_T=\frac{n_{T1}}{n_{T0}+n_{T1}} \]
Then we can calculate the relative risk,
\[ RR=\frac{p_T}{p_C} \tag{9.26}\]
The corresponding standard error can be computed as,
\[ SE_{RR} = \sqrt{\frac{p_T}{n_T} + \frac{p_C}{n_C}} \tag{9.27}\]
The confidence intervals can be computed from the standard error,
\[ CI_{RR} = RR\pm 1.96\times SE_{RR} \tag{9.28}\]
To compute relative risk, we can simply use the equations above in base R.
# Example:
# Treatment Group: 10 diseased, 43 healthy, 53 total
# Control Group: 24 diseased, 41 healthy, 65 total
<- 10/(43+10)
pT <- 24/(41+24)
pC <- 53
nT <- 65
nC
<- pT / pC
RR
<- sqrt(pT/nT + pC/nC)
SE
<- RR - 1.96*SE
RRlow <- RR + 1.96*SE
RRhigh
# print pearson correlation and confidence intervals
data.frame(RR = MOTE::apa(RR),
RRlow = MOTE::apa(RRlow),
RRhigh = MOTE::apa(RRhigh))
RR RRlow RRhigh
1 0.511 0.323 0.699