Lesson 5: Regression, Linearity & Dose-Response

EOH 635 — Environmental Epidemiology Methods

Author

Leah Schinasi

Published

October 14, 2026

In Lesson 4 you decided how to code your exposure. This lesson puts that exposure into a regression model: choosing a model based on your outcome, interpreting the estimate it produces, and checking whether the assumption that the exposure is related to the outcome “in a straight line” actually holds. It has two parts:

Choosing, fitting, and interpreting a crude model — starting from a two-by-two table you can calculate by hand, then making the same comparison with a logistic and a log-binomial regression model, then entering the exposure as a continuous term. Worked example: tree canopy cover and high blood pressure, the course’s primary exposure–outcome relationship.
Checking the linearity assumption — categorizing the exposure, fitting a model with one indicator term per category, and plotting the estimates to see whether the relationship really is a straight line.

Both parts end with a 🧠 Your Turn for the exposure and outcome in your own project. Together they produce the crude model and the evaluation of functional form that Milestone 2 asks for.

How to use this tutorial

The gray boxes below are live R editors — the code actually runs in your browser (no install needed). Click the green Run Code arrow in the top-left corner of a box to run it. The output will appear directly underneath.

You can also edit the code first — click inside the box, change something, then click Run Code to run your edited version. Nothing you do here can break anything; if you make a mess, just refresh the page to start that box over.

Sections marked 🧠 Your Turn give you a task and a blank or partially blank code editor. Try writing the code yourself before revealing the solution, hidden in a collapsed box underneath. Boxes labeled 💡 Answers work the same way: think through the questions first, then click to check your reasoning.

A note on browser vs. RStudio code: chunks labeled #| autorun: true fetch data from the web in the background so it’s available in your browser. When you write your own scripts in RStudio, you will not include this fetching code. Immediately after each one, a green box labeled “In your RStudio script” shows the simpler code you’ll actually use.

Doing this in RStudio instead of (or in addition to) the browser?

This tutorial uses class_data_processed.csv, the same merged file from Lessons 1, 3, and 4. If you saved your Lesson 4 exposure coding in your project script, run that code first so your categorical and binary exposure variables exist. See the Data Codebook for variable definitions.

Part 1: Choosing, Fitting, and Interpreting a Crude Model

In RStudio:

data <- read.csv(file.path("data_processed", "class_data_processed.csv"))

In the browser:

💻 In your RStudio script

data <- read.csv(file.path("data_processed", "class_data_processed.csv"))

A reminder: how is tree canopy coded?

Before modeling an exposure, remind yourself what its values mean and what its distribution looks like — you characterized this in Lesson 3, and it matters again here.

A few things to keep in mind about hrlctree:

It is a proportion, from 0 to 1: a value of 0.12 means 12% of the participant’s census tract is covered by tree canopy. In our data it ranges from about 0.02 (2%) to 0.63 (63%), with a median of about 0.12 (12%).
It is measured at the census-tract level, not for each person: every participant who lives in the same tract has the same value. The 2,000 participants live in 363 tracts.
It is right-skewed: most participants live in tracts with 5–20% canopy, and a small number live in much greener tracts. Keep this in mind in Part 2 — it affects how wide each quartile is.

Step 1: What kind of outcome do you have — and is it rare or common?

The type of outcome determines the model. high_blood_pressure is binary (1 = diagnosed, 0 = not), so a linear regression model is not appropriate — we need logistic or log-binomial regression. Which of those two to prefer depends on whether the outcome is rare or common.

Questions to consider:

What proportion of participants have high blood pressure? Would you call that rare or common?
How many participants are missing this outcome?

💡 Answers — click to reveal after you’ve thought about it

About 30% of participants with data (606 of 1,988) have high blood pressure. That is a common outcome — a rough rule of thumb is that an outcome is “rare” if it affects fewer than about 10% of people.
12 participants are missing high_blood_pressure. R will drop them automatically from every model below, so all models are fit on the 1,988 participants with data.

Step 2: A two-by-two table you can calculate by hand

Before fitting any model, let’s calculate the association the old-fashioned way, so you can see exactly what the regression model will reproduce. We’ll use tree canopy dichotomized at the median — the same kind of binary exposure you created in Lesson 4.

What each step of the code does:

