12.5 – Effect size for ANOVA

Introduction
Effect size estimation
R code
Questions
Quiz
Data set
Chapter 12 contents

Introduction

As noted in the t-test chapter and our discussion of statistical power, an effect size is a measure of the strength of a phenomenon. An effect size calculated from data is a descriptive statistic that indicates how large (or small) the difference is between two or more samples. Effect size measures help to provide context between statistical significance, ie, p-values. Among other uses, effect size measures play an important role in meta-analysis studies biological significance (ie, importance of the measured difference).

Estimation of effect size

eta-squared (η²)

$\begin{align*} \eta^2 = \frac{SS_{B}}{SS_{T}} \end{align*}$

where SS_B and SS_T are from the one-way ANOVA and refer to the sum of squared among (between) groups and the total sum of squares, respectively. $\eta ^2$ measures the proportion of the variation in the response variable that can be explained by membership in one of the groups.

For example, $\eta ^2$ of 0.15 or 15% is interpreted to mean that just 15% of the variation in the response variable can be attributed to membership in the groups (ie, whether a subject was in control group or treatment group, only 15% of the differences between individuals can be attributed to having received control).

Note 1: η² is just the unadjusted coefficient of determination, R².

omega-squared (ω²)

A similar, but “less biased” effect size estimate is given by omega-squared (ω²) (Okada 2013). (The bias comes in because sample estimates were used.) Interpretation of $\omega ^2$ is the same as $\eta ^2$ : measures the proportion of variation explained by membership in the groups. $\omega ^2$ is given as

$\begin{align*} \omega^2 = \frac{SS_{B}-\left (k-1 \right ) MS_{W}}{SS_{T} + MS_{W}} \end{align*}$

The bias for $\eta ^2$ is more pronounced with small sample size, so omega-squared is preferred. $\omega ^2$ will always be less than or equal to $\eta ^2$ (Okada 2013).

R code

If you have not already done so, please install the package effectsize.

We’ll use the example one-way ANOVA problem from Chapter 12.2. I’ve added the data as example.ch12 and some R script to load the data, scroll to bottom of this page or click on the link).

Get $\eta^2$

First, get the ANOVA model. There are several ways to do this within R (and Rcmdr), the most general is to use the general linear model function, lm(),

#Get the model
AnovaModel.4 <- lm(Difference~Group, data=example.ch12)
#turn it into format that the next command needs
aov_fit <- anova(AnovaModel.4)
#get eta
effectsize(aov_fit)

R output

Parameter | Eta2 (partial) | 9e+01% CI
-----------------------------------------
Group | 0.95 | [0.88, 0.97]

$\eta^2$ is 0.95 for the Difference group.

Get $\omega ^2$

R code

omega_squared(aov_fit, partial = FALSE, ci = 0.95)

R output

Parameter | Omega2 | 1e+02% CI
---------------------------------
Group | 0.93 | [0.82, 0.97]

This is a case of a large effect due to group membership. The differences between A, B, and C means account for 93% of the variation. Note that both $\eta^2$ and $\omega^2$ are close; note also that $\eta^2$ is greater than $\omega^2$ .

Questions

In Chapter 7.4 we introduced the concept of number needed to treat, NNT. Discuss the concept of “important differences” between sample means of a control group and a treatment group to NNT.
This question refers to hypothetical data presented in Chapter 12.2. For the “No difference” group, calculate $\eta^2$ and $\omega^2$ . Use the boxplots shown in Figure 1, together with effect size statistics to discuss the relationship between statistical significance (p-value) and important difference.

Quiz Chapter 12.5

Effect size for ANOVA

Data set

