Before we begin

Let’s use the SBP dataset again. We will run a few lines when each topic comes up in class. The order is: normal distribution, confidence intervals, then the t-distribution and small samples.

Open BTE3207_Advanced_Biostatistics.Rproj and run the chunks from the top. Keep the CSV in the dataset folder. The code uses base R; no extra analysis package is needed.

dataset_sbp <- read.csv("dataset/sbp_dataset_korea_2013-2014.csv")
head(dataset_sbp)

dim(dataset_sbp)

## [1] 1000000       7

table(dataset_sbp$SEX, useNA = "ifany")

## 
##      1      2 
## 510227 489773

table(dataset_sbp$DIS, useNA = "ifany")

## 
##      1      2      3      4 
##  53398 162826  43114 740662

colSums(is.na(dataset_sbp))

##   SEX BTH_G   SBP   DBP   FBS   DIS   BMI 
##     0     0     0     0     0     0     0

Changing 1 and 2 to 0 and 1

We want Female = 1, Male = 0. The variable name female tells us what 1 means. == asks whether two values are equal; <- saves a result.

dataset_sbp$female <- ifelse(dataset_sbp$SEX == 2, 1, 0)
table(dataset_sbp$SEX, dataset_sbp$female)

##    
##          0      1
##   1 510227      0
##   2      0 489773

ifelse(condition, yes, no) works like IF in Excel. For each row, R asks whether SEX is 2. If yes, put 1; otherwise, put 0. The table checks that 1 became 0 and 2 became 1.

mean(dataset_sbp$female)

## [1] 0.489773

mean(dataset_sbp$female) * 100

## [1] 48.9773

Voila! The mean of a 0/1 variable is the proportion coded 1. Here, it is the proportion of females: 48.98%. Reversing the coding would give the proportion of males instead.

Try this: run head(dataset_sbp$SEX - 1). Subtracting 1 also works here because the only original codes are 1 and 2. It would not be a general solution for other category codes. Do not overwrite SEX.

dataset_sbp$htn_history <- ifelse(dataset_sbp$DIS == 1 |
                                  dataset_sbp$DIS == 2, 1, 0)
dataset_sbp$dm_history <- ifelse(dataset_sbp$DIS == 1 |
                                 dataset_sbp$DIS == 3, 1, 0)
table(dataset_sbp$DIS, dataset_sbp$htn_history)

##    
##          0      1
##   1      0  53398
##   2      0 162826
##   3  43114      0
##   4 740662      0

table(dataset_sbp$DIS, dataset_sbp$dm_history)

##    
##          0      1
##   1      0  53398
##   2 162826      0
##   3      0  43114
##   4 740662      0

A small sample for class

For this exercise, treat the CSV as a teaching population and take a random sample. We are practicing sampling uncertainty; these calculations alone do not establish estimates for all Koreans.

set.seed(20260005)
class_sbp <- dataset_sbp[sample(nrow(dataset_sbp), 113), ]
x <- class_sbp$SBP
c(n = length(x), mean = mean(x), SD = sd(x))

##         n      mean        SD 
## 113.00000 122.76106  15.72037

set.seed() makes the random sample reproducible. Run it together with the sampling line when you want the same sample again.

The normal distribution: a short recap

We can use R as a calculator. pnorm() gives the probability to the left of a value under a normal distribution; qnorm() works in the opposite direction. Both functions are in the stats package, loaded with R.

pnorm(1) - pnorm(-1)

## [1] 0.6826895

pnorm(2) - pnorm(-2)

## [1] 0.9544997

qnorm(0.975)

## [1] 1.959964

The middle 95% of a normal distribution is about ±1.96 SD. ±2 SD gives about 95.45%, so “two SD” is a convenient approximation.

z_sbp <- (x - mean(x)) / sd(x)
c(mean = mean(z_sbp), SD = sd(z_sbp))

##          mean            SD 
## -1.385323e-16  1.000000e+00

hist(z_sbp, main = "Standardized SBP", xlab = "z-score")

sbp_01 <- (x - min(x)) / (max(x) - min(x))
range(sbp_01)

## [1] 0 1

mean(x <= 130)

## [1] 0.6637168

pnorm(130, mean = mean(x), sd = sd(x))

## [1] 0.6774147

hist(class_sbp$FBS, main = "FBS in our sample", xlab = "FBS (mg/dL)")

Confidence interval of a mean

SD describes variation between people. SE describes variation between sample means. For independent observations, we estimate SE by \(s / \sqrt{n}\).

n <- length(x)
x_bar <- mean(x)
se_mean <- sd(x) / sqrt(n)
c(mean = x_bar, SD = sd(x), SE = se_mean)

##       mean         SD         SE 
## 122.761062  15.720366   1.478848

An approximate 95% CI is sample mean ± 1.96 × SE. It is centered on the sample mean; the population mean is what we are trying to estimate.

ci_mean_normal <- x_bar + c(-1, 1) * qnorm(0.975) * se_mean
ci_mean_normal

## [1] 119.8626 125.6596

Our sample mean is 122.76 mmHg, and this approximate interval is 119.86 to 125.66 mmHg. Later in this lesson, we will use a t critical value because we estimated the population SD from the sample.

x_bar + c(-1, 1) * qnorm(0.95) * se_mean   # 90% CI

## [1] 120.3286 125.1935

x_bar + c(-1, 1) * qnorm(0.975) * se_mean  # 95% CI

## [1] 119.8626 125.6596

x_bar + c(-1, 1) * qnorm(0.995) * se_mean  # 99% CI