Create a binary exposure. median(hrlctree, na.rm = TRUE) finds the median tree canopy value (about 0.12, or 12%); na.rm = TRUE tells R to ignore any missing values when doing so. ifelse() then checks each participant: if their tract’s canopy is above the median, the new variable tree_high is set to 1 (“high canopy”); otherwise it is set to 0 (“low canopy”). mutate() adds tree_high to the data as a new column, leaving hrlctree unchanged.
Build the two-by-two table. table() counts how many participants fall into each combination of tree_high (the rows) and high_blood_pressure (the columns). The names before the = signs (Tree_high, HBP) are just labels for the printed table. Typing two_by_two on its own line prints it. Participants missing high_blood_pressure are left out of the table.
Save the four counts. Each line stores one cell of the table under a letter — A, B, C, or D — so that we can use it in the formulas that follow. The numbers in square brackets tell R which cell to take: the first is the row and the second is the column. The table below the code shows which group of people each letter represents.
Calculate the prevalence in each group. Prevalence is the number of people with the outcome divided by everyone in that group. Among high-canopy participants that is A / (A + B); among low-canopy participants it is C / (C + D).
Calculate the measures of association. The prevalence ratio divides the prevalence in the high-canopy group by the prevalence in the low-canopy group. The odds ratio is the cross-product (A × D) / (B × C). c() combines the four results into one labeled list, so they print together.

Each of the four cells in the table counts a different group of people:

	High blood pressure (HBP = 1)	No high blood pressure (HBP = 0)
High canopy (`tree_high = 1`)	A = people in high-canopy tracts with high blood pressure	B = people in high-canopy tracts without high blood pressure
Low canopy (`tree_high = 0`)	C = people in low-canopy tracts with high blood pressure	D = people in low-canopy tracts without high blood pressure

Note that R prints the table with the 0 (low canopy, no high blood pressure) row and column first, so the printed table is arranged D, C / B, A — which is why the code picks each count out by position.

Prevalence, not risk

Because the class dataset is a cross-sectional survey — everyone was surveyed at one point in time — we are comparing the prevalence of high blood pressure (who has it now), not the risk of developing it over follow-up. That is why we call this a prevalence ratio rather than a risk ratio.

Questions to consider:

What is the prevalence of high blood pressure in each canopy group?
Is the prevalence ratio (PR) or the odds ratio (OR) further from the null value of 1? Why does that happen with an outcome this common?

💡 Answers — click to reveal after you’ve thought about it

The prevalence is 25.5% in high-canopy tracts (252 of 987) and 35.4% in low-canopy tracts (354 of 1,001).
PR = 0.255 / 0.354 = 0.72; OR = 0.63. The OR is further from the null. Odds are prevalence divided by one minus prevalence, p / (1 − p). When p is small, 1 − p is close to 1, so odds and prevalence are nearly the same. When p is large (here, around 0.3), dividing by 1 − p stretches the odds away from the prevalence, and the OR ends up further from 1 than the PR. Describing the OR of 0.63 as if it were a PR would overstate how much lower high blood pressure is in greener tracts.

Step 3: The same comparison as a regression model

Now, we will make the same comparison, running both a logistic and a log-binomial regression model. We will use the results from these models to calculate the odds ratio and its 95% confidence interval (CI), and the prevalence ratio and its 95% CI.

The two models differ only in their link function — what the model treats as a straight-line function of the exposure:

Logistic regression, binomial(link = "logit"), models the log odds of the outcome.
Log-binomial regression, binomial(link = "log"), models the log prevalence.

Why do we exponentiate the coefficient and the CI?

Because both models work on the log scale, the coefficient R reports for tree_high is not a ratio — it is a difference in logs. For the log-binomial model:

log(prevalence in high-canopy tracts) = β₀ + β₁
log(prevalence in low-canopy tracts) = β₀
the difference, β₁, is log(prevalence high) − log(prevalence low) = log(PR)

To get back from log(PR) to the PR itself, we undo the log by exponentiating: PR = exp(β₁). (The same logic applied to log odds gives OR = exp(β₁) for the logistic model.)

The 95% CI is built on the same log scale, as β₁ ± 1.96 × SE, so we exponentiate its two endpoints the same way: exp(β₁ − 1.96 × SE) and exp(β₁ + 1.96 × SE). That is also why ratio CIs are not symmetric around the estimate — they are symmetric on the log scale, not on the ratio scale.

The extract_ratio() helper below does these steps for us: it pulls out each coefficient and its standard error, then exponentiates the coefficient and both CI endpoints. It is the same idea as the extract_or() function from Lesson 4, just renamed because it now produces prevalence ratios as well as odds ratios.

