Open the project and read the data

Open BTE3207_Advanced_Biostatistics.Rproj first. The CSV should be inside the project’s dataset folder. If R says “cannot open the connection”, check the file location before continuing.

# <- saves the result under a name.
dataset_sbp <- read.csv("dataset/sbp_dataset_korea_2013-2014.csv")
dim(dataset_sbp)

## [1] 1000000       7

names(dataset_sbp)

## [1] "SEX"   "BTH_G" "SBP"   "DBP"   "FBS"   "DIS"   "BMI"

read.csv() reads the file into a data frame (a table in R). dim() gives its numbers of rows and columns. names() lists the column names. This file has 1,000,000 rows, not 100 million.

DIS describes disease history, not a diagnosis made from today’s blood pressure. BTH_G codes are categories; we will not use them as actual ages.

Make today’s sample

We will use 1,000 observations so that the examples stay small. This file has no missing values in the variables used below.

set.seed(20260007)
row_number <- sample(nrow(dataset_sbp), 1000)
class_sbp <- dataset_sbp[row_number, ]

set.seed() makes the random selection reproducible. nrow() counts rows, and sample() selects 1,000 row numbers. Inside [rows, columns], a blank after the comma means keep all columns.

Everyone using this file and seed should get the same results. The sample differs from the one in the lecture screenshots, so the fitted numbers can differ. The interpretation is the same.

Change SEX from 1/2 to 0/1

The lecture uses Male = 0 and Female = 1. Keep the original SEX column and make a new column called female.

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

##    
##       0   1
##   1 502   0
##   2   0 498

Read the first line as: if SEX equals 2, write 1; otherwise, write 0.

• $ selects a column: class_sbp$SEX means the SEX column.
• == asks whether two values are equal. It does not assign a value.
• ifelse(condition, yes, no) chooses between two values.
• table() counts observations. Here it checks the conversion.

A value of female = 0 means Male, not missing. We changed the coding, not the observations.

Normal distribution and one mean

The four normal-distribution functions do different jobs. Let’s use a standard normal distribution (mean 0, SD 1).

dnorm(0)          # Height of the density curve at 0

## [1] 0.3989423

pnorm(1.96)       # Probability to the left of 1.96

## [1] 0.9750021

qnorm(0.975)      # Value with 97.5% of the area to its left

## [1] 1.959964

rnorm(5)         # Five random values; these are not the SBP data

## [1] -1.7807003  1.1036969  1.9532459 -0.1801066  0.1717042

dnorm() gives a density, not the probability of one exact value. pnorm() gives an area, and qnorm() goes from an area back to a value.

Now return to our SBP sample. As a classroom example, compare its population mean with 120 mmHg. This is a test reference, not a diagnostic cutoff.

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

## mean of x 
##   122.744

one_mean$conf.int

## [1] 121.817 123.671
## attr(,"conf.level")
## [1] 0.95

one_mean$p.value

## [1] 8.462256e-09

mu gives the mean under H0. The saved result contains several parts; $estimate, $conf.int, and $p.value select the mean estimate, 95% CI, and p-value. The p-value assumes H0 and the test assumptions; it is not the probability that H0 is true.

A proportion and a two-by-two table

For a 0/1 variable, the mean is the proportion of ones. Suppose H0 says that the female proportion is 0.5. We can calculate the normal test statistic directly.

p0 <- 0.5
p_hat <- mean(class_sbp$female)
se0 <- sqrt(p0 * (1 - p0) / nrow(class_sbp))
z <- (p_hat - p0) / se0
2 * pnorm(-abs(z))

## [1] 0.8993432

sqrt() takes a square root; abs() removes a sign. The last line is the two-sided p-value. se0 uses the null proportion and a square root: a variance is not a standard deviation. The normal approximation needs sufficiently large expected counts in both categories; here they are 500 and 500 under H0. This is the idea behind the slide’s z-test, without installing another package.

For the table example, make a 0/1 hypertension-history variable.

class_sbp$htn_history <- ifelse(class_sbp$DIS %in% c(1, 2), 1, 0)
history_table <- table(class_sbp$female, class_sbp$htn_history)
history_table

##    
##       0   1
##   0 413  89
##   1 376 122