example.ch12 <- read.table(header=TRUE, sep=",",text="
Group, No.difference, Difference 
A, 12.04822, 11.336161 
A, 12.67584, 13.476142 
A, 12.99568, 12.961210 
A, 12.01745, 11.746712 
A, 12.34854, 11.275492 
B, 12.17643, 7.167262 
B, 12.77201, 5.136788 
B, 12.07137, 6.820242 
B, 12.94258, 5.318743 
B, 12.07670, 7.153992 
C, 12.58212, 3.344218 
C, 12.69263, 3.792337 
C, 12.60226, 2.444438 
C, 12.02534, 2.576014 
C, 12.60420, 4.575672")
#check the dataframe
head(example.ch12)

Chapter 12 contents

Introduction
The need for ANOVA
One way ANOVA
Fixed effects, random effects, and agreement
ANOVA from “sufficient statistics”
Effect size for ANOVA
ANOVA posthoc tests
Many tests one model
References and suggested readings

7.6 – Confidence intervals

Quantifying uncertainty
Confidence interval for proportions
Confidence interval for ARR
Confidence interval for NNT
Odds ratio Standard error and 95% confidence interval
P-value from confidence intervals
Questions
Quiz
Chapter 7 contents

Quantifying uncertainty.

Although I’ve already presented the concept (eg, Chapter 3.4), and equations for confidence intervals of risk estimates (Chapter 7.4, Chapter 7.5), here, we expand on the idea of confidence intervals. Confidence intervals are a central part of meeting one of the main objectives of statistics, that is, estimation and communicating uncertainty with any reported estimate.

There are three components of statistical analysis

Estimation
Inference
Modeling

Inference refers to statistical hypothesis testing — we ask questions of observations — do men (Rice et al 1999) and women (Fisher et al 2012) differ for blood glucose levels following a bout of aerobic exercise? T-tests, analysis of variance (ANOVA), chi-square, correlation, regression are types of statistical procedures used to do statistical inference. Modeling on the other hand refers to procedures used to relate cause and effect or for use in prediction of future values given new data. Many of the statistical procedures one uses for inference are also used to build statistical models (ANOVA, regression). Studies may intend to either test some hypothesis (inference) or to provide a predictive equation (modeling). But most studies that relate observations gathered from an experiment are obliged to also report statistics, and this is the realm of estimation. Estimates of the mean and standard deviation, for example, would be typical statistics one expects to find in a report. We call these descriptive statistics, and together with graphics, descriptive statistics are the chief way we describe our results.

Confidence intervals are crucial in statistics because they provide a range of values where the true population parameter is likely to fall, offering a more nuanced understanding of the estimate beyond just a single point value. Similarly, for statistical model model building, confidence intervals quantify the uncertainty inherent in data, providing a range of plausible values for model parameters rather than a single point estimate.

We will review how to calculate confidence intervals for proportions and for NNT. These intervals are available in epiR package and automatically returned in RcmdrPlugin.EBM.

Confidence interval for proportions

A proportion is the fraction of individuals in a population with some characteristic. The characteristic might be HIV positive for example. This would be called the population proportion and it would be a parameter of interest. In reality, we calculate a sample proportion and therefore estimate the population proportion with error. We can calculate the confidence interval (CI) of the proportion to communicate the precision of our estimate. For proportions, we use the binomial distribution — either a sample has the characteristic of interest or it does not — there are only two possibilities. There are a variety of ways to go here, and the simplest is to use a normal approximation. This will work well provided the sample size was reasonably large and the proportion is not close to zero or one, that is, we invoke the Central Limit Theorem here. Although the outcomes are binomial, the error is assumed to be normally distributed. The Wald confidence interval for p is

$\begin{align*} \% 95 \ CI = \hat{p} \pm z\sqrt{\frac{1}{n}\cdot \hat{p}\left ( 1- \hat{p} \right )} \end{align*}$

where $\hat{p}$ is the proportion of individuals with the characteristic (also called successes), z is the percentile from the normal distribution that corresponds to 1 – 1/2α. For 95% CI then α = 0.05 and z would be 1.965. (See standard normal table.) Of course, if making the normal approximation for the binomial is not appropriate, the CI is less than ideal. The binomial after all is a discrete distribution whereas the normal distribution is continuous, so errors will enter particularly for low sample numbers.

Other approaches may be used to get better estimates of CI for proportions, including Wilson score intervals and Jeffrey Intervals (Agresti and Coull 1998). See R package propCIs.

Because a statistic like the mean or a calculation of absolute or relative risk reduction are calculated from samples drawn from a population, the estimate comes with error. The error is basically this – if we calculate a statistic like number needed to treat (NNT) or its converse, the number needed to harm (NNH), we need to communicate to the reader how precise our estimate is. Estimation has to do with accuracy, error, and precision.

Confidence interval for ARR

The ARR is simply $\frac{c}{n_{2}}-\frac{a}{n_{1}}$

where $n_{1}$ is the number of treated or exposed individuals for which the event occurred and $n_{2}$ is the number of untreated or unexposed individuals which the event occurred.

	Event happened	Event did not happen
Treated or Exposed	a	b
Control or Not exposed	c	d

Our data from example (Ch 7.5), were a = 612, b = 2192, c = 375, and d = 2543.

$\begin{align*} SE_{ARR}=\sqrt{\left (\frac{a}{n_{1}} \cdot \frac{1-\frac{a}{n_{1}}}{n_{1}} \right )+\left (\frac{c}{n_{2}} \cdot \frac{1-\frac{c}{n_{2}}}{n_{2}} \right )} \end{align*}$

and the 95% confidence interval is then $ARR \pm 2 \cdot SE_{ARR}$

approximately. The “2” is only approximate; you need to use Z = 1.965, the value at probability value = 0.9725 (which comes from the Normal Table).

Confidence interval for NNT

For a sample of 100 people drawn at random from a population (which may number in the millions), then repeat the NNT calculation for a different sample of 100 people, do we expect the first and second NNT estimates to be EXACTLY the same number? No, but we do expect them to be close and we can define what we mean by close as we expect each estimate to be within certain limits. While we expect the second calculation to be close to the first estimate, we would be surprised if it was EXACTLY the same. And so, which is the correct estimate, the first or the second? They both are, in the sense that they both estimate the parameter NNT (a property of a population). But we can do better than two estimates. Confidence Intervals (CI) allow us to assign a probability to how certain we are about the statistic and whether it is likely to be close to the true value. We will

For CI of NNT, we need sample size for control and treatment groups; like all confidence intervals, we need to calculate the standard error of the statistic, in this, case, the standard error (SE) for NNT.

SE = sqrt(risk placebo * (1 - risk placebo) / (# in placebo group) + risk treatment * (1 - risk treatment) / (# in treatment group))

where SE is the standard error for NNT

The CI is approximately then $NNT \pm 2 X SE$

* – the “2” is only approximate; you need to use Z = 1.965, the value at probability value = 0.9725 (by multiplying this, which comes from the Standard Normal Table).

Odds ratio Standard error and 95% confidence interval

Like any statistic we can calculate, an estimate of odds ratio should be accompanied by the confidence limit. The standard error may be calculated with the following formula

$\begin{align*} SE\left \{ ln\left ( OR \right ) \right \}=\sqrt{\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d}} \end{align*}$

R code

seOdds <- sqrt(sum(1/612, 1/2192,1/375,1/2543))

where ln refers to the natural logarithm. An estimate for the 95% confidence interval is

$\begin{align*} lower \ limit = exp\left (ln\left ( OR \right ) \right )-1.96\cdot SE\left \{ ln\left ( OR \right ) \right \} \end{align*}$

where exp is the exponential function

$e^{x}$

and ln is natural logarithm.

R code

 exp(log(1.89,base=exp(1)) - 1.96*seOdds)

Our example the lower limit was 1.64. For the upper limit

$\begin{align*} upper \ limit = exp\left (ln\left ( OR \right ) \right )+1.96\cdot SE\left \{ ln\left ( OR \right ) \right \} \end{align*}$

Our example the upper limit was 2.19.

R code

 exp(log(1.89,base=exp(1)) + 1.96*seOdds)

Thus, our estimate was 1.89 and the 95% confidence interval was $\left ( 1.64,2.19\right )$

which does not include one. Therefore, we conclude statistically different.

Note 1. A reminder, log() without specifying the base returns the natural logarithm. A better coding habit might suggest always specifying the base:

log(x, base = exp(1))

P-values from confidence intervals

While we expect certain reporting criteria for published, it is not uncommon that one or more elements are missing. For example, while we show how to obtain confidence intervals, the p-values were not reported. We discuss statistical inference and p-values in detail in Chapter 8, but for now, a low p-value (typically < 0.05) indicates that the odds ratio (OR) is statistically significant, meaning the observed relationship between a factor and an outcome is unlikely to be due to random chance. An OR of the difference may not be reported, or, an inexact p-value is reported like “< 0.05,” but we “need” an exact p-value for our meta-analysis. A little math, and these missing statistics can be calculated.

Note 2. We illustrated how to calculate confidence intervals from p-values in Chapter 3.4 – Estimating parameters. See two helpful “how to” articles by Bland and Altman in 2011 BMJ (references listed in Chapter 3.6 – References).

We can use the OR example. We need to take the logarithm of ratios like OR or NNT or RRR. From above, we have the confidence intervals.

We need

The OR estimate, $myEST$ , which was 1.89

lower and upper confidence intervals: 1.64 and 2.19, respectively.

the standard error, which can be calculated as $mySE = (upper − lower)/(2×1.96)=0.07377$

the test statistic, $myZ <- abs(myEST/mySE)$

and then the p-value, $exp(-.717*myZ-0.416*myZ^2) = 0.4044$

we would then report that the OR was not statistically different, p = 0.4044).

Questions

1. Instead of 95% confidence interval, obtain the 99% confidence interval for odds ratio 1.89.

Quiz Chapter 7.6

Confidence intervals

Chapter 7 contents

7.4 – Epidemiology: Relative risk and absolute risk, explained

Risk communication
An example of ARR and RRR risk calculations using natural numbers
Confidence interval for NNT is derived from the Confidence interval for ARR
Software: epiR
Improve output from epiR
Software: epitools
Software: OpenEpi
Software: RcmdrPlugin.EBM
Questions
Quiz
Chapter 7 contents

Risk communication.

Epidemiology is the study of patterns of health and illness of populations. An important task in an epidemiology study is to identify risks associated with disease. Epidemiology is a crucial discipline used to inform about possible effective treatment approaches, health policy, and about the etiology of disease.

Please review terms presented in section 7.1 before proceeding. RR and AR are appropriate for cohort-control and cross-sectional studies (see 2.4 and 5.4) where base rates of exposure and unexposed or numbers of affected and non-affected individuals (prevalence) are available. Calculations of relative risk (RR) and relative risk reduction (RRR) are specific to the sampled groups under study whereas absolute risk (AR) and absolute risk reduction (ARR) pertain to the reference population. Relative risks are specific to the study, absolute risks are generalized to the population. Number needed to treat (NNT) is a way to communicate absolute risk reductions.

An example of ARR and RRR risk calculations using natural numbers.

Clinical trials are perhaps the essential research approach (Sibbald and Roland 1998; Sylvester et al 2017); they are often characterized with a binary outcome. Subjects either get better or they do not. There are many ways to represent risk of a particular outcome, but where possible, using natural numbers is generally preferred as a means of communication, particularly to the general public. Consider the following example (pp 34-35, Gigerenzer 2015): What is the benefit of taking a cholesterol-lowering drug, Pravastatin, on the risk of deaths by heart attacks and other causes of mortality? Press releases (eg, Maugh 1995), from the study stated the following:

“… the drug pravastatin reduced … deaths from all causes 22%”

A subsequent report (Skolbekken 1998) presented the following numbers (Table 1).

Table 1. Reduction in total mortality (5 year study) for people who took Pravastatin compared to those who took placebo.

Deaths per 1000 people with high cholesterol
(> 240 mg/dL) No deaths Cumulative
incidence

Treatment

	Deaths per 1000 people with high cholesterol (> 240 mg/dL)	No deaths	Cumulative incidence
Treatment	Pravastatin (n = 3302)	`a` `= 32`	`b` `= 3270`	CIe
Placebo (n = 3293)	`c` `= 41`	`d` `= 3252`	CIu

Pravastatin

(n = 3302)

a
= 32

b
= 3270

CIe

Placebo

(n = 3293)

c
= 41

d
= 3252

CIu

where cumulative incidence refers to the number of new events or cases of disease divided by the total number of individuals in the population at risk.

Do the calculations of risk

The risk reduction (RR), or the number of people who die without treatment (placebo) minus those who die with treatment (Pravastatin), $41-32=9$

$\begin{align*} RR = \frac{\frac{a}{a+b}}{\frac{c}{c+d}} = 0.91 \end{align*}$

and the cumulative incidence in the exposed (treated) group, CIe, is $32/(32+3270) = 0.0097$ , and cumulative incidence in the unexposed (control) group, CIu, is $41/(41+3252) = 0.01245$ . We can calculate another statistic called the risk ratio,

$\begin{align*} RR = \frac{CI_{e}}{CI_{u}} = 0.78 \end{align*}$

because the risk ratio is less than one, we interpret that statins reduce the risk of mortality from heart attack. In other words, statins lowered the risk by 0.78.

But is this risk reduction meaningful?

Now, consider the absolute risk reduction (ARR) is $0.9 \% = 100 \% X \frac{9}{1000}$

Relative risk reduction, or the absolute risk reduction divided by the proportion of patients who die without treatment, is $22 \% = 100 \% X 9 \div 41$

Conclusion: high cholesterol may contribute to increased risk of mortality, but the rate is very low in the population as a whole (the ARR).

Another useful way to communicate benefit is to calculate the Number Needed to Treat (NNT), or the number of people who must receive the treatment to save (benefit) one person. The ideal NNT is a value of one (1), which would be interpreted as everyone improves who receives the treatment. By definition, NNT must be positive; however, a resulting negative NNT would suggest the treatment may cause harm, ie, number needed to harm (NNH).

For this example, the NNT is

$\begin{align*} \frac{1}{\frac{9}{1000}} = 111 \end{align*}$

therefore, to benefit one person, 111 need to be treated. The flip side of the implications of NNT, although one person may benefit by taking the treatment, 110 ( $=111 - 1$ ) will take the treatment, but will NOT RECEIVE THE BENEFIT, but do potentially get any side effect of the treatment.

Confidence interval for NNT is derived from the Confidence interval for ARR.

For a sample of 100 people drawn at random from a population (which may number in the millions), then repeat the NNT calculation for a different sample of 100 people, do we expect the first and second NNT estimates to be exactly the same number? No, but we do expect them to be close and we can define what we mean by close as we expect each estimate to be within certain limits. While we expect the second calculation to be close to the first estimate, we would be surprised if it was exactly the same. And so, which is the correct estimate, the first or the second? They both are, in the sense that they both estimate the parameter NNT (a property of a population).

We use confidence interval to communicate where we believe the true estimate for NNT to be. Confidence Intervals (CI) allow us to assign a probability to how certain we are about the statistic and whether it is likely to be close to the true value (Altman 1998, Bender 2001). We will calculate the 95% CI for the ARR using the Wald method, then take the inverse of these estimates for our 95% CI. The Wald method assumes normality.

For CI of ARR, we need sample size for control and treatment groups; like all confidence intervals, we need to calculate the standard error of the statistic, in this, case, the standard error (SE) for ARR is approximately

$\begin{align*} SE_{\left ( p_{1}-p_{2} \right )}=\sqrt{\frac{p_{1}\left ( 1-p_{1} \right )}{n_{1}}+\frac{p_{2}\left ( 1-p_{2} \right )}{n_{2}}} \end{align*}$

where SE is the standard error for ARR. For our example, we have

$\begin{align*} SE_{\left ( p_{1}-p_{2} \right )}=\sqrt{\frac{0.041\left ( 1-0.041 \right )}{1000}+\frac{0.032\left ( 1-0.032 \right )}{1000}} = 0.00838 \end{align*}$

The 95% CI for ARR is approximately $ARR\pm 2 X SE_{\left ( p_{1}-p_{2} \right )}$

For the Wald estimate, replace the 2 with $z = 1.965$ , which comes from the normal table for z at $\frac{0.95}{2}$

(why the 2 in the equation? Because it is plus or minus so we divide the frequency 0.95 in half) and for our example, we have $0.009 \pm 2 X SE_{\left ( p_{1}-p_{2} \right )} = \left (-0.0078,0.0258 \right )$

and the inverse for NNT CI is $\left (-128, 38 \right )$

Our example exemplifies the limitation of the Wald approach (cf. Altman 1998): our confidence interval includes zero, and doesn’t even include our best estimate of NNT (111).

Note 1: By now I trust you can see differences for results by direct input of the numbers into R and what you get by the natural numbers approach. In part this is because we round in our natural number calculations — remember, while it makes more sense to communicate about whole numbers (people) and not fractions (fractions of people!), rounding through the calculations adds error to the final value. As long as you know the difference and the relevance between approximate and exact solutions, this shouldn’t cause concern.

Software: epiR.

R has many epidemiology packages, epiR and epitools are two. Most of the code presented stems from epiR.

We need to know about our study design in order to tell the functions which statistics are appropriate to estimate. For 2X2 tables, we will use epi.2by2() command, which allows four kinds of studies (Table 2).

Table 2. The method argument in epi.2by2 chooses how to analyze your data based on your study design.

Method	Description	Notes
cohort.count	Use this when you follow a set number of people over time and count how many get sick.	This design directly calculates Relative Risk Reduction (RRR) and Absolute Risk Reduction (ARR) because you track baseline risk over time.
cohort.time	Use this when you follow people over time and track exactly how long each person was monitored before getting sick.	This calculates rate reductions instead of ARR or RRR because it deals with person-time data rather than raw counts.
case.control	Use this when you start with sick people and healthy people and look backwards at their past exposures.	You cannot calculate ARR or RRR here because you do not know the actual baseline risk of the population.
cross.sectional	Use this when you take a snapshot of a population at a single point in time to see who is currently exposed and sick.	This calculates prevalence reductions rather than true ARR or RRR because you are not tracking changes over time.

For our statin example, the design was prospective cohort (ie, cohort.count in epiR package language), not case-control or cross-sectional (review in Chapter 5.4).

library(epiR)
Table1 <- matrix(c(32,3270,41,3252), 2, 2, byrow=TRUE, dimnames = list(c("Statin", "Placebo"), c("Died", "Lived")))
Table1
        Died Lived
Statin    32  3270
Placebo   41  3252
epi.2by2(Table1, method="cohort.count", outcome = "as.columns")

R output

           Outcome + Outcome - Total           Inc risk *
Exposed +         32      3270  3302 0.97  (0.66 to 1.37)
Exposed -         41      3252  3293 1.25  (0.89 to 1.69)
Total             73      6522  6595 1.11  (0.87 to 1.39)

Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio                          0.78 (0.49, 1.23)
Inc odds ratio                          0.78 (0.49, 1.24)
Attrib risk in the exposed *           -0.28 (-0.78, 0.23)
Attrib fraction in the exposed (%)    -28.48 (-103.47, 18.88)
Attrib risk in the population *        -0.14 (-0.59, 0.32)
Attrib fraction in the population (%) -12.48 (-37.60, 8.05)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 1.147 Pr>chi2 = 0.284
Fisher exact test that OR = 1: Pr>chi2 = 0.292
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units

The risk ratio we calculated by hand is shown in green in the R output, along with other useful statistics (see ?epi2x2 for help with these additional terms) not defined in our presentation.

We explain results of chi-square goodness of fit (Ch 9.1) and Fisher exact (Ch 9.5) tests in Chapter 9. Suffice to say here, we interpret the p-value (Pr) = 0.284 and 0.292 to indicate that there is no association between mortality from heart attacks with or without the statin (ie, the Odds Ratio, OR, not statistically different from one).

Wait! Where’s NNT and other results?

Use another command in epiR package, epi.tests(), to determine the specificity, sensitivity, and positive (or negative) predictive value.

epi.tests(Table1)

and R returns

          Outcome + Outcome - Total
Test +           32      3270  3302
Test -           41      3252  3293
Total            73      6522  6595

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                    0.50 (0.49, 0.51)
True prevalence *                        0.01 (0.01, 0.01)
Sensitivity *                            0.44 (0.32, 0.56)
Specificity *                            0.50 (0.49, 0.51)
Positive predictive value *              0.01 (0.01, 0.01)
Negative predictive value *              0.99 (0.98, 0.99)
Positive likelihood ratio                0.87 (0.67, 1.13)
Negative likelihood ratio                1.13 (0.92, 1.38)
False T+ proportion for true D- *        0.50 (0.49, 0.51)
False T- proportion for true D+ *        0.56 (0.44, 0.68)
False T+ proportion for T+ *             0.99 (0.99, 0.99)
False T- proportion for T- *             0.01 (0.01, 0.02)
Correctly classified proportion *        0.50 (0.49, 0.51)
--------------------------------------------------------------
* Exact CIs

Additional statistics are available by saving the output from epi2x2() or epitests() to an object, then using summary(). For example save output from epi.2by2(Table1, method="cohort.count", outcome = "as.columns") to object myEpi, then

summary(myEpi)

look for NNT in the R output

$massoc.detail$NNT.strata.wald
       est    lower    upper
1 -362.377 -128.038  436.481

Thus, the NNT was 362 (compared to 111 we got by hand) with 95% Confidence interval between -436 and +128 (make it positive because it is a treatment improvement).

Note 2: Strata (L. layers), refer to subgroups, for example, sex or age categories (see discussion in Ch05.4). Our examples here are not presented as subgroup analysis, but epiR reports by name strata.

epiR reports a lot of additional statistics in the output and for clarity, I have not defined each one, just the basic terms we need for BI311. As always, see help pages (eg, ?epi.2x2 or ?epitests)for more information about structure of an R command and the output.

We’re good, but we can work the output to make it more useful to us.

Improve output from epiR.

For starters, if we set interpret=TRUE instead of the default, interpret=FALSE, epiR will return a richer response.

fit <- epi.2by2(dat = as.table(Table1), method = "cohort.count", conf.level = 0.95, units = 100, interpret = TRUE, outcome = "as.columns")
fit

R output. In addition to the table of coefficients (above), interpret=TRUE provides more context, shown below

Measures of association strength:
The outcome incidence risk among the exposed was 0.78 (95% CI 0.49 to 1.23) times less than the outcome incidence risk among the unexposed. 

The outcome incidence odds among the exposed was 0.78 (95% CI 0.49 to 1.24) times less than the outcome incidence odds among the unexposed.

Measures of effect in the exposed:
Exposure changed the outcome incidence risk in the exposed by -0.28 (95% CI -0.78 to 0.23) per 100 population units. -28.5% of outcomes in the exposed were attributable to exposure (95% CI -103.5% to 18.9%).

Number needed to treat for benefit (NNTB) and harm (NNTH):
The number needed to treat for one subject to be harmed (NNTH) is 362 (NNTH 128 to infinity to NNTB 436).

Measures of effect in the population:
Exposure changed the outcome incidence risk in the population by -0.14 (95% CI -0.59 to 0.32) per 100 population units. -12.5% of outcomes in the population were attributable to exposure (95% CI -37.6% to 8.1%).

That’s quite a bit. Another trick is to get at the table of results. We install a package called broom, which includes a number of ways to handle output from R functions, including those in the epiR package. Broom takes from the TidyVerse environment; tables are stored as tibbles.

Note 3: A tibble (package tibble) is a modified R data frame optimized for the tidyverse that enforces strict data handling and cleaner console printing. Unlike standard data frames, tibbles refuse to automatically simplify into vectors when sliced or guess column names using partial matching. Because of these intentional design choices, tibbles often don’t play well with R functions not written explicitly with the TidyVerse environment in mind.

To convert a data.frame to a tibble,

my_tibble <- as_tibble(my_dataframe)

To convert a tibble back to a standard data.frame

my_dataframe <- as.data.frame(my_tibble)

library(broom)

# Test statistics
tidy(fit, parameters = "stat")

R output

# A tibble: 3 × 4
term statistic df p.value
<chr> <dbl> <dbl> <dbl>
1 chi2.strata.uncor 1.15 1 0.284
2 chi2.strata.yates 0.909 1 0.340
3 chi2.strata.fisher NA NA 0.292

We can convert the tibbles into our familiar data.frame format, and then select only the statistics we want.

# Measures of association
fitD <- as.data.frame(tidy(fit, parameters = "moa")); fitD

R output, all 15 measures of association!

                  term     estimate     conf.low    conf.high
1       RR.strata.wald    0.7783605    0.4914679   1.23272564
2     RR.strata.taylor    0.7783605    0.4914679   1.23272564
3      RR.strata.score    0.8742994    0.6584540   1.10340173
4       OR.strata.wald    0.7761915    0.4876209   1.23553616
5     OR.strata.cfield    0.7761915           NA           NA
6      OR.strata.score    0.7761915    0.4894450   1.23093168
7        OR.strata.mle    0.7762234    0.4718655   1.26668220
8    ARisk.strata.wald   -0.2759557   -0.7810162   0.22910484
9   ARisk.strata.score   -0.2759557   -0.8000574   0.23482532
10     NNT.strata.wald -362.3770579 -128.0383246 436.48140194
11    NNT.strata.score -362.3770579 -124.9910314 425.84844829
12  AFRisk.strata.wald   -0.2847517   -1.0347210   0.18878949
13  PARisk.strata.wald   -0.1381661   -0.5933541   0.31702189
14  PARisk.strata.piri   -0.1381661   -0.3910629   0.11473067
15 PAFRisk.strata.wald   -0.1248227   -0.3760279   0.08052298

We can call out just the statistics we want from this table by calling to the specific elements in the data.frame (rows, columns).

fitD[c(1,4,7,9,12),]

R output

                 term   estimate   conf.low  conf.high
1      RR.strata.wald  0.7783605  0.4914679  1.2327256
4      OR.strata.wald  0.7761915  0.4876209  1.2355362
7       OR.strata.mle  0.7762234  0.4718655  1.2666822
9  ARisk.strata.score -0.2759557 -0.8000574  0.2348253
12 AFRisk.strata.wald -0.2847517 -1.0347210  0.1887895

Software: epitools.

Another useful R package for epidemiology is epitools, but it comes with it’s own idiosyncrasies. We have introduced the standard 2X2 format, with a, b, c, and d cells defined as in Table 1, above. However, epitools does it differently, and we need to update the matrix. By default, epitools has the unexposed group (control) in the first row and the non-outcome (no disease) is in the first column. To match our a,b,c, and d matrix, use the epitools command to change this arrangement with the rev() argument. Now, the analysis will use the contingency table on the right where the exposed group (treatment) is in the first row and the outcome (disease) is in the first column (h/t M. Bounthavong 2021). Once that’s accomplished, epitools returns what you would expect.

Calculate relative risk

risk1 <- 32 / (3270 + 32)
risk2 <- 41 / (3525 + 41)
risk1 - risk2

and R returns

-0.00180638

odds ratio

library(epitools)
oddsratio.wald(Table1, rev = c("both"))

and R returns

$data
              Outcome
Predictor Disease2 Disease1 Total
 Exposed2      517       36   553
 Exposed1      518       11   529
Total         1035       47  1082

$measure
           odds ratio with 95% C.I.
Predictor   estimate      lower      upper
 Exposed2  1.0000000        NA         NA
 Exposed1  0.3049657 0.1535563  0.6056675

$p.value
two-sided
Predictor  midp.exact  fisher.exact    chi.square
Exposed2           NA            NA            NA
Exposed1 0.0002954494  0.0003001641  0.0003517007

odds ratio highlighted in green.

Software: OpenEpi.

R is fully capable of delivering the calculations you need, but sometimes you just want a quick answer. Online, the OpenEpi tools at https://www.openepi.com/ can be used for homework problems. For example, working with count data in 2 X 2 format, select Counts > 2X2 table from the side menu to bring up the data form (Fig 1).

Data entry for 2X2 table at openepi.com.

Figure 1. Data entry for 2X2 table at openepi.com.

Once the data are entered, click on the Calculate button to return a suite of results (Fig 2).

Results for 2X2 table at openepi.com

Figure 2. Results for 2X2 table at openepi.com.

Software: RcmdrPlugin.EBM

Note 4: Fall 2023 — I have not been able to run the EBM plugin successfully! The package is still active at CRAN, but looks like it hasn’t been updated since 2015. Simply returns an error message — — on data sets which have in the past performed perfectly. Thus, until further notice, do not use the EBM plugin. Instead, use commands in the epiR package.

This isn’t the place nor can I be the author to discuss what evidence based medicine (EBM) entails (cf. Masic et al. 2008), or what its shortcomings may be (Djulbegovic and Guyatt 2017). Rcmdr has a nice plugin, based on the epiR package, that will calculate ARR, RRR and NNT as well as other statistics. The plugin is called RcmdrPlugin.EBM

install.packages("RcmdrPlugin.EBM", dependencies=TRUE)

After acquiring the package, proceed to install the plug-in. Restart Rcmdr, then select Tools and Rcmdr Plugins (Fig 3).

To install an Rcmdr plugin, first go to Rcmdr → Tools → Load Rcmdr plug-in(s)...

Figure 3. Rcmdr: Tools → Load Rcmdr plugins…

Find the EBM plug-in, then proceed to load the package (Fig 4).

Select the Rcmdr plugin, then click the "OK" button to proceed.

Figure 4. Rcmdr plug-ins available (after first download the files from an R mirror site).

Restart Rcmdr again and the menu “EBM” should be visible in the menu bar. We’re going to enter some data, so choose the Enter two-way table… option in the EBM plug-in (Fig 5)

Select "Enter two-way table...".

Figure 5. R Commander EBM plug-in, enter 2X2 table menus.

To review, we have the following problem, illustrated with natural numbers and probability tree (Fig 6).

probability tree Illustration of the statin problem

Figure 6. Illustration of probability tree for the statin problem.

Now, let’s enter the data into the EBM plugin. For the data above I entered the counts as

	Lived	Died
Statin	468	32
Placebo	459	41

and selected the “Therapy” medical indicator (Fig 7).

EBM plugin with two-way table completed for the statin problem.

Figure 7. EBM plugin with two-way table completed for the statin problem.

The output from EBM plugin was as follows. I’ve added index numbers in brackets so that we can point to the output that is relevant for our worked example here.

(1) .Table <- matrix(c(468,32,459,41), 2, 2, byrow=TRUE, dimnames = list(c('Drug', 'Placebo'), c('Lived', 'Dead')))
(2) fncEBMCrossTab(.table=.Table, .x='', .y='', .ylab='', .xlab='', .percents='none', .chisq='1', .expected='0', .chisqComp='0', .fisher='0', .indicators='th', .decimals=2)

R output begins by repeating the commands used, here marked by lines (1) and (2). The statistics we want follow in the next several lines of output.

(3) Pearson's Chi-squared test data: .Table X-squared = 1.197, df = 1, p-value = 0.2739
(4) # Notations for calculations Event + Event -Treatment "a" "b" Control "c" "d"
(5)# Absolute risk reduction (ARR) = -1.8 (95% CI -5.02 - 1.42) %. Computed using formula: [c / (c + d)] - [a / (a + b)]
(6)# Relative risk = 1.02 (95% CI 0.98 - 1.06) %. Computed using formula: [c / (c + d)] / [a / (a + b)]
(7)# Odds ratio = 1.31 (95% CI 0.81 - 2.11). Computed using formula: (a / b) / (c / d)
(8) # Number needed to treat = -55.56 (95% CI 70.29 - Inf). Computed using formula: 1 / ARR
9)# Relative risk reduction = -1.96 (95% CI -5.57 - 1.53) %. Computed using formula: { [c / (c + d)] - [a / (a + b)] } / [c / (c + d)]
(10)# To find more about the results, and about how confidence intervals were computed, type ?epi.2by2 . The confidence limits for NNT were computed as 1/ARR confidence limits. The confidence limits for RRR were computed as 1 - RR confidence limits.