Questions to consider:

Compare each estimate to the one you calculated by hand in Step 2. Do they match?
Write a one-sentence interpretation of the PR and its 95% CI.

💡 Answers — click to reveal after you’ve written your own

Yes — the logistic model’s OR (0.63) and the log-binomial model’s PR (0.72) match the two-by-two table exactly. A crude model with a single binary exposure is doing the same calculation as the table; what the model adds is a confidence interval, and (from next week) the ability to adjust for other variables.
Example interpretation: The prevalence of high blood pressure among participants living in census tracts with tree canopy above the median was 0.72 times the prevalence among those living in tracts at or below the median (95% CI: 0.63, 0.83) — about 28% lower. Because the whole CI is below 1, the data are not consistent with no association (a PR of 1).

Step 4: Entering the exposure as a continuous term

Dichotomizing at the median throws away information. By creating a binary variable in this way, we are in essence treating someone who lives in a census tract with 13% tree canopy coverage (just greater than the median) the same as someone living in a tract with ~63% coverage (the maximum of the distribution). Entering hrlctree as a continuous term uses all the information available. But recall that hrlctree is a proportion (0 to 1), so a “one-unit increase” would mean going from 0% to 100% canopy — a meaningless comparison. We rescale by multiplying by 100, so that the exposure is a percentage (0 to 100) and a one-unit increase means a one-percentage-point increase in tree canopy cover (for example, from 12% to 13%).

Questions to consider:

Interpret the PR per one-percentage-point increase in tree canopy.
This model assumes that every one-point increase — from 5% to 6%, or from 40% to 41% — multiplies prevalence by the same amount. Is that a reasonable assumption?

💡 Answers — click to reveal after you’ve thought about it

Each one-percentage-point increase in census-tract tree canopy cover was associated with a 2.7% lower prevalence of high blood pressure (PR = 0.973; 95% CI: 0.964, 0.983). Because the PR for a single percentage point is close to 1, report it to three decimal places so the estimate and CI are distinguishable.
A one-point change is small relative to the range of canopy (about 2% to 63%). Because the model multiplies the PR for each additional point, you can also describe a larger contrast by raising the PR to that power: for a 10-point increase, 0.973¹⁰ ≈ 0.76, or about 24% lower prevalence. (In R: exp(10 * coef(model_cont)[2]).) The model and its assumptions are the same; only the size of the comparison you report changes.
We can’t tell from this model — a single continuous term will return one PR no matter what the true shape of the relationship is. That is exactly what Part 2 checks.

💡 What if my log-binomial model won’t converge?

Log-binomial models sometimes fail with an error like “no valid set of coefficients has been found: please supply starting values”. This happens because predicted prevalences must stay below 1, which the fitting algorithm can struggle with — especially when an exposure has a few extremely high values. Two common fixes:

Supply starting values. Set the intercept to the log of the overall prevalence and the other coefficients to 0, e.g. start = c(log(0.3), 0) for a model with one exposure term.
Use a modified Poisson model — a Poisson model with robust (“sandwich”) standard errors, which also estimates a prevalence ratio.

Code: running a modified Poisson model (in RStudio)

# Packages for robust ("sandwich") standard errors -- install once with install.packages()
library(sandwich)
library(lmtest)

# Fit a Poisson model with a log link (same outcome and exposure as the log-binomial model)
model_pois <- glm(high_blood_pressure ~ hrlctree_pct, data = data,
                  family = poisson(link = "log"))

# Replace the usual standard errors with robust ones
robust_se <- vcovHC(model_pois, type = "HC0")

coeftest(model_pois, vcov = robust_se)      # coefficients with robust SEs
exp(coef(model_pois))                       # PR point estimate
exp(coefci(model_pois, vcov = robust_se))   # PR 95% CI

🧠 Your Turn: Fit a crude model for your own project

Using the exposure and outcome from your project (and the exposure coding you decided on in Lesson 4):

Check whether your outcome is binary or continuous and, if binary, how common it is.
Fit the appropriate crude model: log-binomial (or logistic, if you can justify it) for a binary outcome, or lm() for a continuous one. If your exposure is continuous and a one-unit change isn’t meaningful (e.g., it’s a proportion), rescale it first.
Use extract_ratio() to get your estimate and 95% CI (skip the exponentiating for a linear model — use coef() and confint() instead), and write a one-sentence interpretation.