c(1, 2) makes a vector. %in% means “is one of these values”. Rows are female codes (0, 1); columns are history codes (0 = No, 1 = Yes). For two independent groups, testing equal history proportions is also a test of independence in this two-by-two table.

chisq.test(history_table, correct = FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  history_table
## X-squared = 6.8804, df = 1, p-value = 0.008715

fisher.test(history_table)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  history_table
## p-value = 0.01044
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.095568 2.073220
## sample estimates:
## odds ratio 
##    1.50504

The chi-square test uses an approximation; inspect expected counts when they may be small. Fisher’s test is also available for a two-by-two table. Its odds ratio depends on the row and column order. We are reviewing the commands here; RR/OR is supplementary in this lecture.

Paired t-tests require related observations, such as two measurements from the same person. This CSV does not supply paired visits; do not pair unrelated rows. Use the before/after example from week 6 for that review.

Fit the model

Our question is: how does mean SBP differ between the two groups?

sex_model <- lm(SBP ~ female, data = class_sbp)
coef(sex_model)

## (Intercept)      female 
##  124.286853   -3.098098

lm() means linear model. The formula SBP ~ female reads as “describe SBP using female”. The outcome is on the left of ~; the predictor is on the right. data = class_sbp tells R where to find these columns. coef() extracts the estimated coefficients.

We save the model as sex_model so we can inspect it again without refitting it each time.

\[\widehat{SBP} = b_0 + b_1\,female\]

• (Intercept) is \(b_0\): the estimated mean when female = 0 (Male).
• female is \(b_1\): Female mean minus Male mean.
• When female = 1, the estimated mean is \(b_0+b_1\).

Let’s check that interpretation using the two sample means.

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

##     Male   Female 
## 124.2869 121.1888

mean(female_sbp) - mean(male_sbp)

## [1] -3.098098

Inside [...], the condition keeps only the matching observations. The estimated male mean is 124.29 mmHg, and the female-minus-male difference is -3.1 mmHg. Keep the minus sign. A negative coefficient here means lower mean SBP in the female group.

Is this the same as a t-test?

summary(sex_model)

## 
## Call:
## lm(formula = SBP ~ female, data = class_sbp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.189 -11.189  -1.189   8.811  68.811 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 124.2869     0.6635 187.325  < 2e-16 ***
## female       -3.0981     0.9402  -3.295  0.00102 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.87 on 998 degrees of freedom
## Multiple R-squared:  0.01076,    Adjusted R-squared:  0.009772 
## F-statistic: 10.86 on 1 and 998 DF,  p-value: 0.001018

confint(sex_model)

##                 2.5 %     97.5 %
## (Intercept) 122.98487 125.588833
## female       -4.94307  -1.253125

summary() prints the model results. For now, read the female row:

confint() gives the 95% confidence intervals by default. The row named female is the interval for our mean difference. We will return to residual standard error and R² later.

# Keep the same direction: Female minus Male.
t.test(female_sbp, male_sbp, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  female_sbp and male_sbp
## t = -3.2952, df = 998, p-value = 0.001018
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.943070 -1.253125
## sample estimates:
## mean of x mean of y 
##  121.1888  124.2869

Voila! For this binary predictor, ordinary lm() and the equal-variance two-sample t-test give the same difference, CI, and two-sided p-value. The default t.test() is Welch’s test, so its CI and p-value can differ. Using t.test(SBP ~ female, data = class_sbp) also reverses the comparison to group 0 minus group 1. Always check the group order.

This comparison assumes independent observations; the pooled t-test and usual model standard errors assume equal error variances. The groups can differ in age or other characteristics, so this unadjusted difference does not establish a causal effect.

Try this: calculate the fitted mean for female = 1 from the two coefficients. Does it match mean(female_sbp)?

Tell R that DIS contains categories

The numbers 1, 2, 3, and 4 are labels for disease-history groups. A change from 1 to 2 is not the same kind of measurement as a 1 mmHg increase. Do not fit SBP ~ DIS while DIS is numeric.

class_sbp$history <- factor(class_sbp$DIS, levels = 1:4)
levels(class_sbp$history)

## [1] "1" "2" "3" "4"

table(class_sbp$history)

## 
##   1   2   3   4 
##  56 155  46 743

factor() tells R that a column contains categories. 1:4 means c(1, 2, 3, 4). We kept numeric-looking labels so that they match the slides, but history is now a factor. The first level, 1 = both hypertension and diabetes histories, is the reference category.

history_model <- lm(SBP ~ history, data = class_sbp)
coef(history_model)

## (Intercept)    history2    history3    history4 
##  131.553571    1.710945   -7.597050  -11.743343

There are four groups, so R creates three comparisons with the reference. We do not need to create the three dummy columns by hand.

tapply(class_sbp$SBP, class_sbp$history, mean)

##        1        2        3        4 
## 131.5536 133.2645 123.9565 119.8102

coef(history_model)["(Intercept)"] +
  coef(history_model)["history4"]

## (Intercept) 
##    119.8102

tapply(values, groups, function) applies a function within each group. Here it calculates four means. The second command adds the reference mean and the group-4 difference, giving the group-4 mean. A name inside ["..."] selects that coefficient. The + at the end of the line tells R that the command continues on the next line.

Try this: find the mean SBP for DIS = 2 from the coefficients. The coefficient alone is a difference, not the group mean.

Fit a line and read its coefficients

Now the outcome is DBP, and the predictor is SBP. Keep the variables in this order to match the lecture.

bp_model <- lm(DBP ~ SBP, data = class_sbp)
coef(bp_model)

## (Intercept)         SBP 
##  15.8784649   0.4911648

\[\widehat{DBP} = b_0 + b_1\,SBP\]

The estimated slope is 0.491. People who differ by 1 mmHg in SBP differ by about that many mmHg in estimated mean DBP in this model. This association across people does not say that changing one person’s SBP causes their DBP to change.

The intercept is 15.88, the line’s fitted DBP at SBP = 0. That SBP value is outside a meaningful range for this example, so the intercept is needed to position the line but is not a useful clinical prediction.

Draw the same relationship

plot(DBP ~ SBP, data = class_sbp,
     xlab = "SBP (mmHg)", ylab = "DBP (mmHg)",
     pch = 16, col = "grey60")
abline(bp_model, col = "red", lwd = 2)

plot() draws a scatter plot: each point is one observation. With DBP ~ SBP, SBP is the horizontal axis and DBP is the vertical axis. xlab and ylab label the axes. pch = 16 selects filled circles, and col selects their color. abline() adds the fitted straight line; lwd = 2 makes it thicker. These options change appearance, not the fit.

Read the fitted mean at SBP = 120

new_sbp <- data.frame(SBP = c(120, 140))
predict(bp_model, newdata = new_sbp)

##        1        2 
## 74.81824 84.64154

data.frame() makes a small table. Its column must be named SBP, matching the predictor in our model. predict() evaluates the fitted line at the supplied values. The two results correspond to SBP 120 and 140, in that order. They are estimated mean DBPs, not the actual DBP of every person with that SBP.

How far is a point from the line?

A residual is observed DBP minus fitted DBP. fitted() gives the values on the line for the observations used to fit the model; residuals() gives their differences from the observed outcomes.

class_sbp$fitted_dbp <- fitted(bp_model)
class_sbp$residual_dbp <- residuals(bp_model)
head(class_sbp[c("DBP", "fitted_dbp", "residual_dbp")])

head() shows only the first six rows. Supplying column names inside class_sbp[...] selects those columns. Pick one row and check: DBP − fitted_dbp = residual_dbp. A positive residual is above the line; a negative residual is below it.

sse <- sum(class_sbp$residual_dbp^2)
sse

## [1] 39943.16

^2 squares each residual and sum() adds them. The least-squares line is the line that minimizes this sum of squared errors (SSE). Squaring prevents positive and negative residuals from cancelling.

Read the uncertainty in the slope

summary(bp_model)$coefficients

##               Estimate Std. Error   t value      Pr(>|t|)
## (Intercept) 15.8784649  1.6567216  9.584268  7.133565e-21
## SBP          0.4911648  0.0133986 36.657911 4.820774e-187

confint(bp_model)

##                  2.5 %     97.5 %
## (Intercept) 12.6274074 19.1295224
## SBP          0.4648722  0.5174575

Read the SBP row. Its estimate is the slope; the p-value tests H0: population slope = 0. The 95% CI gives a range of plausible population slopes under the model assumptions. confint() uses the appropriate t critical value, so we do not need to round it to 2 by hand.

For the usual regression CI and p-value, observations should be independent, the mean relationship should be linear, and the error variance should be constant. Normal errors give exact small-sample inference. For today’s practice we focus on interpreting the model, not on adding new diagnostic procedures.

Compare SBP 120 with SBP 140

Our direction is 120 minus 140, so the multiplier is −20.

difference <- (120 - 140) * coef(bp_model)["SBP"]
slope_ci <- confint(bp_model)["SBP", ]
difference_ci <- range((120 - 140) * slope_ci)
difference

##       SBP 
## -9.823296

difference_ci

## [1] -10.349150  -9.297443

For the CI matrix, ["SBP", ] selects the SBP row and both columns. Multiplying an interval by a negative number reverses its endpoints. range() returns the smaller and larger values in the correct order.

The estimated mean DBP difference is -9.82 mmHg, with a 95% CI from -10.35 to -9.3 mmHg. This is an interval for a difference in modelled means, not for two individual people’s DBPs.

Residual standard error

sqrt(sse / (nrow(class_sbp) - 2))

## [1] 6.32639

summary(bp_model)$sigma

## [1] 6.32639

These give the same residual standard error, in DBP’s units (mmHg). The denominator is n − 2 because this simple line estimates two coefficients: an intercept and a slope.

Calculate R²

First compare two ways to describe DBP: one overall mean for everyone, or a fitted mean that changes with SBP.

sst <- sum((class_sbp$DBP - mean(class_sbp$DBP))^2)
1 - sse / sst

## [1] 0.5738325

summary(bp_model)$r.squared

## [1] 0.5738325

sst is the total squared variation around the overall DBP mean. sse is the squared variation left around our fitted line. Thus, \(R^2 = 1 - SSE/SST\). The hand calculation and model output agree.

In this sample, the line accounts for about 57.4% of DBP variation. This is an in-sample description, not a percentage caused by SBP or a guarantee of accurate predictions for new people.

Compare the lecture’s other examples

bmi_model <- lm(BMI ~ SBP, data = class_sbp)
dbp_sex_model <- lm(DBP ~ female, data = class_sbp)

c(DBP_from_SBP = summary(bp_model)$r.squared,
  BMI_from_SBP = summary(bmi_model)$r.squared,
  DBP_from_female = summary(dbp_sex_model)$r.squared)

##    DBP_from_SBP    BMI_from_SBP DBP_from_female 
##      0.57383254      0.11058836      0.02303494

Read each as the fraction of variation explained in its own outcome. A small R² can coexist with a small p-value. The slope describes the size and direction of an association in the variables’ units; R² describes the fraction of variation explained. A larger slope does not automatically mean a stronger fit.

Pearson: a linear relationship

r <- cor(class_sbp$SBP, class_sbp$DBP)
r

## [1] 0.7575173

r^2

## [1] 0.5738325

summary(bp_model)$r.squared

## [1] 0.5738325

cor(x, y) computes Pearson’s correlation by default. It ranges from −1 to 1. Its sign gives the direction; its magnitude describes the strength of a linear relationship. Unlike the slope, it has no units.

Here, \(r^2=R^2\) because we fitted an ordinary simple linear regression with an intercept to the same observations, using one numeric predictor. Do not apply this shortcut to the four-level DIS model. For reference, \(\sqrt{0.583}\) is about 0.764, not 0.73.

Spearman: use the ranks

cor(class_sbp$SBP, class_sbp$DBP, method = "spearman")

## [1] 0.7360043

method = "spearman" changes the calculation to correlation of ranks. It describes whether larger values of one variable tend to go with larger (or smaller) values of the other in a monotonic relationship, which need not be a straight line. Ties receive average ranks. It does not measure every possible nonlinear relationship, and correlation alone does not establish causation.

Try this: compare the Pearson and Spearman results. Are their signs the same? They need not give identical numbers.