Hypothesis Testing in R for Health and Social Science

Multiple worked examples of three t-tests, one-way ANOVA, and chi-square tests using NHANES and social-science data

Why hypothesis testing matters

Health and social-science research often moves from description to comparison. We may ask whether average systolic blood pressure differs between groups, whether body mass index differs from a clinically meaningful reference value, whether average sleep differs across demographic groups, or whether smoking status is associated with sex or education.

This chapter uses real health and social-science sample data rather than invented observations. The main dataset is NHANES from the R package of the same name. It contains variables derived from the US National Health and Nutrition Examination Survey. We also use MASS::survey, a student survey dataset that is useful for social-science teaching.

The purpose is not only to run tests. For every method we will:

1. state a realistic research question;
2. write the null and alternative hypotheses;
3. inspect the data;
4. calculate descriptive statistics;
5. check the assumptions that matter;
6. run the test in R;
7. extract useful results;
8. discuss effect size where appropriate;
9. interpret the result in plain language.

Important: The examples below are conventional unweighted teaching analyses. NHANES has a complex survey design. For substantive population inference from NHANES, analysts should use the appropriate sampling weights, strata, clusters, eligible subsamples and survey-design methods. The simple tests here are intended to teach the mechanics and interpretation of hypothesis testing.

Setup

Install these packages once if needed.

install.packages(c("NHANES", "dplyr", "tidyr", "ggplot2", "knitr", "broom", "effectsize"))

library(NHANES)
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)

options(dplyr.summarise.inform = FALSE)

Load the datasets.

data(NHANES, package = "NHANES")
data(survey, package = "MASS")

glimpse(NHANES)

## Rows: 10,000
## Columns: 76
## $ ID               <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646, 5164…
## $ SurveyYr         <fct> 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 2009_10,…
## $ Gender           <fct> male, male, male, male, female, male, male, female, f…
## $ Age              <int> 34, 34, 34, 4, 49, 9, 8, 45, 45, 45, 66, 58, 54, 10, …
## $ AgeDecade        <fct>  30-39,  30-39,  30-39,  0-9,  40-49,  0-9,  0-9,  40…
## $ AgeMonths        <int> 409, 409, 409, 49, 596, 115, 101, 541, 541, 541, 795,…
## $ Race1            <fct> White, White, White, Other, White, White, White, Whit…
## $ Race3            <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ Education        <fct> High School, High School, High School, NA, Some Colle…
## $ MaritalStatus    <fct> Married, Married, Married, NA, LivePartner, NA, NA, M…
## $ HHIncome         <fct> 25000-34999, 25000-34999, 25000-34999, 20000-24999, 3…
## $ HHIncomeMid      <int> 30000, 30000, 30000, 22500, 40000, 87500, 60000, 8750…
## $ Poverty          <dbl> 1.36, 1.36, 1.36, 1.07, 1.91, 1.84, 2.33, 5.00, 5.00,…
## $ HomeRooms        <int> 6, 6, 6, 9, 5, 6, 7, 6, 6, 6, 5, 10, 6, 10, 10, 4, 3,…
## $ HomeOwn          <fct> Own, Own, Own, Own, Rent, Rent, Own, Own, Own, Own, O…
## $ Work             <fct> NotWorking, NotWorking, NotWorking, NA, NotWorking, N…
## $ Weight           <dbl> 87.4, 87.4, 87.4, 17.0, 86.7, 29.8, 35.2, 75.7, 75.7,…
## $ Length           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ HeadCirc         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ Height           <dbl> 164.7, 164.7, 164.7, 105.4, 168.4, 133.1, 130.6, 166.…
## $ BMI              <dbl> 32.22, 32.22, 32.22, 15.30, 30.57, 16.82, 20.64, 27.2…
## $ BMICatUnder20yrs <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ BMI_WHO          <fct> 30.0_plus, 30.0_plus, 30.0_plus, 12.0_18.5, 30.0_plus…
## $ Pulse            <int> 70, 70, 70, NA, 86, 82, 72, 62, 62, 62, 60, 62, 76, 8…
## $ BPSysAve         <int> 113, 113, 113, NA, 112, 86, 107, 118, 118, 118, 111, …
## $ BPDiaAve         <int> 85, 85, 85, NA, 75, 47, 37, 64, 64, 64, 63, 74, 85, 6…
## $ BPSys1           <int> 114, 114, 114, NA, 118, 84, 114, 106, 106, 106, 124, …
## $ BPDia1           <int> 88, 88, 88, NA, 82, 50, 46, 62, 62, 62, 64, 76, 86, 6…
## $ BPSys2           <int> 114, 114, 114, NA, 108, 84, 108, 118, 118, 118, 108, …
## $ BPDia2           <int> 88, 88, 88, NA, 74, 50, 36, 68, 68, 68, 62, 72, 88, 6…
## $ BPSys3           <int> 112, 112, 112, NA, 116, 88, 106, 118, 118, 118, 114, …
## $ BPDia3           <int> 82, 82, 82, NA, 76, 44, 38, 60, 60, 60, 64, 76, 82, 7…
## $ Testosterone     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ DirectChol       <dbl> 1.29, 1.29, 1.29, NA, 1.16, 1.34, 1.55, 2.12, 2.12, 2…
## $ TotChol          <dbl> 3.49, 3.49, 3.49, NA, 6.70, 4.86, 4.09, 5.82, 5.82, 5…
## $ UrineVol1        <int> 352, 352, 352, NA, 77, 123, 238, 106, 106, 106, 113, …
## $ UrineFlow1       <dbl> NA, NA, NA, NA, 0.094, 1.538, 1.322, 1.116, 1.116, 1.…
## $ UrineVol2        <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ UrineFlow2       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ Diabetes         <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N…
## $ DiabetesAge      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ HealthGen        <fct> Good, Good, Good, NA, Good, NA, NA, Vgood, Vgood, Vgo…
## $ DaysPhysHlthBad  <int> 0, 0, 0, NA, 0, NA, NA, 0, 0, 0, 10, 0, 4, NA, NA, 0,…
## $ DaysMentHlthBad  <int> 15, 15, 15, NA, 10, NA, NA, 3, 3, 3, 0, 0, 0, NA, NA,…
## $ LittleInterest   <fct> Most, Most, Most, NA, Several, NA, NA, None, None, No…
## $ Depressed        <fct> Several, Several, Several, NA, Several, NA, NA, None,…
## $ nPregnancies     <int> NA, NA, NA, NA, 2, NA, NA, 1, 1, 1, NA, NA, NA, NA, N…
## $ nBabies          <int> NA, NA, NA, NA, 2, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ Age1stBaby       <int> NA, NA, NA, NA, 27, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ SleepHrsNight    <int> 4, 4, 4, NA, 8, NA, NA, 8, 8, 8, 7, 5, 4, NA, 5, 7, N…
## $ SleepTrouble     <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No, No, Y…
## $ PhysActive       <fct> No, No, No, NA, No, NA, NA, Yes, Yes, Yes, Yes, Yes, …
## $ PhysActiveDays   <int> NA, NA, NA, NA, NA, NA, NA, 5, 5, 5, 7, 5, 1, NA, 2, …
## $ TVHrsDay         <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ CompHrsDay       <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TVHrsDayChild    <int> NA, NA, NA, 4, NA, 5, 1, NA, NA, NA, NA, NA, NA, 4, N…
## $ CompHrsDayChild  <int> NA, NA, NA, 1, NA, 0, 6, NA, NA, NA, NA, NA, NA, 3, N…
## $ Alcohol12PlusYr  <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, Yes, Y…
## $ AlcoholDay       <int> NA, NA, NA, NA, 2, NA, NA, 3, 3, 3, 1, 2, 6, NA, NA, …
## $ AlcoholYear      <int> 0, 0, 0, NA, 20, NA, NA, 52, 52, 52, 100, 104, 364, N…
## $ SmokeNow         <fct> No, No, No, NA, Yes, NA, NA, NA, NA, NA, No, NA, NA, …
## $ Smoke100         <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, Yes, No, …
## $ Smoke100n        <fct> Smoker, Smoker, Smoker, NA, Smoker, NA, NA, Non-Smoke…
## $ SmokeAge         <int> 18, 18, 18, NA, 38, NA, NA, NA, NA, NA, 13, NA, NA, N…
## $ Marijuana        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, NA, Ye…
## $ AgeFirstMarij    <int> 17, 17, 17, NA, 18, NA, NA, 13, 13, 13, NA, 19, 15, N…
## $ RegularMarij     <fct> No, No, No, NA, No, NA, NA, No, No, No, NA, Yes, Yes,…
## $ AgeRegMarij      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, 15, N…
## $ HardDrugs        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No, Yes, …
## $ SexEver          <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, Yes, Y…
## $ SexAge           <int> 16, 16, 16, NA, 12, NA, NA, 13, 13, 13, 17, 22, 12, N…
## $ SexNumPartnLife  <int> 8, 8, 8, NA, 10, NA, NA, 20, 20, 20, 15, 7, 100, NA, …
## $ SexNumPartYear   <int> 1, 1, 1, NA, 1, NA, NA, 0, 0, 0, NA, 1, 1, NA, NA, 1,…
## $ SameSex          <fct> No, No, No, NA, Yes, NA, NA, Yes, Yes, Yes, No, No, N…
## $ SexOrientation   <fct> Heterosexual, Heterosexual, Heterosexual, NA, Heteros…
## $ PregnantNow      <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