💡 Example: traffic density (ctb1traff) and asthma — click to reveal

Interpretation: the prevalence of asthma among participants living in census tracts with traffic density above the median was 1.81 times the prevalence among those in tracts at or below the median (95% CI: 1.37, 2.38).

Asthma affects only about 10% of participants, so the OR (1.93) and PR are closer to each other here than they were for high blood pressure — the rarer the outcome, the better the OR approximates the PR.

Why dichotomize rather than enter traffic as a continuous term? Traffic density is very right-skewed — a handful of tracts have traffic many times higher than everyone else — and a log-binomial model with continuous traffic fails to converge for exactly the reason described in the “won’t converge” box above. The quartile approach in Part 2 is another way to handle it.

Part 2: Checking the Linearity Assumption

The continuous model in Step 4 produced one number to summarize the whole relationship. The check below asks whether that one number is a fair summary. The steps:

Categorize the exposure into several groups (quartiles here).
Fit a model with a separate indicator term for each category, compared with a reference category.
Plot the estimates. If they change by roughly the same amount from one category to the next, the linearity assumption is reasonable; if they bend, level off, or jump, it isn’t.

Step 1: Categorize tree canopy into quartiles

ntile(hrlctree, 4) sorts participants from lowest to highest tree canopy and splits them into four groups of (nearly) equal size — 500 participants each:

Q1 — the 25% of participants living in the least green tracts (about 2% to 8% canopy)
Q2 — the next 25% (about 8% to 12% canopy)
Q3 — the next 25% (about 12% to 18% canopy)
Q4 — the 25% of participants living in the greenest tracts (about 18% to 63% canopy)

Notice that the quartiles have equal numbers of people, not equal ranges of canopy: because canopy is right-skewed, Q4 spans a much wider range of canopy (about 45 percentage points) than Q1–Q3 combined (about 16).

Look at the last column: prevalence of high blood pressure is about 35–36% in Q1 and Q2, 31% in Q3, and 20% in Q4.

Step 2: Fit the model with quartiles as indicator terms

Wrapping the quartile variable in factor() is what tells R to create a separate indicator term for each category. Without it, R treats 1, 2, 3, 4 as a continuous number and returns a single coefficient — which would assume linearity instead of testing it.

What “ref = Q1” means. With a categorical exposure, R needs a reference category that every other category is compared with. By default it uses the first level — here, Q1, the least green quartile. Q1 therefore does not get its own row in the output: its PR is 1 by definition. Each PR you see is the prevalence of high blood pressure in that quartile divided by the prevalence in Q1. That is why every interpretation below says “compared with participants in the lowest quartile.”

Interpreting the results:

Quartile	PR (95% CI)	Interpretation
Q1 (2–8% canopy)	1 (reference)	The comparison group for every other row
Q2 (8–12%)	1.03 (0.87, 1.22)	Prevalence essentially the same as Q1 (3% higher). The CI includes 1, so the data are consistent with no difference.
Q3 (12–18%)	0.89 (0.74, 1.06)	Prevalence 11% lower than Q1, but the CI still includes 1, so this could plausibly be no difference.
Q4 (18–63%)	0.58 (0.47, 0.72)	Prevalence 42% lower than Q1. The whole CI is below 1: the data are not consistent with no difference.

In words, for Q4: the prevalence of high blood pressure among participants living in the greenest quarter of tracts was 0.58 times the prevalence among those in the least green quarter (95% CI: 0.47, 0.72).

Step 3: Plot the estimates

We plot each quartile’s PR, with Q1 fixed at 1, against the quartile itself. Each quartile gets the same amount of space on the x-axis, so you can see directly how much the PR changes from one quartile to the next. The canopy range of each quartile is shown under its label.

The y-axis is on a log scale because the log-binomial model is linear on the log scale: if the linearity assumption holds, the PRs should step down (or up) by about the same amount from each quartile to the next.

Questions to consider:

Do Q1, Q2, and Q3 look different from each other? Where does most of the change happen?
Does the line fall at a steady rate across the quartiles, or does it stay flat and then drop?
Now look at the canopy ranges under each label. How does the width of Q4 change your reading of the “drop”?
Given this plot, would you keep tree canopy as a single continuous term, use the quartile categories, or collapse some categories together?

💡 Discussion — click to reveal after you’ve looked at the plot yourself

