BTE3207 week 6

Minsik Kim

2026 Fall

Before we begin

This week, we are going to compare groups. We will look at the estimate, its confidence interval, and the p-value together.

Run a few chunks as we reach each topic in class. The examples use base R, so we do not need another package.

Prepare our class sample

Let’s use the SBP data again. We will take 1,000 complete observations for today’s practice, instead of testing all one million rows.

dataset_sbp <- read.csv("dataset/sbp_dataset_korea_2013-2014.csv")
class_sbp <- subset(dataset_sbp,
                    !is.na(SBP) & SEX %in% c(1, 2) & DIS %in% 1:4)
set.seed(20260006)
class_sbp <- class_sbp[sample(nrow(class_sbp), 1000), ]

SEX is coded as 1 = male and 2 = female. Keep the original column and make two new variables.

class_sbp$female <- ifelse(class_sbp$SEX == 2, 1, 0)
class_sbp$sex <- factor(class_sbp$SEX, levels = c(1, 2),
                        labels = c("Male", "Female"))
table(class_sbp$SEX, class_sbp$female)

##    
##       0   1
##   1 485   0
##   2   0 515

DIS is disease history: 1 = both hypertension and diabetes, 2 = hypertension only, 3 = diabetes only, and 4 = neither. Here, 1 means a history of hypertension. We are not diagnosing hypertension from today’s SBP.

class_sbp$htn_history <- ifelse(class_sbp$DIS == 1 |
                                class_sbp$DIS == 2, 1, 0)
table(class_sbp$DIS, class_sbp$htn_history)

##    
##       0   1
##   1   0  48
##   2   0 184
##   3  42   0
##   4 726   0

BTH_G is an age category, not age in years. We will not treat its codes as actual ages.

Null hypothesis, p-value, and CI

Suppose our classroom question is whether the population mean SBP differs from 120 mmHg. This is a reference value for practicing a test, not a diagnostic cutoff.

H0: the population mean is 120.
HA: the population mean is different from 120.

one_sample <- t.test(class_sbp$SBP, mu = 120)
one_sample

## 
##  One Sample t-test
## 
## data:  class_sbp$SBP
## t = 3.9359, df = 999, p-value = 8.862e-05
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
##  120.9036 122.7004
## sample estimates:
## mean of x 
##   121.802

Read the estimated mean first, then the 95% CI and p-value. The p-value is the probability, assuming H0 and the test assumptions, of a result at least as extreme as ours. It is not the probability that H0 is true.

one_sample$p.value < 0.05

## [1] TRUE

one_sample$conf.int

## [1] 120.9036 122.7004
## attr(,"conf.level")
## [1] 0.95

With a two-sided test at alpha = 0.05, the matching 95% CI excludes the null value when we reject H0. If we do not reject H0, this does not prove that H0 is true. The CI helps us see which differences remain plausible.

Paired t-test

Same people, two measurements

For a paired test, each before value must match the after value from the same person. The NHIS file does not provide these paired visits, so we will use the 10-woman oral contraceptive example from the lecture (slide 40).

before <- c(115, 112, 107, 119, 115, 138, 126, 105, 104, 115)
after <- c(128, 115, 106, 128, 122, 145, 132, 109, 102, 117)
change <- after - before
data.frame(before, after, change)

Our comparison is after minus before. A positive change means SBP increased.

mean(change)

## [1] 4.8

sd(change)

## [1] 4.565572

sd(change) / sqrt(length(change))

## [1] 1.443761

A paired t-test is a one-sample t-test on the differences. Let’s compare the two commands.

t.test(after, before, paired = TRUE)

## 
##  Paired t-test
## 
## data:  after and before
## t = 3.3247, df = 9, p-value = 0.008874
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  1.533987 8.066013
## sample estimates:
## mean difference 
##             4.8

t.test(change, mu = 0)