glimpse(survey)

## Rows: 237
## Columns: 12
## $ Sex    <fct> Female, Male, Male, Male, Male, Female, Male, Female, Male, Mal…
## $ Wr.Hnd <dbl> 18.5, 19.5, 18.0, 18.8, 20.0, 18.0, 17.7, 17.0, 20.0, 18.5, 17.…
## $ NW.Hnd <dbl> 18.0, 20.5, 13.3, 18.9, 20.0, 17.7, 17.7, 17.3, 19.5, 18.5, 17.…
## $ W.Hnd  <fct> Right, Left, Right, Right, Right, Right, Right, Right, Right, R…
## $ Fold   <fct> R on L, R on L, L on R, R on L, Neither, L on R, L on R, R on L…
## $ Pulse  <int> 92, 104, 87, NA, 35, 64, 83, 74, 72, 90, 80, 68, NA, 66, 60, NA…
## $ Clap   <fct> Left, Left, Neither, Neither, Right, Right, Right, Right, Right…
## $ Exer   <fct> Some, None, None, None, Some, Some, Freq, Freq, Some, Some, Fre…
## $ Smoke  <fct> Never, Regul, Occas, Never, Never, Never, Never, Never, Never, …
## $ Height <dbl> 173.00, 177.80, NA, 160.00, 165.00, 172.72, 182.88, 157.00, 175…
## $ M.I    <fct> Metric, Imperial, NA, Metric, Metric, Imperial, Imperial, Metri…
## $ Age    <dbl> 18.250, 17.583, 16.917, 20.333, 23.667, 21.000, 18.833, 35.833,…

A practical test-selection guide

Research question	Outcome	Comparison	Test
Does one sample mean differ from a reference value?	Continuous	One sample versus reference	One-sample t-test
Do two independent groups have different means?	Continuous	Two unrelated groups	Independent-samples t-test
Do two measurements from the same people differ?	Continuous	Paired measurements	Paired t-test
Do 3 or more independent groups have different means?	Continuous	3+ groups	One-way ANOVA
Are two categorical variables associated?	Categorical	Cross-tabulation	Chi-square test

Understanding H0, H1 and the p-value

For most tests:

• H0 represents no mean difference, no change, or no association.
• H1 represents a difference, change, or association.

We will use 0.05 as a conventional teaching threshold.

alpha <- 0.05

A p-value is calculated assuming H0. It is not the probability that H0 is true. A statistically significant result also does not automatically imply a large, important, clinically relevant or causal effect.

Part I: One-sample t-tests

A one-sample t-test compares the mean of one numerical variable with a specified reference value.

Example 1: BMI compared with 25 kg/m2

Research question

Among NHANES adults in this teaching subset, does mean BMI differ from 25 kg/m2?

H0: mean BMI = 25
H1: mean BMI != 25

Create a complete adult sample.

bmi_adults <- NHANES |>
  filter(Age >= 18, !is.na(BMI))

bmi_adults |>
  summarise(
    n = n(),
    mean = mean(BMI),
    sd = sd(BMI),
    median = median(BMI),
    q1 = quantile(BMI, .25),
    q3 = quantile(BMI, .75)
  )

## # A tibble: 1 × 6
##       n  mean    sd median    q1    q3
##   <int> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1  7414  28.7  6.70   27.6  24.0  32.1

ggplot(bmi_adults, aes(BMI)) +
  geom_histogram(bins = 30, colour = "white") +
  geom_vline(xintercept = 25, linetype = 2) +
  labs(
    title = "BMI among adults in the NHANES teaching dataset",
    x = "BMI (kg/m2)",
    y = "Count"
  )

Run the test.

bmi_t1 <- t.test(bmi_adults$BMI, mu = 25)
bmi_t1