end R output

In summary, we found no difference between statin and placebo (P-value = 0.2739) and ARR of -1.8%

Questions.

Data from a case-control study on alcohol use and esophageal cancer (Tuyns et al (1977), example from Gerstman 2014). Cases were men diagnosed with esophageal cancer from a region in France. Controls were selected at random from electoral lists from the same geographical region.

	Esophageal cancer
Alcohol grams/day	Cases	Noncases	Total
> 80	96	109	205
< 80	104	666	770
Total	200	775	975

What was the null hypothesis? Be able to write the hypothesis in symbolic form and as a single sentence.
What was the alternative hypothesis? Be able to write the hypothesis in symbolic form and as a single sentence.
What was the observed frequency of subjects with esophageal cancer in this study? And the observed frequency of subjects without esophageal cancer?
Estimate Relative Risk, Absolute Risk, NNT, and Odds ratio?
1. Which is more appropriate, RR or OR? Justify your decision.
The American College of Obstetricians and Gynecologists recommends that women with an average risk of breast cancer (BC) over 40 get an annual mammogram. Nationally, the sensitivity of mammography is about 68% and specificity of mammography is about 75%. Moreover, mammography involves exposure of women to radiation, which is known to cause mutations. Given that the prevalence of BC in women between 40 and 49 is about 0.1%, please evaluate the value of this recommendation by completing your analysis.
A) In this age group, how many women are expected to develop BC?
B) How many False negative would we expect?
C) How many positive mammograms are likely to be true positives?
“Less than 5% of women with screen-detectable cancers have their lives saved,” (quote from BMC Med Inform Decis Mak. 2009 Apr 2;9:18. doi: 10.1186/1472-6947-9-18): Using the information from question 5, What is the Number Needed to Treat for mammography screening?