## 
##  One Sample t-test
## 
## data:  change
## t = 3.3247, df = 9, p-value = 0.008874
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.533987 8.066013
## sample estimates:
## mean of x 
##       4.8

Voila! The results agree: the mean increase is 4.8 mmHg, with a 95% CI of about 1.5 to 8.1 mmHg. These are the results from the listed observations; use them in place of the slightly different interval printed on the lecture slide.

For this small sample, the paired differences should be approximately normally distributed, without extreme outliers. Different people’s pairs should be independent. A before-after change alone does not show that the contraceptive caused the change; there was no control group.

Try this

Change the order to t.test(before, after, paired = TRUE). What happens to the estimate, CI, and two-sided p-value?

Independent groups: Welch’s t-test

Now we return to class_sbp. Male and female observations are separate groups. We will keep the direction Male minus Female throughout these comparisons.

male_sbp <- class_sbp$SBP[class_sbp$SEX == 1]
female_sbp <- class_sbp$SBP[class_sbp$SEX == 2]
c(Male = mean(male_sbp), Female = mean(female_sbp))

##     Male   Female 
## 124.7443 119.0311

mean(male_sbp) - mean(female_sbp)

## [1] 5.713262

Let’s have a quick look before testing.

boxplot(SBP ~ sex, data = class_sbp,
        xlab = "Sex", ylab = "SBP (mmHg)")

H0 is that the population mean difference is zero. t.test() uses Welch’s method by default: we do not need to assume equal population variances.

welch_result <- t.test(male_sbp, female_sbp)
welch_result

## 
##  Welch Two Sample t-test
## 
## data:  male_sbp and female_sbp
## t = 6.3778, df = 996.34, p-value = 2.746e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.955373 7.471151
## sample estimates:
## mean of x mean of y 
##  124.7443  119.0311

The standard error uses the variability and sample size in each group, not a vector of paired differences.

se_difference <- sqrt(var(male_sbp) / length(male_sbp) +
                      var(female_sbp) / length(female_sbp))
(mean(male_sbp) - mean(female_sbp)) / se_difference

## [1] 6.377769

welch_result$statistic

##        t 
## 6.377769

Do not make pairs just because two groups have the same number of rows. Pairing comes from the study design.

Try this

Use conf.level = 0.99 in the Welch test. Does the CI become wider or narrower? Did the estimated mean difference change?

Binary outcomes: RD, RR, and OR with CIs

Let’s compare the proportions with a history of hypertension, Male compared with Female. These are proportions of recorded history, not new cases followed over time. We use the familiar RD/RR calculations; here their interpretation is a difference or ratio of history prevalence.

htn_table <- table(
  Sex = class_sbp$sex,
  History = factor(class_sbp$htn_history, levels = c(1, 0),
                   labels = c("Yes", "No")))
htn_table

##         History
## Sex      Yes  No
##   Male   114 371
##   Female 118 397

prop.table(htn_table, margin = 1)

##         History
## Sex            Yes        No
##   Male   0.2350515 0.7649485
##   Female 0.2291262 0.7708738

The first column is the event, and the first row is Male. Check the labels before doing any calculation!

n_male <- sum(htn_table["Male", ])
n_female <- sum(htn_table["Female", ])
p_male <- htn_table["Male", "Yes"] / n_male
p_female <- htn_table["Female", "Yes"] / n_female

rd <- p_male - p_female
rr <- p_male / p_female
odds_ratio <- (p_male / (1 - p_male)) / (p_female / (1 - p_female))
c(RD = rd, RR = rr, OR = odds_ratio)

##          RD          RR          OR 
## 0.005925333 1.025860563 1.033806935

Multiply RD by 100 to express it in percentage points. The null value is 0 for RD, and 1 for RR and OR.

CI for the difference

Like the mean difference, the difference in proportions has its own SE. This large-sample CI uses the two separately estimated proportions.