## 
##  One Sample t-test
## 
## data:  bmi_adults$BMI
## t = 47.212, df = 7413, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 25
## 95 percent confidence interval:
##  28.52000 28.82497
## sample estimates:
## mean of x 
##  28.67249

Extract a compact result.

tibble(
  n = nrow(bmi_adults),
  mean = mean(bmi_adults$BMI),
  reference = 25,
  mean_difference = mean(bmi_adults$BMI) - 25,
  t = unname(bmi_t1$statistic),
  df = unname(bmi_t1$parameter),
  p_value = bmi_t1$p.value,
  ci_low = bmi_t1$conf.int[1],
  ci_high = bmi_t1$conf.int[2]
) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kable()

n	mean	reference	mean_difference	t	df	p_value	ci_low	ci_high
7414	28.672	25	3.672	47.212	7413	0	28.52	28.825

The test asks whether the observed sample mean is compatible with a population mean of 25 under the assumptions of the one-sample t-test. Interpretation should include the observed mean, its difference from the reference, the confidence interval and the p-value.

Example 2: Average sleep compared with 7 hours

A second health example uses reported hours of sleep.

sleep_adults <- NHANES |>
  filter(Age >= 18, !is.na(SleepHrsNight))

sleep_adults |>
  summarise(
    n = n(),
    mean_sleep = mean(SleepHrsNight),
    sd_sleep = sd(SleepHrsNight),
    median_sleep = median(SleepHrsNight)
  )

## # A tibble: 1 × 4
##       n mean_sleep sd_sleep median_sleep
##   <int>      <dbl>    <dbl>        <dbl>
## 1  7464       6.91     1.34            7

Research question: does mean reported nightly sleep differ from 7 hours?

sleep_t1 <- t.test(sleep_adults$SleepHrsNight, mu = 7)
sleep_t1

## 
##  One Sample t-test
## 
## data:  sleep_adults$SleepHrsNight
## t = -6.094, df = 7463, p-value = 1.156e-09
## alternative hypothesis: true mean is not equal to 7
## 95 percent confidence interval:
##  6.874809 6.935748
## sample estimates:
## mean of x 
##  6.905279

tibble(
  mean_sleep = mean(sleep_adults$SleepHrsNight),
  difference_from_7 = mean(sleep_adults$SleepHrsNight) - 7,
  t = unname(sleep_t1$statistic),
  df = unname(sleep_t1$parameter),
  p = sleep_t1$p.value,
  ci_low = sleep_t1$conf.int[1],
  ci_high = sleep_t1$conf.int[2]
) |>
  mutate(across(where(is.numeric), ~ round(.x, 3)))

## # A tibble: 1 × 7
##   mean_sleep difference_from_7     t    df     p ci_low ci_high
##        <dbl>             <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>
## 1       6.90            -0.095 -6.09  7463     0   6.88    6.94

Example 3: Systolic blood pressure compared with 120 mmHg

sbp_adults <- NHANES |>
  filter(Age >= 18, !is.na(BPSysAve))

sbp_adults |>
  summarise(
    n = n(),
    mean_sbp = mean(BPSysAve),
    sd_sbp = sd(BPSysAve),
    median_sbp = median(BPSysAve)
  )

## # A tibble: 1 × 4
##       n mean_sbp sd_sbp median_sbp
##   <int>    <dbl>  <dbl>      <int>
## 1  7205     121.   17.1        118

sbp_t1 <- t.test(sbp_adults$BPSysAve, mu = 120)
sbp_t1

## 
##  One Sample t-test
## 
## data:  sbp_adults$BPSysAve
## t = 3.6188, df = 7204, p-value = 0.0002981
## alternative hypothesis: true mean is not equal to 120
## 95 percent confidence interval:
##  120.3334 121.1215
## sample estimates:
## mean of x 
##  120.7274

These three examples use exactly the same statistical machinery. What changes is the substantive question and the reference value.

Assumptions for a one-sample t-test

The main object is the distribution of the numerical observations around their mean. Particularly with small samples, strong skewness or influential observations deserve attention.

ggplot(bmi_adults, aes(sample = BMI)) +
  stat_qq() +
  stat_qq_line() +
  labs(title = "Q-Q plot of adult BMI")

With large samples, a t procedure is often relatively robust to moderate non-normality, but large sample size does not make data-quality problems or dependence disappear.

Part II: Independent-samples t-tests

An independent-samples t-test compares the means of two unrelated groups. R’s default t.test() uses Welch’s method, which does not assume equal population variances.

Example 1: Systolic blood pressure by sex

Research question

Does mean systolic blood pressure differ between the two sex categories recorded in this NHANES dataset?

H0: the population means are equal.
H1: the population means differ.

sbp_sex <- NHANES |>
  filter(Age >= 18, !is.na(BPSysAve), !is.na(Gender))

sbp_sex |>
  summarise(
    n = n(),
    mean = mean(BPSysAve),
    sd = sd(BPSysAve),
    median = median(BPSysAve),
    .by = Gender
  ) |>
  kable(caption = "Systolic blood pressure by sex category")

Systolic blood pressure by sex category

Gender	n	mean	sd	median
male	3581	122.6898	15.89424	121
female	3624	118.7884	17.93533	116

ggplot(sbp_sex, aes(Gender, BPSysAve)) +
  geom_boxplot() +
  labs(
    title = "Systolic blood pressure by sex category",
    x = "Sex",
    y = "Average systolic blood pressure"
  )

sbp_sex_test <- t.test(BPSysAve ~ Gender, data = sbp_sex)
sbp_sex_test

## 
##  Welch Two Sample t-test
## 
## data:  BPSysAve by Gender
## t = -9.7746, df = 7119.4, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -4.683820 -3.118972
## sample estimates:
## mean in group female   mean in group male 
##             118.7884             122.6898

Example 2: BMI by current smoking status

First inspect the available smoking categories.

NHANES |>
  count(SmokeNow, .drop = FALSE)

## # A tibble: 3 × 2
##   SmokeNow     n
##   <fct>    <int>
## 1 No        1745
## 2 Yes       1466
## 3 <NA>      6789

Use complete adult observations.

bmi_smoke <- NHANES |>
  filter(Age >= 18, !is.na(BMI), !is.na(SmokeNow))

bmi_smoke |>
  summarise(
    n = n(),
    mean_bmi = mean(BMI),
    sd_bmi = sd(BMI),
    median_bmi = median(BMI),
    .by = SmokeNow
  )

## # A tibble: 2 × 5
##   SmokeNow     n mean_bmi sd_bmi median_bmi
##   <fct>    <int>    <dbl>  <dbl>      <dbl>
## 1 No        1726     29.1   6.19       28.3
## 2 Yes       1458     27.7   6.42       26.6