Q2 is essentially the same as Q1, Q3 is slightly lower, and Q4 is much lower (PR ≈ 0.58). On this plot the line is roughly flat from Q1 to Q3 and then drops sharply to Q4 — not the steady, equal-sized steps we would expect if the linearity assumption held. That is a sign the assumption may not hold.

The canopy ranges add an important caveat. The quartiles are equal-sized groups of people, not equal-width ranges of canopy: Q1–Q3 together cover only about 2% to 18% canopy, while Q4 alone covers about 18% to 63%. Part of the “sudden drop” is simply that Q4 contains a much wider — and much greener — range of tracts. Measured in canopy percentage points, the decline above roughly 10% canopy is more gradual than this plot makes it look. This is why it is always worth checking what range of the exposure each category covers before interpreting the shape.

There is no single right answer about what to do next — reasonable options include:

Keep the continuous term, if you judge the departure from linearity to be small. It is the most efficient and easiest to communicate (one PR per percentage point of canopy).
Collapse categories that look alike. Here, Q1–Q3 could be combined into one reference group and compared with Q4 (the top quartile): fewer coefficients, and the estimate reflects where the change actually appears.
Use a flexible term such as a spline (optional section below).

Whatever you choose, justify it using the plot, your sample size in each category, and what will be most meaningful for your research question.

Collapsing categories that look alike

Question to consider: interpret this PR. How is it different from the PR for tree canopy dichotomized at the median (Part 1, Step 3), and why?

💡 Answer — click to reveal after you’ve thought about it

The prevalence of high blood pressure among participants in the top quartile of tree canopy was 0.60 times the prevalence among participants in the lower three quartiles combined (95% CI: 0.49, 0.72) — about 40% lower.

This is further from the null than the median-dichotomized PR (0.72) because the comparison is different. Dichotomizing at the median puts Q3 in the “high canopy” group, even though Q3’s prevalence is close to Q1 and Q2’s; that dilutes the contrast. Comparing Q4 with Q1–Q3 separates the groups where the difference actually appears.

Optional: relaxing linearity with a spline

You will not be required to use splines in this course, but they are common in environmental epidemiology (for example, temperature–mortality curves). A natural cubic spline lets the relationship bend smoothly instead of forcing a straight line. df sets how flexible the curve is.

Why doesn’t the “linear” model look like a straight line?

In the first plot, the curve labeled “Linear” bends. That’s because the log-binomial model is linear in the log of prevalence, not in prevalence itself: log(P) = β₀ + β₁ × canopy. Undoing the log gives P = exp(β₀ + β₁ × canopy), which is a curve — each additional percentage point of canopy multiplies prevalence by the same amount (the PR), rather than subtracting the same amount. In the second plot, with the y-axis on a log scale, the “Linear” model is a straight line. The linearity assumption we are checking is always about the scale the model works on — here, the log scale.

Question to consider: the spline and the linear model have almost identical AIC values. What does that suggest about how much the linearity assumption actually matters for this particular exposure and outcome?

💡 Answer — click to reveal after you’ve thought about it

The spline allows the curve to bend wherever the data suggest, yet it fits the data no better than a single straight line on the log scale (the AICs differ by less than 1). The two curves differ mainly at the lowest canopy values, where the spline levels off, and slightly at the highest values, where few participants live. So, for tree canopy and high blood pressure, the linearity assumption appears to be a reasonable approximation, and a continuous term (per percentage point of canopy) is a defensible choice. The quartile plot is still valuable: it shows where along the canopy range the association is concentrated.

🧠 Your Turn: Check linearity for your own exposure

Using your own project’s exposure and outcome:

Create quartiles (or quintiles, if your sample size allows) of your exposure with ntile(), wrapped in factor(), and look at the range of your exposure within each category.
Fit your crude model with the categorical exposure and extract the estimates.
Plot the estimates by category, with the reference category fixed at 1 (or at 0 for a linear model — and use a regular, not log, y-axis).
In two or three sentences, state whether the linearity assumption looks reasonable and what coding you will use in your models going forward.

💡 Example: traffic density (ctb1traff) and asthma — click to reveal

Compared with Q1 (the lowest-traffic quarter of participants), asthma prevalence is similar in Q2 (PR 1.09), somewhat higher in Q3 (PR 1.35; CI includes 1), and more than twice as high in Q4 (PR 2.41; 95% CI: 1.65, 3.53). The rise is steepest at the top, so a continuous term may not capture the shape well; the quartile categories, or the top quartile vs. the rest, are reasonable choices.