## [1] 118.9518 126.5703

fbs <- class_sbp$FBS
se_fbs <- sd(fbs) / sqrt(length(fbs))
mean(fbs) + c(-1, 1) * qnorm(0.975) * se_fbs

## [1]  94.46032 104.01756

Confidence interval of a proportion

Now use the 0/1 variable we made earlier. Let the outcome be recorded hypertension history, and estimate its proportion in our sample.

yes <- sum(class_sbp$htn_history)
n_people <- nrow(class_sbp)
p_hat <- yes / n_people
se_p <- sqrt(p_hat * (1 - p_hat) / n_people)
c(yes = yes, n = n_people, proportion = p_hat, SE = se_p)

##          yes            n   proportion           SE 
##  32.00000000 113.00000000   0.28318584   0.04238379

ci_p_normal <- p_hat + c(-1, 1) * qnorm(0.975) * se_p
ci_p_normal

## [1] 0.2001151 0.3662566

ci_p_normal * 100

## [1] 20.01151 36.62566

This is the normal approximation, often called a Wald interval. Its performance can be poor with small counts or a proportion near 0 or 1. We will compare it with an exact interval below.

Try this: use female instead of htn_history. What does the estimated proportion refer to now? Keep the original data unchanged.

What does 95% mean?

A CI estimates a population mean or proportion, not a range containing 95% of individual SBP values. Under the method’s assumptions, about 95% of intervals from repeated samples would contain the fixed true parameter. After we calculate one interval, it either contains that parameter or it does not.

set.seed(20260501)
ci_many <- replicate(100, {
  one_sample <- sample(dataset_sbp$SBP, 113, replace = TRUE)
  mean(one_sample) + c(-1, 1) * qnorm(0.975) *
    sd(one_sample) / sqrt(113)
})

pool_mean <- mean(dataset_sbp$SBP)
covered <- ci_many[1, ] <= pool_mean & pool_mean <= ci_many[2, ]
mean(covered)

## [1] 0.94

plot(NA, xlim = c(1, 100), ylim = range(ci_many),
     xlab = "Sample number", ylab = "95% CI for mean SBP (mmHg)")
segments(1:100, ci_many[1, ], 1:100, ci_many[2, ],
         col = ifelse(covered, "grey50", "red"))
abline(h = pool_mean, col = "blue", lwd = 2)

Comparing two groups: a preview

Let’s look at SBP by recorded sex. Give the categories readable labels.

class_sbp$sex <- factor(class_sbp$SEX, levels = c(1, 2),
                        labels = c("Male", "Female"))
boxplot(SBP ~ sex, data = class_sbp,
        ylab = "SBP (mmHg)", xlab = "Recorded sex")

Overlapping boxes do not answer whether population means differ. Next week, we will estimate and test the difference between means. We do not pair unrelated people by subtracting their rows one by one.

The t-distribution and degrees of freedom

When estimating one mean, the degrees of freedom are \(n - 1\). The t critical value accounts for uncertainty from estimating the population SD. It also applies to large samples, where it approaches the normal critical value.

qnorm(0.975)

## [1] 1.959964

qt(0.975, df = 4)

## [1] 2.776445

qt(0.975, df = n - 1)

## [1] 1.981372

curve(dnorm(x), from = -4, to = 4,
      ylab = "Density", xlab = "Value", lwd = 2)
curve(dt(x, df = 4), add = TRUE, col = "red", lwd = 2)
legend("topright", c("Standard normal", "t, df = 4"),
       col = c("black", "red"), lwd = 2)

The t-distribution has heavier tails. A small sample does not become normal just because we use a t method. The usual t interval is exact for independent normal observations and approximate more generally.

ci_mean_t <- x_bar + c(-1, 1) * qt(0.975, df = n - 1) * se_mean
ci_mean_t

## [1] 119.8309 125.6912

t.test(x)$conf.int

## [1] 119.8309 125.6912
## attr(,"conf.level")
## [1] 0.95

set.seed(20260502)
small_sbp <- sample(dataset_sbp$SBP, 12)
t.test(small_sbp)$conf.int

## [1] 109.2854 127.3812
## attr(,"conf.level")
## [1] 0.95

A small sample proportion

Use a small sample from the same dataset and count people with a history of diabetes. binom.test() takes a count of successes, not a proportion.

set.seed(20260503)
small_data <- dataset_sbp[sample(nrow(dataset_sbp), 20), ]
yes_small <- sum(small_data$dm_history)
c(yes = yes_small, n = nrow(small_data))

## yes   n 
##   0  20

binom.test(x = yes_small, n = nrow(small_data))$conf.int

## [1] 0.0000000 0.1684335
## attr(,"conf.level")
## [1] 0.95

Even when we observe 0 cases, the population proportion need not be 0. Check the upper endpoint before interpreting the result.

This is a Clopper–Pearson exact binomial interval. “Exact” refers to its binomial probability calculation; it does not fix biased sampling or dependent observations. Score intervals are another option, so proportion intervals do not all have to use the same method.

Try this: the lecture’s cold-symptom example has 3 successes out of 20. Run binom.test(3, 20)$conf.int. Compare it with 0.15 + c(-1, 1) * 1.96 * sqrt(0.15 * 0.85 / 20). The normal interval extends below 0, which cannot be a proportion.

Before next class

Can you explain what 1 means in female, htn_history, and dm_history? Can you identify whether an interval refers to a mean or a proportion? Next week, we will keep using the same dataset to compare groups.

Function details: t.test, binom.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.