bmi_smoke_test <- t.test(BMI ~ SmokeNow, data = bmi_smoke)
bmi_smoke_test

## 
##  Welch Two Sample t-test
## 
## data:  BMI by SmokeNow
## t = 6.1667, df = 3052.9, p-value = 7.89e-10
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##  0.9444331 1.8249879
## sample estimates:
##  mean in group No mean in group Yes 
##          29.08616          27.70145

The test compares means. It does not show that smoking caused a BMI difference. Observational comparisons may reflect age, socioeconomic position, disease history, selection and other factors.

Example 3: Sleep hours by marital status, restricted to two groups

Here we compare two selected social groups to illustrate how an independent t-test can be applied to a social-science variable.

NHANES |>
  count(MaritalStatus, sort = TRUE)

## # A tibble: 7 × 2
##   MaritalStatus     n
##   <fct>         <int>
## 1 Married        3945
## 2 <NA>           2769
## 3 NeverMarried   1380
## 4 Divorced        707
## 5 LivePartner     560
## 6 Widowed         456
## 7 Separated       183

sleep_marital <- NHANES |>
  filter(
    Age >= 18,
    MaritalStatus %in% c("Married", "NeverMarried"),
    !is.na(SleepHrsNight)
  ) |>
  droplevels()

sleep_marital |>
  summarise(
    n = n(),
    mean_sleep = mean(SleepHrsNight),
    sd_sleep = sd(SleepHrsNight),
    .by = MaritalStatus
  )

## # A tibble: 2 × 4
##   MaritalStatus     n mean_sleep sd_sleep
##   <fct>         <int>      <dbl>    <dbl>
## 1 Married        3938       6.94     1.27
## 2 NeverMarried   1372       6.88     1.37

sleep_marital_test <- t.test(
  SleepHrsNight ~ MaritalStatus,
  data = sleep_marital
)

sleep_marital_test

## 
##  Welch Two Sample t-test
## 
## data:  SleepHrsNight by MaritalStatus
## t = 1.3122, df = 2252.1, p-value = 0.1896
## alternative hypothesis: true difference in means between group Married and group NeverMarried is not equal to 0
## 95 percent confidence interval:
##  -0.02729019  0.13768721
## sample estimates:
##      mean in group Married mean in group NeverMarried 
##                   6.939309                   6.884111

Again, this is an unadjusted descriptive comparison, not evidence that marital status causes differences in sleep.

Welch versus pooled t-test

Welch’s test is R’s default.