Quiz Chapter 7.4

Epidemiology: Relative risk and absolute risk, explained

Chapter 7 contents

Probability, Risk Analysis
Epidemiology definitions
Epidemiology basics
Conditional Probability and Evidence Based Medicine
Epidemiology: Relative risk and absolute risk, explained
Odds ratio
Confidence intervals
References and suggested readings

7.1 – Epidemiology definitions

Introduction
Definitions
Questions
Quiz
Chapter 7 contents

Introduction

This subchapter may lack for drama, but let’s start by providing a list of key terms with definitions as you start our introduction to epidemiology. More definitions will follow in the sections as well.

Definitions

Absolute risk: The probability that a specified event will occur in a specified population. See Ch07.4 – Epidemiology: Relative risk and absolute risk, explained.

Absolute risk reduction (ARR): the decrease in risk of an event in an exposed (treatment) group compared to an unexposed (control) group. Also called the risk difference. $CER-EER$ , see Contingency table. See Ch07.4 – Epidemiology: Relative risk and absolute risk, explained.

$CER \ \& \ EER$ : Control Event Rate and Experimental Event Rate. In observational studies, transition from “control” and “experimental” to exposure terms. For example, prospective or retrospective cohort studies, $EER$ is Exposed Event Rate and its partner $UER$ , Unexposed Event Rate. Equivalently, you will encounter $I_{e}$ for Incidence in the Exposed and $I_{u}$ forIncidence in the Unexposed. For Case-Control studies, where the direct event (incidence) rates over time are unavailable, then $EER$ becomes Exposure Rate in Cases — the proportion of people with the condition who were exposed to the risk factor and $CER$ is now Exposure Rate in Controls.