se_rd <- sqrt(p_male * (1 - p_male) / n_male +
              p_female * (1 - p_female) / n_female)
rd + c(-1, 1) * qnorm(0.975) * se_rd

## [1] -0.04643514  0.05828581

CI for a ratio

Remember log() and exp()? We calculate a CI on the log scale and transform its endpoints back. These approximate formulas require sufficiently large cell counts and no zero cells.

The following two chunks are extra calculation practice; focus on the estimate and CI before typing the formulas.

se_log_rr <- sqrt(1 / htn_table["Male", "Yes"] - 1 / n_male +
                  1 / htn_table["Female", "Yes"] - 1 / n_female)
exp(log(rr) + c(-1, 1) * qnorm(0.975) * se_log_rr)

## [1] 0.8187172 1.2854132

se_log_or <- sqrt(sum(1 / htn_table))
exp(log(odds_ratio) + c(-1, 1) * qnorm(0.975) * se_log_or)

## [1] 0.7706735 1.3867828

A ratio CI is usually asymmetric on the original scale. Also, an odds ratio is not generally the same as a prevalence ratio.

Comparing two proportions: z-test

H0 says that the two population proportions are equal. For this test, we estimate their common proportion by pooling both groups under H0.

p_pool <- sum(htn_table[, "Yes"]) / sum(htn_table)
se_null <- sqrt(p_pool * (1 - p_pool) * (1 / n_male + 1 / n_female))
z <- (p_male - p_female) / se_null
p_z <- 2 * pnorm(-abs(z))
c(z = z, p_value = p_z)

##         z   p_value 
## 0.2218516 0.8244294

Notice that se_null is different from the SE we used for the RD confidence interval. We are now computing sampling variability under the null hypothesis.

R also has prop.test(). x gives event counts and n gives group totals, in the same Male, Female order.

proportion_test <- prop.test(x = htn_table[, "Yes"],
                            n = rowSums(htn_table), correct = FALSE)
proportion_test

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  htn_table[, "Yes"] out of rowSums(htn_table)
## X-squared = 0.049218, df = 1, p-value = 0.8244
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04643514  0.05828581
## sample estimates:
##    prop 1    prop 2 
## 0.2350515 0.2291262

prop.test() reports a chi-square statistic. For two groups, its uncorrected two-sided p-value matches the pooled z-test above. We set correct = FALSE to compare the same calculation; R otherwise applies a continuity correction when applicable.

RD = 0, RR = 1, and OR = 1 express the same null hypothesis of equal proportions. This does not mean that every possible test or CI method gives identical numbers.

Chi-square test

We can give R the full table of counts.

chi_result <- chisq.test(htn_table, correct = FALSE)
chi_result

## 
##  Pearson's Chi-squared test
## 
## data:  htn_table
## X-squared = 0.049218, df = 1, p-value = 0.8244

chi_result$expected

##         History
## Sex         Yes     No
##   Male   112.52 372.48
##   Female 119.48 395.52

expected contains the counts expected under independence. Check these counts when deciding whether the chi-square approximation is suitable; the total sample size alone is not enough.

c(z_squared = z^2, chi_squared = unname(chi_result$statistic))

##   z_squared chi_squared 
##  0.04921815  0.04921815

c(z_test = p_z, proportion_test = proportion_test$p.value,
  chi_square = chi_result$p.value)

##          z_test proportion_test      chi_square 
##       0.8244294       0.8244294       0.8244294

The statistics and p-values agree! This equivalence is for the two-sided, pooled two-proportion z-test and the uncorrected Pearson chi-square test on a 2 by 2 table. Chi-square tests can also compare tables with more categories.

More than two categories

Keep all four disease-history categories this time. Is their distribution associated with recorded sex? This test asks about the table as a whole.

dis_table <- table(Sex = class_sbp$sex, History = class_sbp$DIS)
dis_table

##         History
## Sex        1   2   3   4
##   Male    27  87  20 351
##   Female  21  97  22 375