t.test(BPSysAve ~ Gender, data = sbp_sex, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  BPSysAve by Gender
## t = -9.7746, df = 7119.4, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -4.683820 -3.118972
## sample estimates:
## mean in group female   mean in group male 
##             118.7884             122.6898

The pooled version assumes equal population variances.

t.test(BPSysAve ~ Gender, data = sbp_sex, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  BPSysAve by Gender
## t = -9.7676, df = 7203, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -4.684381 -3.118412
## sample estimates:
## mean in group female   mean in group male 
##             118.7884             122.6898

Welch’s test is normally a strong default for independent groups because it does not require the equal-variance assumption.

Effect size for two independent means

A standardised mean difference can complement the p-value. If effectsize is installed:

if (requireNamespace("effectsize", quietly = TRUE)) {
  effectsize::cohens_d(BPSysAve ~ Gender, data = sbp_sex)
}

## Cohen's d |         95% CI
## --------------------------
## -0.23     | [-0.28, -0.18]
## 
## - Estimated using pooled SD.

Do not allow an effect-size label such as “small” or “large” to replace subject-matter interpretation.

Part III: Paired t-tests

Paired t-tests require genuinely linked measurements. NHANES contains repeated blood-pressure readings on the same examination participant, making it useful for demonstrating paired analysis.

Example 1: First versus second systolic BP reading

bp_pair12 <- NHANES |>
  filter(Age >= 18, !is.na(BPSys1), !is.na(BPSys2)) |>
  transmute(
    first = BPSys1,
    second = BPSys2,
    difference = second - first
  )

bp_pair12 |>
  summarise(
    pairs = n(),
    mean_first = mean(first),
    mean_second = mean(second),
    mean_difference = mean(difference),
    sd_difference = sd(difference)
  )

## # A tibble: 1 × 5
##   pairs mean_first mean_second mean_difference sd_difference
##   <int>      <dbl>       <dbl>           <dbl>         <dbl>
## 1  6807       122.        121.          -0.835          5.74

bp_pair12_test <- t.test(
  bp_pair12$second,
  bp_pair12$first,
  paired = TRUE
)

bp_pair12_test

## 
##  Paired t-test
## 
## data:  bp_pair12$second and bp_pair12$first
## t = -11.994, df = 6806, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.9715017 -0.6985439
## sample estimates:
## mean difference 
##      -0.8350228

The analysis concerns the mean within-person difference.

ggplot(bp_pair12, aes(difference)) +
  geom_histogram(bins = 30, colour = "white") +
  geom_vline(xintercept = 0, linetype = 2) +
  labs(
    title = "Difference between second and first systolic BP readings",
    x = "Second minus first reading",
    y = "Count"
  )

Example 2: Second versus third systolic BP reading

bp_pair23 <- NHANES |>
  filter(Age >= 18, !is.na(BPSys2), !is.na(BPSys3)) |>
  transmute(
    second = BPSys2,
    third = BPSys3,
    difference = third - second
  )

bp_pair23 |>
  summarise(
    pairs = n(),
    mean_second = mean(second),
    mean_third = mean(third),
    mean_difference = mean(difference),
    sd_difference = sd(difference)
  )

## # A tibble: 1 × 5
##   pairs mean_second mean_third mean_difference sd_difference
##   <int>       <dbl>      <dbl>           <dbl>         <dbl>
## 1  6956        121.       120.          -0.656          5.40

bp_pair23_test <- t.test(
  bp_pair23$third,
  bp_pair23$second,
  paired = TRUE
)

bp_pair23_test

## 
##  Paired t-test
## 
## data:  bp_pair23$third and bp_pair23$second
## t = -10.144, df = 6955, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.7832660 -0.5295574
## sample estimates:
## mean difference 
##      -0.6564117

Example 3: First versus third systolic BP reading

bp_pair13 <- NHANES |>
  filter(Age >= 18, !is.na(BPSys1), !is.na(BPSys3)) |>
  transmute(
    first = BPSys1,
    third = BPSys3,
    difference = third - first
  )

bp_pair13_test <- t.test(
  bp_pair13$third,
  bp_pair13$first,
  paired = TRUE
)

bp_pair13_test

## 
##  Paired t-test
## 
## data:  bp_pair13$third and bp_pair13$first
## t = -18.762, df = 6832, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -1.599906 -1.297211
## sample estimates:
## mean difference 
##       -1.448558

Why the paired design matters

For the first-versus-second comparison, compare the correct paired analysis with an independent analysis purely as a teaching demonstration.

t.test(bp_pair12$second, bp_pair12$first, paired = TRUE)

## 
##  Paired t-test
## 
## data:  bp_pair12$second and bp_pair12$first
## t = -11.994, df = 6806, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.9715017 -0.6985439
## sample estimates:
## mean difference 
##      -0.8350228

t.test(bp_pair12$second, bp_pair12$first, paired = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  bp_pair12$second and bp_pair12$first
## t = -2.8465, df = 13609, p-value = 0.004426
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.4100233 -0.2600223
## sample estimates:
## mean of x mean of y 
##  120.8947  121.7297

The second analysis ignores the fact that both readings belong to the same person. Equal vector lengths do not make an independent analysis appropriate; the study design determines dependence.

Assumptions for paired t-tests

The distribution that matters most is the distribution of the paired differences.

ggplot(bp_pair12, aes(sample = difference)) +
  stat_qq() +
  stat_qq_line() +
  labs(title = "Q-Q plot of within-person BP differences")

Part IV: One-way ANOVA

One-way ANOVA compares the means of a numerical outcome across three or more independent groups.

Example 1: BMI across education categories

Research question

Does mean BMI differ across education categories among adults with complete data?

bmi_education <- NHANES |>
  filter(Age >= 20, !is.na(BMI), !is.na(Education)) |>
  droplevels()

bmi_education |>
  summarise(
    n = n(),
    mean_bmi = mean(BMI),
    sd_bmi = sd(BMI),
    median_bmi = median(BMI),
    .by = Education
  ) |>
  arrange(mean_bmi) |>
  kable(caption = "BMI by educational attainment")

BMI by educational attainment

Education	n	mean_bmi	sd_bmi	median_bmi
College Grad	2086	27.59868	5.734922	26.500
9 - 11th Grade	880	29.13134	7.177855	27.955
Some College	2247	29.16932	6.891499	28.400
High School	1500	29.40979	7.098449	28.270
8th Grade	445	29.48836	6.105724	28.570

ggplot(bmi_education, aes(Education, BMI)) +
  geom_boxplot() +
  labs(
    title = "BMI across education categories",
    x = "Education",
    y = "BMI (kg/m2)"
  ) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Fit the ANOVA.

bmi_education_aov <- aov(BMI ~ Education, data = bmi_education)
summary(bmi_education_aov)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Education      4   4178  1044.6    23.9 <2e-16 ***
## Residuals   7153 312614    43.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The overall F-test evaluates H0 that all population group means are equal.

Post-hoc comparisons

TukeyHSD(bmi_education_aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = BMI ~ Education, data = bmi_education)
## 
## $Education
##                                    diff        lwr        upr     p adj
## 9 - 11th Grade-8th Grade    -0.35701864 -1.4062379  0.6922006 0.8859677
## High School-8th Grade       -0.07857288 -1.0522474  0.8951016 0.9994793
## Some College-8th Grade      -0.31903601 -1.2549493  0.6168773 0.8853018
## College Grad-8th Grade      -1.88967786 -2.8315426 -0.9478131 0.0000005
## High School-9 - 11th Grade   0.27844576 -0.4874703  1.0443618 0.8590992
## Some College-9 - 11th Grade  0.03798263 -0.6793172  0.7552824 0.9999022
## College Grad-9 - 11th Grade -1.53265922 -2.2577072 -0.8076113 0.0000001
## Some College-High School    -0.24046312 -0.8418782  0.3609520 0.8113330
## College Grad-High School    -1.81110498 -2.4217404 -1.2004695 0.0000000
## College Grad-Some College   -1.57064186 -2.1190644 -1.0222193 0.0000000

A significant overall test does not mean that every pair of education categories differs.

Example 2: Sleep hours across education categories

sleep_education <- NHANES |>
  filter(Age >= 20, !is.na(SleepHrsNight), !is.na(Education)) |>
  droplevels()

sleep_education |>
  summarise(
    n = n(),
    mean_sleep = mean(SleepHrsNight),
    sd_sleep = sd(SleepHrsNight),
    .by = Education
  )

## # A tibble: 5 × 4
##   Education          n mean_sleep sd_sleep
##   <fct>          <int>      <dbl>    <dbl>
## 1 High School     1515       6.78     1.41
## 2 Some College    2261       6.87     1.35
## 3 College Grad    2093       7.04     1.13
## 4 9 - 11th Grade   885       6.77     1.50
## 5 8th Grade        450       6.89     1.52

sleep_education_aov <- aov(
  SleepHrsNight ~ Education,
  data = sleep_education
)

summary(sleep_education_aov)

##               Df Sum Sq Mean Sq F value  Pr(>F)    
## Education      4     78  19.507   10.99 6.9e-09 ***
## Residuals   7199  12774   1.774                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(sleep_education_aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = SleepHrsNight ~ Education, data = sleep_education)
## 
## $Education
##                                     diff         lwr        upr     p adj
## 9 - 11th Grade-8th Grade    -0.119322034 -0.32975180 0.09110773 0.5316741
## High School-8th Grade       -0.115775578 -0.31090068 0.07934952 0.4852080
## Some College-8th Grade      -0.024248857 -0.21185742 0.16335971 0.9966915
## College Grad-8th Grade       0.144889314 -0.04396493 0.33374356 0.2229776
## High School-9 - 11th Grade   0.003546456 -0.15022372 0.15731663 0.9999965
## Some College-9 - 11th Grade  0.095073176 -0.04903940 0.23918575 0.3735587
## College Grad-9 - 11th Grade  0.264211347  0.11848082 0.40994187 0.0000076
## Some College-High School     0.091526720 -0.02914454 0.21219798 0.2334938
## College Grad-High School     0.260664891  0.13806594 0.38326385 0.0000001
## College Grad-Some College    0.169138171  0.05889448 0.27938186 0.0002777

Example 3: Systolic blood pressure across race/ethnicity categories

sbp_race <- NHANES |>
  filter(Age >= 18, !is.na(BPSysAve), !is.na(Race1)) |>
  droplevels()

sbp_race |>
  summarise(
    n = n(),
    mean_sbp = mean(BPSysAve),
    sd_sbp = sd(BPSysAve),
    .by = Race1
  ) |>
  kable(caption = "Systolic blood pressure across race/ethnicity categories")

Systolic blood pressure across race/ethnicity categories

Race1	n	mean_sbp	sd_sbp
White	4832	121.2132	16.90188
Mexican	619	118.0210	15.78582
Hispanic	406	118.4335	17.65959
Other	545	117.1615	16.59087
Black	803	123.4707	18.22333

sbp_race_aov <- aov(BPSysAve ~ Race1, data = sbp_race)
summary(sbp_race_aov)

##               Df  Sum Sq Mean Sq F value   Pr(>F)    
## Race1          4   20784    5196   18.02 9.84e-15 ***
## Residuals   7200 2076469     288                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is a descriptive unadjusted comparison. Race/ethnicity is a social classification, and differences in health outcomes must not be interpreted as intrinsic biological effects of group membership. Relevant structural, social, environmental and clinical factors require substantive analysis.

Example 4: Social-science example using MASS::survey

The MASS::survey dataset contains responses from university students.

Inspect its variables.

summary(survey)

##      Sex          Wr.Hnd          NW.Hnd        W.Hnd          Fold    
##  Female:118   Min.   :13.00   Min.   :12.50   Left : 18   L on R : 99  
##  Male  :118   1st Qu.:17.50   1st Qu.:17.50   Right:218   Neither: 18  
##  NA's  :  1   Median :18.50   Median :18.50   NA's :  1   R on L :120  
##               Mean   :18.67   Mean   :18.58                            
##               3rd Qu.:19.80   3rd Qu.:19.73                            
##               Max.   :23.20   Max.   :23.50                            
##               NA's   :1       NA's   :1                                
##      Pulse             Clap       Exer       Smoke         Height     
##  Min.   : 35.00   Left   : 39   Freq:115   Heavy: 11   Min.   :150.0  
##  1st Qu.: 66.00   Neither: 50   None: 24   Never:189   1st Qu.:165.0  
##  Median : 72.50   Right  :147   Some: 98   Occas: 19   Median :171.0  
##  Mean   : 74.15   NA's   :  1              Regul: 17   Mean   :172.4  
##  3rd Qu.: 80.00                            NA's :  1   3rd Qu.:180.0  
##  Max.   :104.00                                        Max.   :200.0  
##  NA's   :45                                            NA's   :28     
##        M.I           Age       
##  Imperial: 68   Min.   :16.75  
##  Metric  :141   1st Qu.:17.67  
##  NA's    : 28   Median :18.58  
##                 Mean   :20.37  
##                 3rd Qu.:20.17  
##                 Max.   :73.00  
## 

Compare student height across exercise categories.

height_exercise <- survey |>
  filter(!is.na(Height), !is.na(Exer)) |>
  droplevels()

height_exercise |>
  summarise(
    n = n(),
    mean_height = mean(Height),
    sd_height = sd(Height),
    .by = Exer
  )

##   Exer   n mean_height sd_height
## 1 Some  84    170.3969  9.471735
## 2 None  20    169.0280  9.470654
## 3 Freq 105    174.6067  9.778547

height_exercise_aov <- aov(Height ~ Exer, data = height_exercise)
summary(height_exercise_aov)

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## Exer          2   1076   537.8   5.802 0.00354 **
## Residuals   206  19095    92.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(height_exercise_aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Height ~ Exer, data = height_exercise)
## 
## $Exer
##                diff        lwr         upr     p adj
## None-Freq -5.578667 -11.123977 -0.03335612 0.0482672
## Some-Freq -4.209762  -7.536948 -0.88257557 0.0088321
## Some-None  1.368905  -4.286225  7.02403410 0.8354672

The example demonstrates ANOVA mechanics. A difference in student height across self-reported exercise categories should not be interpreted causally.

ANOVA diagnostics

Inspect residuals rather than treating ANOVA as a single command.

par(mfrow = c(2, 2))
plot(bmi_education_aov)

par(mfrow = c(1, 1))

Compare within-group sample sizes, means, SDs and variances.

bmi_education |>
  summarise(
    n = n(),
    mean = mean(BMI),
    sd = sd(BMI),
    variance = var(BMI),
    .by = Education
  )

## # A tibble: 5 × 5
##   Education          n  mean    sd variance
##   <fct>          <int> <dbl> <dbl>    <dbl>
## 1 High School     1500  29.4  7.10     50.4
## 2 Some College    2247  29.2  6.89     47.5
## 3 College Grad    2086  27.6  5.73     32.9
## 4 9 - 11th Grade   880  29.1  7.18     51.5
## 5 8th Grade        445  29.5  6.11     37.3

ANOVA effect size

Eta-squared can be calculated from the ANOVA sums of squares.

eta_squared_aov <- function(model) {
  tab <- summary(model)[[1]]
  ss_effect <- tab[1, "Sum Sq"]
  ss_total <- sum(tab[, "Sum Sq"])
  ss_effect / ss_total
}

eta_squared_aov(bmi_education_aov)

## [1] 0.01318979

eta_squared_aov(sleep_education_aov)

## [1] 0.006071288

Effect size adds magnitude information that the p-value does not provide.

Part V: Chi-square tests

The Pearson chi-square test of independence examines whether two categorical variables are associated.

Example 1: Current smoking and sex

smoke_gender <- NHANES |>
  filter(Age >= 18, !is.na(SmokeNow), !is.na(Gender)) |>
  with(table(Gender, SmokeNow))

smoke_gender

##         SmokeNow
## Gender    No Yes
##   female 764 638
##   male   981 828

addmargins(smoke_gender)

##         SmokeNow
## Gender     No  Yes  Sum
##   female  764  638 1402
##   male    981  828 1809
##   Sum    1745 1466 3211

Row percentages help describe the association before testing.

round(prop.table(smoke_gender, margin = 1) * 100, 1)

##         SmokeNow
## Gender     No  Yes
##   female 54.5 45.5
##   male   54.2 45.8

chi_smoke_gender <- chisq.test(smoke_gender)
chi_smoke_gender

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  smoke_gender
## X-squared = 0.012915, df = 1, p-value = 0.9095

Inspect expected counts.

round(chi_smoke_gender$expected, 2)

##         SmokeNow
## Gender       No    Yes
##   female 761.91 640.09
##   male   983.09 825.91

Inspect Pearson residuals.

round(chi_smoke_gender$residuals, 2)

##         SmokeNow
## Gender      No   Yes
##   female  0.08 -0.08
##   male   -0.07  0.07

Example 2: Diabetes status and sex

diabetes_gender <- NHANES |>
  filter(Age >= 18, !is.na(Diabetes), !is.na(Gender)) |>
  with(table(Gender, Diabetes))

diabetes_gender

##         Diabetes
## Gender     No  Yes
##   female 3449  344
##   male   3287  399

round(prop.table(diabetes_gender, margin = 1) * 100, 1)

##         Diabetes
## Gender     No  Yes
##   female 90.9  9.1
##   male   89.2 10.8

chi_diabetes_gender <- chisq.test(diabetes_gender)
chi_diabetes_gender

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  diabetes_gender
## X-squared = 6.2432, df = 1, p-value = 0.01247

chi_diabetes_gender$expected

##         Diabetes
## Gender         No      Yes
##   female 3416.185 376.8149
##   male   3319.815 366.1851

Example 3: Current smoking and education

smoke_education <- NHANES |>
  filter(Age >= 20, !is.na(SmokeNow), !is.na(Education)) |>
  with(table(Education, SmokeNow))

smoke_education

##                 SmokeNow
## Education         No Yes
##   8th Grade      122  95
##   9 - 11th Grade 198 328
##   High School    367 388
##   Some College   552 476
##   College Grad   506 178

Column or row percentages answer different questions. Here we use education-specific row percentages.

round(prop.table(smoke_education, margin = 1) * 100, 1)

##                 SmokeNow
## Education          No  Yes
##   8th Grade      56.2 43.8
##   9 - 11th Grade 37.6 62.4
##   High School    48.6 51.4
##   Some College   53.7 46.3
##   College Grad   74.0 26.0

chi_smoke_education <- chisq.test(smoke_education)
chi_smoke_education

## 
##  Pearson's Chi-squared test
## 
## data:  smoke_education
## X-squared = 175.89, df = 4, p-value < 2.2e-16

round(chi_smoke_education$expected, 2)

##                 SmokeNow
## Education            No    Yes
##   8th Grade      117.96  99.04
##   9 - 11th Grade 285.94 240.06
##   High School    410.43 344.57
##   Some College   558.83 469.17
##   College Grad   371.83 312.17

Example 4: Health insurance and education

insurance_education <- NHANES |>
  filter(Age >= 20, !is.na(HealthGen), !is.na(Education))

# Inspect categories before constructing a social/health cross-tabulation.
insurance_education |>
  count(Education, HealthGen)

## # A tibble: 25 × 3
##    Education      HealthGen     n
##    <fct>          <fct>     <int>
##  1 8th Grade      Excellent    24
##  2 8th Grade      Vgood        41
##  3 8th Grade      Good        129
##  4 8th Grade      Fair        141
##  5 8th Grade      Poor         34
##  6 9 - 11th Grade Excellent    45
##  7 9 - 11th Grade Vgood       162
##  8 9 - 11th Grade Good        337
##  9 9 - 11th Grade Fair        199
## 10 9 - 11th Grade Poor         37
## # ℹ 15 more rows

Use general self-rated health as the health variable.

health_education_tab <- with(
  insurance_education,
  table(Education, HealthGen)
)

health_education_tab

##                 HealthGen
## Education        Excellent Vgood Good Fair Poor
##   8th Grade             24    41  129  141   34
##   9 - 11th Grade        45   162  337  199   37
##   High School          113   372  603  247   36
##   Some College         196   684  885  229   54
##   College Grad         360   825  588  105   19

chi_health_education <- chisq.test(health_education_tab)
chi_health_education

## 
##  Pearson's Chi-squared test
## 
## data:  health_education_tab
## X-squared = 799.82, df = 16, p-value < 2.2e-16

round(chi_health_education$expected, 2)

##                 HealthGen
## Education        Excellent  Vgood   Good   Fair  Poor
##   8th Grade          42.12 118.95 145.09  52.57 10.27
##   9 - 11th Grade     89.04 251.43 306.69 111.12 21.72
##   High School       156.50 441.94 539.07 195.31 38.17
##   Some College      233.79 660.18 805.26 291.76 57.02
##   College Grad      216.55 611.50 745.89 270.25 52.82

Example 5: Social-science chi-square using MASS::survey

Examine smoking category and exercise frequency among students.

student_smoke_exer <- survey |>
  filter(!is.na(Smoke), !is.na(Exer)) |>
  with(table(Smoke, Exer))

student_smoke_exer

##        Exer
## Smoke   Freq None Some
##   Heavy    7    1    3
##   Never   87   18   84
##   Occas   12    3    4
##   Regul    9    1    7

round(prop.table(student_smoke_exer, margin = 1) * 100, 1)

##        Exer
## Smoke   Freq None Some
##   Heavy 63.6  9.1 27.3
##   Never 46.0  9.5 44.4
##   Occas 63.2 15.8 21.1
##   Regul 52.9  5.9 41.2

chi_student_smoke_exer <- chisq.test(student_smoke_exer)

## Warning in chisq.test(student_smoke_exer): Chi-squared approximation may be
## incorrect

chi_student_smoke_exer

## 
##  Pearson's Chi-squared test
## 
## data:  student_smoke_exer
## X-squared = 5.4885, df = 6, p-value = 0.4828

round(chi_student_smoke_exer$expected, 2)

##        Exer
## Smoke    Freq  None  Some
##   Heavy  5.36  1.07  4.57
##   Never 92.10 18.42 78.48
##   Occas  9.26  1.85  7.89
##   Regul  8.28  1.66  7.06

If expected counts are sparse, the ordinary chi-square approximation needs careful consideration.

Cramer’s V for chi-square effect size

cramers_v <- function(tab) {
  test <- suppressWarnings(chisq.test(tab))
  n <- sum(tab)
  r <- nrow(tab)
  c <- ncol(tab)
  sqrt(unname(test$statistic) / (n * min(r - 1, c - 1)))
}

cramers_v(smoke_gender)

## [1] 0.002005559

cramers_v(smoke_education)

## [1] 0.2340838

cramers_v(health_education_tab)

## [1] 0.1758663

Cramer’s V is a standardised measure of association. Its substantive importance should be interpreted in context.

Fisher’s exact test when counts are sparse

For a 2 x 2 table with small expected counts, Fisher’s exact test is commonly useful.

# General form:
# fisher.test(my_2_by_2_table)

For any chi-square example, inspect expected counts first.

min(chi_smoke_gender$expected)

## [1] 640.0909

min(chi_diabetes_gender$expected)

## [1] 366.1851

A reusable t-test reporting function

t_test_summary <- function(test_object) {
  tibble(
    statistic = unname(test_object$statistic),
    df = unname(test_object$parameter),
    p_value = test_object$p.value,
    ci_low = test_object$conf.int[1],
    ci_high = test_object$conf.int[2]
  )
}

t_test_summary(bmi_t1)

## # A tibble: 1 × 5
##   statistic    df p_value ci_low ci_high
##       <dbl> <dbl>   <dbl>  <dbl>   <dbl>
## 1      47.2  7413       0   28.5    28.8

t_test_summary(sbp_sex_test)

## # A tibble: 1 × 5
##   statistic    df  p_value ci_low ci_high
##       <dbl> <dbl>    <dbl>  <dbl>   <dbl>
## 1     -9.77 7119. 2.00e-22  -4.68   -3.12

t_test_summary(bp_pair12_test)

## # A tibble: 1 × 5
##   statistic    df  p_value ci_low ci_high
##       <dbl> <dbl>    <dbl>  <dbl>   <dbl>
## 1     -12.0  6806 8.18e-33 -0.972  -0.699

A reusable chi-square reporting function

chi_summary <- function(tab) {
  test <- suppressWarnings(chisq.test(tab))

  tibble(
    n = sum(tab),
    chi_square = unname(test$statistic),
    df = unname(test$parameter),
    p_value = test$p.value,
    minimum_expected = min(test$expected),
    cramers_v = cramers_v(tab)
  )
}

chi_summary(smoke_gender)

## # A tibble: 1 × 6
##       n chi_square    df p_value minimum_expected cramers_v
##   <int>      <dbl> <int>   <dbl>            <dbl>     <dbl>
## 1  3211     0.0129     1   0.910             640.   0.00201

chi_summary(smoke_education)

## # A tibble: 1 × 6
##       n chi_square    df  p_value minimum_expected cramers_v
##   <int>      <dbl> <int>    <dbl>            <dbl>     <dbl>
## 1  3210       176.     4 5.68e-37             99.0     0.234

How to write the result

One-sample t-test

A good report contains:

• the sample mean and SD;
• the reference mean;
• the mean difference;
• t statistic and degrees of freedom;
• p-value;
• confidence interval.

Avoid writing only “significant” or “not significant.”

Independent-samples t-test

Report the sample size, mean and SD for each group, followed by the Welch t statistic, degrees of freedom, p-value and confidence interval for the mean difference. Add an effect-size estimate when it helps the research question.

Paired t-test

Report the mean at each measurement, the mean within-person difference, SD of the differences, confidence interval, t statistic, degrees of freedom and p-value.

ANOVA

Report group descriptive statistics, F statistic, numerator and denominator degrees of freedom, and p-value. If post-hoc comparisons are performed, identify the correction or procedure used. An effect size such as eta-squared can help communicate magnitude.

Chi-square

Report the contingency-table counts and useful percentages, chi-square statistic, degrees of freedom and p-value. Check expected cell counts. Add Cramer’s V when a standardised association measure is useful.

Common mistakes in health and social-science applications

Test selection based only on variable names

A variable called group does not automatically imply an independent t-test. Ask how many groups exist and whether observations are independent.

Treating repeated blood-pressure measurements as independent

Multiple readings from the same participant are paired. The paired t-test uses within-person differences.

Converting a continuous health outcome into categories only to use chi-square

Categorising blood pressure, BMI or age may discard substantial information. Do it only when the categories have a clear analytical or reporting purpose.

Ignoring survey design in NHANES

The demonstrations in this chapter are deliberately simple. National population inference from NHANES requires the correct survey design, weights, strata, clusters and analytic subsample. A conventional unweighted t-test or ANOVA should not automatically be described as a nationally representative estimate.

Using p < 0.05 as the entire conclusion

A p-value does not tell us how large a difference is, whether the result matters for policy or health, or whether an observed association is causal.

Interpreting demographic associations as biological causation

Associations involving sex, race/ethnicity, education, marital status or other social categories require careful substantive interpretation. Social, environmental, structural and clinical pathways may be important.

Running many pairwise t-tests after a multi-group comparison

With three or more groups, an omnibus ANOVA followed by an appropriate multiple-comparison procedure is usually preferable to a large collection of unadjusted t-tests.

Ignoring expected counts in chi-square tables

Always inspect chisq.test(tab)$expected. Sparse cells can make the asymptotic approximation unreliable.

Quick practice exercises

Exercise 1: One-sample t-test

Using adults in NHANES, test whether mean Pulse differs from a reference value chosen before examining the sample result. Report n, mean, SD, t, df, p and the confidence interval.

Exercise 2: Independent t-test

Among adults with complete observations, compare Pulse between two categories of Gender. Start with descriptive statistics and a boxplot.

Exercise 3: Paired t-test

Compare BPDia1 with BPDia2 among participants who have both measurements. Create the paired difference first and examine it graphically.

Exercise 4: ANOVA

Compare BMI across categories of MaritalStatus, restricting the analysis to adults with complete observations. If the overall ANOVA is significant, use Tukey’s procedure and interpret selected pairwise contrasts.

Exercise 5: Chi-square

Create a table of Diabetes by Education among adults with complete data. Calculate row percentages, run the chi-square test, inspect expected counts and calculate Cramer’s V.

Practical reference

Situation	Main R command	Main quantity tested
One mean versus reference	t.test(x, mu = value)	Mean minus reference
Two independent means	t.test(y ~ group)	Difference in population means
Two linked measurements	t.test(x, y, paired = TRUE)	Mean paired difference
Three or more independent means	aov(y ~ group)	Equality of population means
Follow-up after one-way ANOVA	TukeyHSD(model)	Adjusted pairwise mean contrasts
Two categorical variables	chisq.test(table(x, y))	Independence/association
Sparse 2 x 2 table	fisher.test(table(x, y))	Exact categorical association

Final decision map

One continuous variable compared with a prespecified value
Use a one-sample t-test.

One continuous outcome and two unrelated groups
Use a Welch independent-samples t-test as a sensible default.

Two measurements from the same participant
Use a paired t-test and analyse within-person differences.

One continuous outcome and three or more independent groups
Use one-way ANOVA for the omnibus comparison, then a justified post-hoc method if specific pairwise comparisons are needed.

Two categorical variables
Use a chi-square test of independence when expected counts support the approximation. Consider Fisher’s exact test for an appropriate sparse 2 x 2 table.

Conclusion

Hypothesis testing is most useful when the statistical test follows the study design and research question. In health and social science, the same R function can answer very different substantive questions, so the analysis should never stop at obtaining a p-value.

The examples in this chapter use health and social-science datasets rather than fabricated records. NHANES provides useful examples involving BMI, blood pressure, sleep, smoking, education, marital status and self-rated health, while MASS::survey provides additional student social-science examples. Together they show why good hypothesis testing combines descriptive statistics, graphics, assumptions, effect size, uncertainty and careful interpretation.