Contingency table, also called cross tabulation or crosstab, is a display of counts of variables in a matrix format. We briefly introduced the contingency table in Chapter 5.1. In epidemiology, rows of contingency table represent treatment or exposure groups, columns represent outcomes.

Table 1. A 2 X 2 contingency table

	Outcome
	Yes	No
Treatment or exposed group	a	b
Control or nonexposed group	c	d

The 2X2 contingency table is referred to frequently in this chapter and again in Chapter 9.

R code

a = 4; b = 46; c = 5; d = 45

Table1 <- matrix(c(a,b,c,d), 2, 2, byrow=TRUE, dimnames = list(c("Treatment", "Control"), c("Yes", "No"))); Table1

R output

          Yes  No
Treatment   4  46
Control     5  45

In statistical classification, a field dedicated to establishing algorithms for correctly identifying and predicting classes of data, eg, correct disease diagnosis or , the 2X2 contingency table is often called a confusion or error matrix (Greene 2003). Machine learning methods are often built on developments from statistical classification findings (Michie et al 1995). We expand on contingency tables in Chapter 7.4 and Chapter 9.2.

Note 1: First discussed in Chapter 5.1, by convention, the standard format for an epidemiological 2×2 contingency table organizes exposure on the rows and disease outcome on the columns as follows:

Rows (Exposure): The top row represents Exposed / Treatment (Yes), and the bottom row represents Non-exposed / Control (No).
Columns (Outcome): The left column represents Disease / Case (Yes), and the right column represents No Disease / Control (No). See Chapter 7.

Control event rate (CER): How often an event occurs in the control group. $\frac{c}{c+d}$ , see Contingency table.

Diagnosis: identification of the nature of a disease or condition.

Event: From probability theory, an event is a set of one or more outcomes to which a probability is assigned to each outcome.

Experimental event rate (EER): How often an event occurs in the treatment group. $\frac{a}{a+b}$ , see Contingency table.

Exposure risk: the potential for harm from contact with or vulnerability to a hazard or risk factor.

Hazard: anything that can cause harm.

Incidence: the number of newly diagnosed individuals in a population having a condition, disease or other characteristic. Compare to prevalence.

Likelihood: In every-day English, the term likelihood and probability mean the same thing — the chance that an event will occur. In statistics, however, probability refers to possible outcomes in the future — a prediction for an event to occur, whereas likelihood is associated with the possible explanation (hypothesis) for the outcome that has already occurred, the past. Put another way, likelihood quantifies how likely is the data we observe, given the hypothesis (cf. excellent discussion by Gallistel 2015)? We will look at likelihood methods in practice in subsequent chapters.