chisq.test(dis_table)

## 
##  Pearson's Chi-squared test
## 
## data:  dis_table
## X-squared = 1.2833, df = 3, p-value = 0.7331

chisq.test(dis_table)$expected

##         History
## Sex          1     2     3      4
##   Male   23.28 89.24 20.37 352.11
##   Female 24.72 94.76 21.63 373.89

Fisher’s exact test

What if some expected counts are small? Fisher’s exact test is useful for a small or sparse 2 by 2 table. It can also be used with our class sample.

fisher_result <- fisher.test(htn_table)
fisher_result

## 
##  Fisher's Exact Test for Count Data
## 
## data:  htn_table
## p-value = 0.8809
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7621301 1.4019413
## sample estimates:
## odds ratio 
##   1.033786

c(chi_square = chi_result$p.value, Fisher = fisher_result$p.value)

## chi_square     Fisher 
##  0.8244294  0.8808693

Fisher’s test conditions on the table margins. Its p-value need not equal the chi-square p-value, even if they round to the same displayed number. The reported OR is a conditional estimate, so it may also differ slightly from our direct calculation.

Try this: a smaller SBP sample

The first 30 rows of our already randomized class sample give us a smaller table. Keep the same Male/Female order and Yes/No labels.

small_data <- class_sbp[1:30, ]
small_table <- table(
  Sex = small_data$sex,
  History = factor(small_data$htn_history,
                   levels = c(1, 0), labels = c("Yes", "No")))
small_table

##         History
## Sex      Yes No
##   Male     3 12
##   Female   3 12

fisher.test(small_table)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  small_table
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.1103964 9.0582665
## sample estimates:
## odds ratio 
##          1

Read the OR, its CI, and the p-value. Does a small table give a precise estimate? A small p-value and a narrow CI are different ideas.

Two-sided and one-sided tests

Our tests so far asked whether the groups differ in either direction. A one-sided test asks a directional question. For t.test(male_sbp, female_sbp), "greater" means that the male population mean is greater than the female population mean.

p_greater <- t.test(male_sbp, female_sbp,
                    alternative = "greater")$p.value
p_less <- t.test(male_sbp, female_sbp, alternative = "less")$p.value
c(two_sided = welch_result$p.value,
  greater = p_greater, less = p_less)

##    two_sided      greater         less 
## 2.745681e-10 1.372840e-10 1.000000e+00

Choose the direction before looking at the results. For this continuous, symmetric t-test, the one-sided p-value is half the two-sided value only when the observed difference is in the specified direction. It is not a rule to apply automatically to every test.

Errors and interpretation

Our decision can be wrong even when we use the test correctly.

Reality	Reject H0	Do not reject H0
H0 is true	Type I error	Correct decision
H0 is false	Correct decision	Type II error

Alpha controls the Type I error probability under the test assumptions. For a specified alternative, beta is the Type II error probability and power is 1 - beta. Power depends on the effect size, variability, sample size, and test.

Before reporting a result, ask:

How large is the difference, and what does its CI tell us?
Are observations independent, and do the assumptions fit the study design?
Who does this sample represent? Sampling 1,000 rows does not make the original NHIS file representative of every population.
Could age or another variable explain part of the association?

A small p-value does not prove causation or practical importance. A large p-value does not prove that the groups are identical. Our simple comparisons do not adjust for confounding or account for a complex survey design.

Function details: t.test, prop.test, chisq.test, fisher.test.

Bibliography

R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Xie Y (2025). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.50, https://yihui.org/knitr/.

Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/.

Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.

Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2025). rmarkdown: Dynamic Documents for R. R package version 2.30, https://github.com/rstudio/rmarkdown.

Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown.

Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.

Barnier J (2022). rmdformats: HTML Output Formats and Templates for ‘rmarkdown’ Documents. R package version 1.0.4, https://CRAN.R-project.org/package=rmdformats.