Statistics on categorical variables

2026 fall

Hypothesis Testing Recap

Hypothesis Testing Recap

  1. Specify your outcome and population parameter of interest
    • Figure out which mathematical distribution best represents your population data.
      • which parameter of this distribution are you interested in?
  2. Set up your null and alternative hypotheses
    • What would the population parameters be under the conditions of the null hypothesis?
  3. Collect data
  4. Choose an appropriate statistical test; assess limitations and assumptions
    • how can you use the sample observations to estimate the parameters of the population distribution?
      • point estimator
        • Use calculus to derive the mathematical distribution of your sample point estimator.
        • use the tried and true methods
    • are all the prerequisite conditions for using this test satisfied?
  5. Compute likelihood of observed result assuming the null hypothesis is true
    • If you know the distribution of the point estimator, you can calculate probabilities mathematically.
      • how likely is it to see an observation as extreme as what we’ve seen in our sample?
        • p-value
  6. Conclusion
    • If the p-value < the test’s significance level \(\alpha\) then reject \(H_0\) in favor of \(H_1\).

Hypothesis Test Recap (example #1)

Q: A school nurse believes the average weight of students is normally distributed and has a mean greater than 100 pounds. A random sample of 30 students has a mean weight of 105 pounds. Assume the population standard deviation is known and equals 15 pounds.

At the 5% significance level, test whether the average student weight is greater than 100 pounds.

A:

1. Specify your outcome and population parameter of interest

2. Set up your null and alternative hypotheses

We want to see if the average student weight is greater than 100. This is the alternative hypothesis. The null hypothesis is that average is not greater than 100. We’re trying to find evidence against the null in favor of the alternative.

3. Collect data

We have been provided the mean weights from a sample of 30 students: 105 lbs.

4. Choose an appropriate statistical test; assess limitations and assumptions

5. Compute likelihood of observed result assuming the null hypothesis is true

Here’s what the distribution of the sample mean looks like under our null hypothesis.

# Parameters
mu <- 100
sigma <- 15 / sqrt(30)

# Plot the normal curve
curve(dnorm(x, mean = mu, sd = sigma), from = 90, to = 110, 
      xlab = "x", ylab = "Density", main = "Normal Distribution with Shaded Area (x > 105)")

# Shade area to the right of 105
x_shade <- seq(105, 110, length = 100)
y_shade <- dnorm(x_shade, mean = mu, sd = sigma)
polygon(c(105, x_shade, 110), c(0, y_shade, 0), col = "skyblue", border = NA)

# Add vertical line at 105
abline(v = 105, col = "red", lty = 2)

#we want the area of the curve to the right of the critical value
#this is equal to 1 - the area to the left of the critical value
1 - pnorm(105, mean=100, sd=15/sqrt(30))
## [1] 0.03394458
#equivalently, we can transform the sample mean and compare it to the standard normal distribution
z = (105- 100)/(15/sqrt(30))
1 - pnorm(z, mean=0, sd=1)
## [1] 0.03394458

6. Conclusion

The p-value is less than 0.05 so we reject the null in favor of the alternative.

.

Categorical variables

The focus of this lecture is how to deal with nominal qualitative data, specifically categorical variables that have no meaningful ordering.

.

Statistical Terms

Many hypothesis tests on categorical variable uses the binomial distribution to model population parameters. Here’s what the distribution looks like under different values of n and p.

# 4x4 grid of binomial pmfs

# parameter grids
ns <- c(10, 20, 50, 100)
ps <- c(0.1, 0.3, 0.5, 0.7)

# prepare plotting area
par(mfrow = c(4, 4), mar = c(2, 2, 2, 1), oma = c(4, 4, 4, 2))

# iterate rows: n values
for (n in ns) {
  # iterate columns: p values
  for (p in ps) {
    x <- 0:n
    y <- dbinom(x, size = n, prob = p)
    plot(x, y, type = "h",
         lwd = 2, col = "steelblue",
         main = paste0("n=", n, " p=", p),
         xlab = "", ylab = "")
  }
}

# add common axis labels
mtext("Number of successes (x)", side = 1, outer = TRUE, line = 2)
mtext("Probability P(X = x)", side = 2, outer = TRUE, line = 2)
mtext("Binomial PMFs for selected n & p", side = 3, outer = TRUE, line = 1.5, cex = 1.2)

.

Tests on One Categorical Variable

Dealing with One Categorical Variable(exact).

Answer:

The boils down to a question on the counts appearing in two categories.
Did the patient survive 5 years? (Yes/No).

First think about our point estimate for the survival rate: 6/52 = 0.1153846.
If the 5 year survival rate really was 20%, then we’d expect this proportion to be close to 0.20. The proportion we see is less than 0.20, but is this due to the random fluctuations in our sampling procedure, or should the survival rate actually be lower?

We set up a hypothesis test to find out:

We can think of this in terms of coin flips

x = 0:52
y <- dbinom(x, size = 52, prob = 0.2)
plot(x, y, type = "h", lwd = 2)
points(0:6, y[1:7], col="red", type = "h", lwd = 2)

#sum of the probabilties of all the values that are as or more extreme than what we've observed
sum(dbinom(0:6, size=52, prob=0.20))
## [1] 0.08228342
#cumulative distribution function
pbinom(6, size=52, prob=0.20)
## [1] 0.08228342
#use the built-in test
binom.test(6, n=52,p=0.20, alternative="less")
## 
##  Exact binomial test
## 
## data:  6 and 52
## number of successes = 6, number of trials = 52, p-value = 0.08228
## alternative hypothesis: true probability of success is less than 0.2
## 95 percent confidence interval:
##  0.0000000 0.2150871
## sample estimates:
## probability of success 
##              0.1153846

The pval is 0.08228342, so we would not reject the null hypothesis. We can’t conclude that the actual survival rate is lower than 20% based on this data alone.

You’ll notice the built-in test also has 95% confidence intervals for the point estimate \(\hat{p} = \frac{x}{n}\) This interval contains the survival rate we were checking, 0.20, so there’s no evidence to reject the null here either.

.

Dealing with One Categorical Variable (approximate).

When n is large, the distribution of the sample proportion \(\hat p\) can be approximated with a Normal distribution:
\(N(\mu = p, \sigma=\sqrt{p(1-p)/n})\)

But this is not a very good approximation unless:

When n is large enough (#successes & #failures are both more than 5), we can set up a hypothesis test using the approximate normal distribution.

#run the hypothesis test
phat = 6/52
p0 = 0.20
n = 52

#compute the value
z = (phat - p0)/(sqrt(p0*(1-p0)/n))
z
## [1] -1.525426
#whats the probability of seeing a value more extreme than what we've observed?
pval = pnorm(z, mean=0, sd=1)
pval
## [1] 0.0635765
#or you can just use R's built-in function
prop.test(6, 52,p=.20, alternative="less", correct = F)
## 
##  1-sample proportions test without continuity correction
## 
## data:  6 out of 52, null probability 0.2
## X-squared = 2.3269, df = 1, p-value = 0.06358
## alternative hypothesis: true p is less than 0.2
## 95 percent confidence interval:
##  0.0000000 0.2079583
## sample estimates:
##         p 
## 0.1153846
#compute the wald intervals
alpha = 0.05
qhat = 1-phat
z2 = qnorm(alpha/2,mean=0, sd=1)
lower = phat + z2 *sqrt(phat*qhat/n)
upper = phat - z2 *sqrt(phat*qhat/n)
lower
## [1] 0.02854905
upper
## [1] 0.2022202

Some notes on using R.

Continuity Correction

By default, R does a procedure called “Yate’s continuity correction”, to address the issue of low counts in a categorical comparison. This procedure artificially adds fractional counts to the cells of a contingency table in order to get a better approximation to the normal distribution. We can disable continuity correction by adding the , correct=F parameter to the prop.test function call. If we don’t do this, we’ll get slightly different results from what we computed above.

#if we enable continuity correct, we won't get the same result we did by hand
prop.test(6, 52,p=.20, alternative="less")
## 
##  1-sample proportions test with continuity correction
## 
## data:  6 out of 52, null probability 0.2
## X-squared = 1.8281, df = 1, p-value = 0.08817
## alternative hypothesis: true p is less than 0.2
## 95 percent confidence interval:
##  0.0000000 0.2193964
## sample estimates:
##         p 
## 0.1153846

CI methods

You’ll also notice that the confidence intervals from prop.test are slightly different from the values we calculated. This is because prop.test uses a different procedure called the “Wilson method” which is more accurate but harder to implement than the “Wald method” we ran ourselves. Base R doesn’t have a built in method to do Wald CI’s for proportions but the “binom” package has it implemented, and we get matching results to what we computed ourselves.

library(binom)
#in this package, wald intervals are called "asymptotic"
binom.confint(x = 6, n = 52, methods = "asymptotic")
##       method x  n      mean      lower     upper
## 1 asymptotic 6 52 0.1153846 0.02854905 0.2022202

1-sided vs 2-sided CIs

Though two sided CIs are much more commonly reported, it is also possible to calculate one sided versions.

For the Wald method, this would involve using \(z_{\alpha}\) instead of \(z_{\alpha/2}\) in the interval calculation

And since it is a one sided interval, the range only applies to a single direction
(the other direction for the proportion is either a 0 or 1).

To find the one sided 95% Wald confidence interval corresponding to our 1 sided hypothesis test:

alpha = 0.05
z1 = qnorm(alpha,mean=0, sd=1)
lower = phat + z1 *sqrt(phat*qhat/52)
upper = phat - z1 *sqrt(phat*qhat/52)

lower
## [1] 0.04250991
upper
## [1] 0.1882593

The one sided Wald CI depends on the direction you are interested in in the hypothesis test. It will be one of the following:

Notice that the parameter under the null hypothesis (\(H_0: p >= 0.20\)) lies outside the range of the one-sided Wald CI [0, 0.1882593]!
This gives us seemingly conflicting results compared to the one-sided hypothesis tests we ran above.

How can we explain and resolve this?

Summary For Tests of Proportion on 1 Sample

.

Tests on Two Categorical Variables

Dealing with Two Categorical Variables (associations)

When we have more than one categorical variable, we’re often interested in the association between them.

In the example below, we want to know if the chances of a patient progressing to a disease is related to whether they received either a placebo or drug treatment. In this case, we have categorical labels (ctrl/drug) for each patient telling us which cohort they belong to and an indicator (y/n) of whether or not they got the disease.

How would you check for association?

table(x$cohort)
## 
## ctrl drug 
##   58   42
table(x$disease)
## 
##  no yes 
##  51  49
tab1 = table(x$disease, x$cohort)
tab1
##      
##       ctrl drug
##   no    26   25
##   yes   32   17
addmargins(tab1)
##      
##       ctrl drug Sum
##   no    26   25  51
##   yes   32   17  49
##   Sum   58   42 100

.

Statistical Terms

Notation

  • proportions (risk)
    • \(p_1 = \frac{a}{a+c}\): proportion of disease from the exposed
    • \(p_2 = \frac{b}{b+d}\): proportion of disease from the non exposed
  • relative risk (RR)
    • \(p_1/p_2\)
  • risk difference (RD)
    • \(p_1 - p_2\)
  • odds
    • \(\frac{p_1}{1-p_1} = \frac{a}{c}\)
    • \(\frac{p_2}{1-p_2} = \frac{b}{d}\)
  • odds ratio (OR)
    • \(\frac{p_1/(1-p_1)}{p_2/(1-p_2} = \frac{\frac{a}{c}}{\frac{b}{d}} = a d / b c\)

Confidence Intervals

If all cells in the 2 by 2 tables are at least 5, we can use the central limit theorem to estimate confidence intervals.

  • Notes:
    • \(p_1\) is the proportion in the first population
    • \(\hat{p_1}\) is the sample proportion from the first sample
    • \(\hat{q_1} = 1 - \hat{p_1}\)
    • if \(\alpha\) = 0.05, 100 (1 - \(\alpha\))% CI is another way of saying 95% confidence interval.
  • Relative Risk(RR)
    • 100 (1 - \(\alpha\))% CI for log(\(p_1/p_2\)) is
    • \(\text{log}(\hat{p_1} / \hat{p_2}) \pm z_{\alpha/2} \sqrt{\frac{\hat{q_1}}{n_1 \hat{p_1}} + \frac{\hat{q_2}}{n_2 \hat{p_2}}}\)
  • Risk Difference(RD)
    • 100 (1 - \(\alpha\))% CI for \(p_1 - p_2\) is
    • \((\hat{p_1} - \hat{p_2}) \pm z_{\alpha/2} \sqrt{\frac{\hat{p_1}\hat{q_2}}{n_1} + \frac{\hat{p_2}\hat{q_2}}{n_2}}\)
  • Odds Ratio(OR)
    • 100 (1 - \(\alpha\))% CI for log(OR) is
    • \(\text{log}(\frac{ad}{bc}) \pm z_{\alpha/2}\sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d} }\)

Important Notes

  • In cohort studies:
    • Start with a group classified by exposure status (placebo vs drug)
    • Then follow them forward to observe whether they develop the outcome (disease)
    • RR, RD, OR can all be estimated
  • What about case-control studies?
    • Start by selecting based on outcome (cases = diseased, controls = not diseased)
    • Then look backward to see who was exposed
    • Because you’ve artificially fixed the number of cases and controls, you can’t estimate the actual risk (or incidence) of disease in exposed vs unexposed.
    • Only OR can be estimated!
  • OR doesn’t change when rows and cols are switched and you can always use it to check for associations.
##      
##       ctrl drug Sum
##   no    26   25  51
##   yes   32   17  49
##   Sum   58   42 100

.

Statistical Terms

There are R packages that can find all these numbers for you, but you have to be sure you understand how they represent the conditions and groupings. Transposing the rows/cols can give completely different meanings to the summary statistics and results. If you’re going to use these packages, it’s always a good idea to compute some of the summary stats yourself to make sure they match up with what’s being reported by the package functions.

library("epiR")
print(tab1)
##      
##       ctrl drug
##   no    26   25
##   yes   32   17
epi.2by2(tab1, method = "cohort.count")
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+          26          25         51     50.98 (36.60 to 65.25)
## Exposure-          32          17         49     65.31 (50.36 to 78.33)
## Total              58          42        100     58.00 (47.71 to 67.80)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 0.78 (0.56, 1.09)
## Inc odds ratio                                 0.55 (0.25, 1.24)
## Attrib risk in the exposed *                   -14.33 (-33.45, 4.80)
## Attrib fraction in the exposed (%)            -28.10 (-82.00, 8.39)
## Attrib risk in the population *                -7.31 (-23.77, 9.16)
## Attrib fraction in the population (%)         -12.60 (-15.52, -5.55)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 2.105 Pr>chi2 = 0.147
## Fisher exact test that OR = 1: Pr>chi2 = 0.162
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units
print(t(tab1))
##       
##        no yes
##   ctrl 26  32
##   drug 25  17
epi.2by2(t(tab1), method = "cohort.count")
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+          26          32         58     44.83 (31.74 to 58.46)
## Exposure-          25          17         42     59.52 (43.28 to 74.37)
## Total              51          49        100     51.00 (40.80 to 61.14)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 0.75 (0.52, 1.10)
## Inc odds ratio                                 0.55 (0.25, 1.24)
## Attrib risk in the exposed *                   -14.70 (-34.30, 4.90)
## Attrib fraction in the exposed (%)            -32.78 (-94.62, 9.70)
## Attrib risk in the population *                -8.52 (-26.31, 9.26)
## Attrib fraction in the population (%)         -16.71 (-21.65, -6.07)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 2.105 Pr>chi2 = 0.147
## Fisher exact test that OR = 1: Pr>chi2 = 0.162
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units
tab2 = tab1[2:1,2:1]
print(tab2)
##      
##       drug ctrl
##   yes   17   32
##   no    25   26
epi.2by2(tab2, method = "cohort.count")
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+          17          32         49     34.69 (21.67 to 49.64)
## Exposure-          25          26         51     49.02 (34.75 to 63.40)
## Total              42          58        100     42.00 (32.20 to 52.29)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 0.71 (0.44, 1.14)
## Inc odds ratio                                 0.55 (0.25, 1.24)
## Attrib risk in the exposed *                   -14.33 (-33.45, 4.80)
## Attrib fraction in the exposed (%)            -41.29 (-129.69, 11.21)
## Attrib risk in the population *                -7.02 (-23.81, 9.77)
## Attrib fraction in the population (%)         -16.71 (-21.26, -7.93)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 2.105 Pr>chi2 = 0.147
## Fisher exact test that OR = 1: Pr>chi2 = 0.162
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units
print(t(tab2))
##       
##        yes no
##   drug  17 25
##   ctrl  32 26
epi.2by2(t(tab2), method = "cohort.count")
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+          17          25         42     40.48 (25.63 to 56.72)
## Exposure-          32          26         58     55.17 (41.54 to 68.26)
## Total              49          51        100     49.00 (38.86 to 59.20)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 0.73 (0.48, 1.13)
## Inc odds ratio                                 0.55 (0.25, 1.24)
## Attrib risk in the exposed *                   -14.70 (-34.30, 4.90)
## Attrib fraction in the exposed (%)            -36.31 (-114.82, 9.92)
## Attrib risk in the population *                -6.17 (-22.29, 9.95)
## Attrib fraction in the population (%)         -12.60 (-15.30, -6.89)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 2.105 Pr>chi2 = 0.147
## Fisher exact test that OR = 1: Pr>chi2 = 0.162
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

This last way of looking at the data is what we want. The columns are ordered by Outcome(+ then -), and the rows are ordered by exposure (+ then -).

“epiR” package Terminology

  • Inc risk: incidence risk. The risk of getting a ‘+’ outcome given an exposure level. It’s reported as a % ranging from 0 to 100.
    • 40.48 comes from 17/(17+25) * 100.
    • 55.17 comes from 32 / (32+26) * 100.
  • Inc risk ratio: Risk Ratio or Relative Risk(RR).
    • 0.73 comes from 40.48 / 55.17
    • RR > 1 means exposure is associated with higher risk.
    • RR < 1 means protective effect.
  • Inc odds ratio: odds ratio (OR)
    • 0.55 comes from (17 * 26) / ( 25 * 32)
    • Odds of outcome in exposed vs. unexposed.
    • OR = 1 means no association (odds are the same in both groups)
    • OR > 1 means odds of the outcome are higher in the exposed group
    • OR < 1 means odds of the outcome are lower in the exposed group (possible protective effect)
    • For cohort studies of rare diseases, OR and RR are usually close.
  • Attrib risk in the exposed *: Risk difference (RD)
    • tells you how much of the outcome among the exposed group is due to the exposure (assuming a causal relationship)
    • -14.70 comes from 40.48 - 55.17
print(t(tab2))
##       
##        yes no
##   drug  17 25
##   ctrl  32 26
epi.2by2(t(tab2), method = "cohort.count")
##              Outcome+    Outcome-      Total                 Inc risk *
## Exposure+          17          25         42     40.48 (25.63 to 56.72)
## Exposure-          32          26         58     55.17 (41.54 to 68.26)
## Total              49          51        100     49.00 (38.86 to 59.20)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 0.73 (0.48, 1.13)
## Inc odds ratio                                 0.55 (0.25, 1.24)
## Attrib risk in the exposed *                   -14.70 (-34.30, 4.90)
## Attrib fraction in the exposed (%)            -36.31 (-114.82, 9.92)
## Attrib risk in the population *                -6.17 (-22.29, 9.95)
## Attrib fraction in the population (%)         -12.60 (-15.30, -6.89)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 2.105 Pr>chi2 = 0.147
## Fisher exact test that OR = 1: Pr>chi2 = 0.162
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

.

What to look for in terms of association, proportions and odds

##       
##        yes no
##   drug  17 25
##   ctrl  32 26

Testing for equality of Proportions

In our example, we want to see if the risk of getting the disease(yes/no) is different in the two patient groups(ctrl/drug). We’re checking for association between disease and cohort using the proportions.

addmargins(tab3)
##       
##        yes  no Sum
##   drug  17  25  42
##   ctrl  32  26  58
##   Sum   49  51 100

Set up the hypothesis test.

row_totals <- rowSums(tab3)
col_totals <- colSums(tab3)
grand_total <- sum(tab3)

# Outer product to get expected counts
tab3.exp <- outer(row_totals, col_totals) / grand_total
addmargins(tab3.exp)
##        yes    no Sum
## drug 20.58 21.42  42
## ctrl 28.42 29.58  58
## Sum  49.00 51.00 100
stat = sum((tab3 - tab3.exp)^2/ tab3.exp)
stat
## [1] 2.105341
# Set up values
df <- 1
x <- seq(0, 6, length.out = 500)
y <- dchisq(x, df)

# Plot chi-squared density
plot(x, y, type = "l", lwd = 2, col = "steelblue",
     ylab = "Density", xlab = expression(chi^2),
     main = expression(paste("Chi-squared Distribution with 1 df")))

# Shade the area to the right of stat
x_shade <- seq(stat, max(x), length.out = 200)
y_shade <- dchisq(x_shade, df)
polygon(c(stat, x_shade, max(x)), c(0, y_shade, 0), col = "orange", border = NA)

# Add vertical line for threshold
abline(v = stat, col = "red", lty = 2)
text(stat, dchisq(stat, df) + 0.02, paste0("Stat = ", stat), pos = 4, col = "red")

#we want the right hand tail of the distribution
1 - pchisq(stat, df=1)
## [1] 0.1467856

In the past, we had to do all this by hand on paper with calculators and pdf lookup tables.
Just now, to illustrate the procedure, we did everything ourselves in R. In practice, this can all be done with a built-in R function.

#turn off continuity correction
chisq.test(tab3, correct=F)
## 
##  Pearson's Chi-squared test
## 
## data:  tab3
## X-squared = 2.1053, df = 1, p-value = 0.1468

We get the same answer doing it ‘by hand’ or using the built-in test.

.

Testing for Odds ratio != 1

Another way to test for association between disease and cohort in our example dataset is to test if the odds ratio close to 1 using a Fisher’s test.
Testing on the odds can be done in case-control studies where it’s not possible to observe the true risks. It’s an exact test and does not have any requirements on the number of success/failures like some of the other tests we covered.

fixed margins

If you assume the row/col counts of a 2x2 table are fixed, then all you need is a single cell number, and you can fill in the rest of the values.

##       
##        yes no Sum
##   drug 17  x  42 
##   ctrl x   x  58 
##   Sum  49  51 100

This scenario can be modeled with a hypergeometric distribution.

hypergeometric distribution

How do we plug in the data from our example into the hypergeometric distribution to calculate probabilities?

An urn contains 49 white balls, and 51 black balls. If we random scooped up 42 balls, what’s the probability of observing 17 white?

#parameters

m = 49
n = 51
k = 42
x = 17

#probability
dhyper(x, m,n,k)
## [1] 0.05702839

How do we use this to run a hypothesis test?

  1. Set up the null. Assume there’s no signal
    • the margins of the contingency table are fixed
    • there’s no association between rows and columns
    • the counts of the rest of the table are filled in randomly.
  2. calculate the probabilty of seeing a sample as extreme as what we’ve observed in our data
    • what are the chances that 17 or fewer of the 42 drug patients got the disease?
    • what are the chances of seeing 17 or fewer white balls in a scoop of 42 balls?
  3. If this probability is small, reject the null.
    • if pval < 0.05 reject \(H_0\)
#visualize the distribution
x_max <- min(k, m)  # maximum possible successes

# Probability mass function
x_vals <- 0:x_max
probs <- dhyper(x_vals, m = m, n = n, k = k)

# Color bars: blue if x <= 17, gray otherwise
colors <- ifelse(x_vals <= 17, "skyblue", "lightgray")

# Plot
barplot(probs, names.arg = x_vals, col = colors,
        main = "P(X ≤ 17) in Hypergeometric(m=49, n=51, k=42)",
        xlab = "x", ylab = "Probability")

#calculate the area under the curve
probs = dhyper(0:x, m,n,k)
probs
##  [1] 1.076589e-19 2.215620e-17 1.981973e-15 1.035030e-13 3.570855e-12
##  [6] 8.723088e-11 1.577910e-09 2.180897e-08 2.357293e-07 2.028436e-06
## [11] 1.409229e-05 7.994172e-05 3.736958e-04 1.450358e-03 4.702403e-03
## [16] 1.280098e-02 2.937826e-02 5.702839e-02
#this is the area under the curve
sum(probs)
## [1] 0.1058304
#or equivalently we can use the cumulative distribution function
phyper(x, m,n,k)
## [1] 0.1058304

We can do all this ourselves, or we use the built-in tools

tab3
##       
##        yes no
##   drug  17 25
##   ctrl  32 26
fisher.test(tab3, alternative="less")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab3
## p-value = 0.1058
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##  0.000000 1.171882
## sample estimates:
## odds ratio 
##  0.5558203

We get the same answer.

important note

the ‘epi.2by2’ function in the ‘epiR’ package can only run a two sided fisher test. If you want to run a one sided test you’ll have use the builtin fisher.test as above.

.

Summary For Tests on Two Categorical Variables

Tests on Two Paired Categorical Variables

Dealing with Two Paired Categorical Variables (format)

As an example, we wanted to test two different foot cream treatments to see which one had better results at clearing up an infection.

wide format

If your data is set up so that you have 1 row for each patient, and 2 columns for each of the conditions of the paired experiment, you don’t need to do anything extra. The data is already in a paired format and you can make the paired contingency table directly.

##          fungacream
## pedacream  0  1
##         0 82 16
##         1 37  9

long format

If your data has multiple rows for each patient, each row a different condition of the paired experiment, then you have to reformat it to be able to create the paired contingency tables.

htmltools::div(
    style = "width: 400px; text-align: left;",
    datatable(y.foot, rownames=F)
)
#this is an unpaired contingency table!  It's not the right way to handle this data
table(treatment=y.foot$Treatment, result=y.foot$Result)
##             result
## treatment      0   1
##   fungacream 119  25
##   pedacream   98  46
library("tidyr")
#first get the data into wide format
x.foot2 = pivot_wider(y.foot, names_from = Treatment, values_from = Result)
htmltools::div(
    style = "width: 400px; text-align: left;",
    datatable(x.foot2, rownames=F)
)
#now you can make the paired contingency table
table(pedacream=x.foot2$pedacream, fungacream=x.foot2$fungacream)
##          fungacream
## pedacream  0  1
##         0 82 16
##         1 37  9

.

Dealing with Two Paired Categorical Variables (exact)

Q: Does the treatment(pedacream/fungicream) have an effect on the outcome(0/1)?

Is there a significant difference in paired proportions?

How do we set up a hypothesis test to check this?

addmargins(table(pedacream=x.foot$pedacream, fungacream=x.foot$fungacream))
##          fungacream
## pedacream   0   1 Sum
##       0    82  16  98
##       1    37   9  46
##       Sum 119  25 144
  1. Set up the null. Assume there’s no signal
    • the rate for the discordant pairs should be the same
      • an even coin flip
    • we tossed a coin r+s times, and saw 16 heads and 37 tails.
    • binomial(x, n=r+s, p=0.5)
  2. calculate the probabilty of seeing a sample as extreme as what we’ve observed in our data
    • If you toss a coin r + s times, whats the chance of seeing 16 or fewer heads or 37 or greater heads?
      • we need to check both directions
  3. If this probability is small, reject the null.
    • if pval < 0.05 reject \(H_0\)
# Critical values
s <- 16
r <- 37
n = s + r
x <- 0:n
prob <- dbinom(x, size = n, prob = 0.5)


#area of the left tail
a1 = pbinom(s, size=n, prob=0.5)
#area of the right tail.  Pay attention to the tricky calculation
a2 = 1 - pbinom(r-1, size=n, prob=0.5)
#pval
pval = a1+a2
pval 
## [1] 0.005486345
# Plot
barplot(prob, names.arg = x, col = ifelse(x <= s | x >= r, "red", "gray"),
        main = paste0("Binomial(n = 53, p = 0.5)\nTwo-tailed p-value = ", round(pval, 7)),
        xlab = "Number of Successes", ylab = "Probability")

or you can use the built-in test, and ignore everything except the p-value

binom.test(16, 16+37)
## 
##  Exact binomial test
## 
## data:  16 and 16 + 37
## number of successes = 16, number of trials = 53, p-value = 0.005486
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.183402 0.443434
## sample estimates:
## probability of success 
##              0.3018868

We get the same result.

.

Dealing with Two Paired Categorical Variables (approx)

Another way to see if there’s a difference in the treatments with with a McNemar’s test. This tests whether the proportion of discordant pairs is different from 0.50 by looking at the minimum(r,s). It requires a large sample(both r & s must be at least 5) and uses a normal approximation to the sample proportion.

phat = 16/(16+37)
z = (phat-.50)/(sqrt(.5*.5/53))
z
## [1] -2.884572
pval = 2*pnorm(z)
pval
## [1] 0.003919463
tab2 = table(pedacream=x.foot$pedacream, fungacream=x.foot$fungacream)
mcnemar.test(tab2, correct=F)
## 
##  McNemar's Chi-squared test
## 
## data:  tab2
## McNemar's chi-squared = 8.3208, df = 1, p-value = 0.003919

.

Summary For Tests on Two Paired Categorical Variables

Which test to use?

.

Next Lecture