Negative predictive value of a test (NPV), defined as the probability that a negative test result identifies a person who truly does not have the disease. Calculated as the total number of individuals without the disease divided by the total that tested negative. $NPV = \frac{TN}{TN+FN}$

Number needed to treat (NNT): the inverse of the absolute risk reduction. $NNT = \frac{1}{ARR}$ . See Ch07.4 – Epidemiology: Relative risk and absolute risk, explained.

Odds, the ratio (OR) of two probabilities: the probability of getting a one on throwing a dice is $1\6$ , the probability of not getting a one is $5\6$ , therefore the odds of getting a one are 1 to 5. $OR = \frac{a \cdot d}{b \cdot c}$ . See Ch07.5 – Odds ratio.

Outcome, the single result from a random experiment. For example, result of a single coin toss is “heads.” Results of a birth, the outcome is the child born is a singleton, single birth.

Per capita rate, Latin phrase, for each head, meaning per person.

Positive predictive value of a test (PPV), defined as the probability that a positive test result identifies a person who truly has the disease. Calculated as the total number of individuals with the disease divided by the total that tested positive. $PPV = \frac{TP}{TP+FP}$

Posttest probability refers to the probability that the patient has the disease after the results of the test are known.

Pretest probability is the prevalence of the disease, ie, the chance that the a randomly selected person from the population has the disease.