Quick-Start Script for Your Own Project

💻 In your RStudio script

# ---- Setup ----
library(tidyverse)

# ---- Read in your processed data ----
data <- read.csv(file.path("data_processed", "class_data_processed.csv"))

# ---- Helper: exponentiate coefficients and 95% CIs ----
extract_ratio <- function(model, label) {
  b  <- coef(model)[-1]
  se <- sqrt(diag(vcov(model)))[-1]
  data.frame(Model = label, Term = names(b), Estimate = exp(b),
             CI_Lower = exp(b - 1.96 * se), CI_Upper = exp(b + 1.96 * se),
             row.names = NULL)
}

# ============================================================
# REPLACE: swap in your own exposure and outcome variable names
# ============================================================

# ---- Outcome: binary or continuous? rare or common? ----
table(data$YOUR_OUTCOME, useNA = "ifany")
prop.table(table(data$YOUR_OUTCOME))

# ---- Look at your exposure's coding and distribution ----
summary(data$YOUR_EXPOSURE)
hist(data$YOUR_EXPOSURE)

# ---- OPTIONAL: rescale a continuous exposure ----
# Only needed if a one-unit change isn't meaningful for your exposure. For a proportion
# (0-1, like tree canopy), multiplying by 100 converts it to a percentage, so one unit =
# one percentage point.
# Choose the scaling based on your exposure's units and distribution -- or skip this step
# and use YOUR_EXPOSURE directly below.
data <- data %>% mutate(YOUR_EXPOSURE_scaled = YOUR_EXPOSURE * 100)   # adjust the scaling

# ---- Crude model, continuous exposure ----
model_cont <- glm(YOUR_OUTCOME ~ YOUR_EXPOSURE_scaled, data = data,
                  family = binomial(link = "log"))   # continuous outcome: lm(YOUR_OUTCOME ~ ...)
extract_ratio(model_cont, "Continuous")

# ---- Linearity check: quartiles as indicator terms ----
data <- data %>% mutate(YOUR_EXPOSURE_quartile = factor(ntile(YOUR_EXPOSURE, 4)))

data %>%                                   # exposure range in each quartile
  group_by(YOUR_EXPOSURE_quartile) %>%
  summarise(n = n(), min = min(YOUR_EXPOSURE, na.rm = TRUE), max = max(YOUR_EXPOSURE, na.rm = TRUE))

model_quart <- glm(YOUR_OUTCOME ~ YOUR_EXPOSURE_quartile, data = data,
                   family = binomial(link = "log"))

results <- extract_ratio(model_quart, "Quartiles") %>%
  mutate(YOUR_EXPOSURE_quartile = factor(2:4)) %>%
  bind_rows(data.frame(YOUR_EXPOSURE_quartile = factor(1), Estimate = 1,
                       CI_Lower = 1, CI_Upper = 1)) %>%
  mutate(YOUR_EXPOSURE_quartile = factor(YOUR_EXPOSURE_quartile, levels = 1:4)) %>%   # put Q1 first
  arrange(YOUR_EXPOSURE_quartile)

ggplot(results, aes(x = YOUR_EXPOSURE_quartile, y = Estimate, group = 1)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.15) +
  geom_point(size = 3) +
  geom_line() +
  scale_x_discrete(labels = c("Q1", "Q2", "Q3", "Q4")) +
  scale_y_log10() +
  labs(x = "[Exposure] quartile", y = "Prevalence ratio (95% CI)",
       title = "[Outcome] by [exposure] quartile (reference = Q1)") +
  theme_minimal()

Wrap Up

This tutorial covered:

Choosing a model from the outcome type, and preferring the log-binomial model (prevalence ratio) over logistic regression (odds ratio) when a binary outcome is common
Reproducing a hand-calculated two-by-two table estimate with a regression model, and why we exponentiate coefficients and confidence intervals to get ratio measures
Rescaling a continuous exposure so that “one unit” is a meaningful comparison
Checking the linearity assumption with categorical indicator terms, interpreting each category’s estimate relative to the reference category, and deciding whether to keep a continuous term, collapse categories, or use a spline

The crude model and linearity check you completed for your own project are the core of Milestone 2.

Next week: Confounding and causal diagrams — every model in this tutorial was crude (unadjusted). Next week you’ll use a DAG to decide which variables belong in the model alongside your exposure, and see how adjusting for them changes your estimate.