Prevalence: The proportion of individuals in a population having a condition, disease, or characteristic. Compare to incidence.

Probability, defined in Chapter 6.1.

Prognosis, how a disease plays out.

Relative risk: Ratio of the risk of an event among those exposed to the risk factor to the risk among those not exposed to the risk factor. See Ch07.4 – Epidemiology: Relative risk and absolute risk, explained.

Relative risk reduction (RRR): is a measure calculated by dividing the absolute risk reduction by the control event rate. See Ch07.4 – Epidemiology: Relative risk and absolute risk, explained.

Risk: Probability of an event. Risk is not restricted to just bad events, but refers to the uncertainty of a particular event (eg, the risk that a child will be born male seems a melodramatic statement, but it is accurate as far as this definition goes).

Therapy, treatment intended to treat, relieve, or cure a disorder of condition.

Questions

Compare and contrast ARR and RRR.
What’s the difference between event, hazard, and risk?
What’s the difference between incidence and prevalence?
What’s the difference between diagnosis and prognosis?
What’s the implication of a NNT greater than 100 in terms of the utility of a proposed therapy or treatment?

Quiz

Quiz Chapter 7.1

Epidemiology definitions

Chapter 7 contents

Introduction
Epidemiology definitions
Epidemiology basics
Conditional Probability and Evidence Based Medicine
Epidemiology: Relative risk and absolute risk, explained
Odds ratio
Confidence intervals
References and suggested readings