1 Why Logistic Regression Appears in the Lectures

The causal inference lectures use logistic regression for two different but related tasks.

  1. Treatment model, or propensity score model: estimate \(e(X)=P(T=1\mid X)\), the probability of being a smoker given baseline covariates.
  2. Outcome model: estimate \(m(t,X)=P(Y=1\mid T=t,X)\), the probability of CVD under a fixed smoking status and baseline covariates.

Both are logistic regressions because both targets are probabilities for binary variables.

Tutorial Goal

This note explains logistic regression as a probability model, connects its coefficients to odds ratios, and shows how the fitted model is used in the lectures for propensity scores and outcome-regression standardization. The focus is the classroom NHANES complete-case analysis; NHANES survey weights are not used here.

2 Binary Outcomes Need Probability Models

For a binary outcome \(Y\in\{0,1\}\), the conditional mean is a probability:

\[ E(Y\mid X=x)=P(Y=1\mid X=x). \]

If the outcome is CVD, then \(Y=1\) means that CVD was reported before the interview and \(Y=0\) means it was not. A regression model for this outcome must return values between 0 and 1.

A plain linear model,

\[ P(Y=1\mid X=x)\approx \alpha+\beta x, \]

can predict probabilities below 0 or above 1. Logistic regression avoids this by modeling the log odds, which can take any real value, and then mapping back to a probability.

3 Odds, Log Odds, and the Logistic Curve

For a probability \(p\), the odds are

\[ \text{odds}=\frac{p}{1-p}. \]

The logit transform is the log odds:

\[ \operatorname{logit}(p)=\log\left(\frac{p}{1-p}\right). \]

The inverse transform is the logistic, or expit, function:

\[ p = \operatorname{expit}(\eta) = \frac{\exp(\eta)}{1+\exp(\eta)} = \frac{1}{1+\exp(-\eta)}. \]

prob_grid <- c(0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90)

odds_table <- data.frame(
  Probability = prob_grid,
  Odds = round(prob_grid / (1 - prob_grid), 3),
  Log_odds = round(logit(prob_grid), 3)
)

knitr::kable(odds_table)
Probability Odds Log_odds
0.01 0.010 -4.595
0.05 0.053 -2.944
0.10 0.111 -2.197
0.25 0.333 -1.099
0.50 1.000 0.000
0.75 3.000 1.099
0.90 9.000 2.197

The probability scale is bounded, while the log-odds scale is unbounded. This is why logistic regression can use a linear predictor without producing impossible probabilities.

4 The Logistic Regression Model

A logistic regression model has the form

\[ \operatorname{logit}\{P(Y=1\mid X)\} = \beta_0+\beta_1X_1+\cdots+\beta_pX_p. \]

Equivalently,

\[ P(Y=1\mid X) = \operatorname{expit}(\beta_0+\beta_1X_1+\cdots+\beta_pX_p). \]

The coefficient \(\beta_j\) is a change in conditional log odds. Its exponential, \(\exp(\beta_j)\), is a conditional odds ratio for a one-unit increase in \(X_j\), holding the other model covariates fixed.

Important Distinction

A logistic regression coefficient is not automatically a marginal causal effect. In the lectures, logistic regression is mainly used to estimate nuisance functions: \(e(X)\) for IPW and \(m(t,X)\) for outcome regression. The causal targets are marginal risks such as \(E\{Y(1)\}\), \(E\{Y(0)\}\), risk differences, and risk ratios after standardization or weighting.

5 Load the Classroom NHANES Data

The tutorial uses the same classroom NHANES variables as Lectures 1-6:

Symbol Dataset variable Meaning
\(Y\) cvd_indicator CVD indicator
\(T\) smoker_indicator Lifetime smoker indicator
\(X\) baseline variables age, gender, race, education, income-to-poverty ratio, BMI 
data_candidates <- c(
  file.path("..", "Data", "NHANES_data.csv"),
  file.path("Data", "NHANES_data.csv"),
  "C:/Users/seyoo/OneDrive/GitHub_repository/Causal_Inference_Survival_Analysis/Causal_Inference_Lecture/Data/NHANES_data.csv"
)

data_path <- data_candidates[file.exists(data_candidates)][1]
if (is.na(data_path)) {
  stop("Could not find NHANES_data.csv. Check the data path.")
}

nhanes <- read.csv(data_path, stringsAsFactors = FALSE)

baseline_vars <- c(
  "age_yr",
  "gender",
  "race",
  "educ_lvl",
  "inc_to_pov_ratio",
  "bmi"
)

analysis_vars <- c("smoker_indicator", "cvd_indicator", baseline_vars)

tutorial_data <- nhanes[complete.cases(nhanes[, analysis_vars]), analysis_vars]

tutorial_data$smoker_indicator <- as.integer(tutorial_data$smoker_indicator)
tutorial_data$cvd_indicator <- as.integer(tutorial_data$cvd_indicator)

tutorial_data$smoking_group <- ifelse(
  tutorial_data$smoker_indicator == 1,
  "Smoker",
  "Non-smoker"
)
tutorial_data$smoking_group <- factor(
  tutorial_data$smoking_group,
  levels = c("Non-smoker", "Smoker")
)

tutorial_data$gender_label <- factor(
  tutorial_data$gender,
  levels = c("F", "M"),
  labels = c("Female", "Male")
)

tutorial_data$race_label <- factor(
  tutorial_data$race,
  levels = c(
    "mex_american",
    "other_hispanic",
    "nh_white",
    "nh_black",
    "nh_asian"
  ),
  labels = c(
    "Mexican American",
    "Other Hispanic",
    "Non-Hispanic White",
    "Non-Hispanic Black",
    "Non-Hispanic Asian"
  )
)

tutorial_data$educ_label <- factor(
  tutorial_data$educ_lvl,
  levels = c(
    "lt_9th_grade",
    "9_to_12th_grade_no_diploma",
    "hs_grad_or_ged",
    "start_college_to_aa",
    "college_grad_or_above"
  ),
  labels = c(
    "Less than 9th grade",
    "9th-12th, no diploma",
    "High school or GED",
    "Some college or AA",
    "College graduate or above"
  )
)

nrow(tutorial_data)
## [1] 6299
sample_summary <- data.frame(
  Quantity = c(
    "Complete-case sample size",
    "Number of smokers",
    "Number of non-smokers",
    "Number with CVD",
    "Observed smoker risk (%)",
    "Observed non-smoker risk (%)"
  ),
  Value = c(
    nrow(tutorial_data),
    sum(tutorial_data$smoker_indicator == 1),
    sum(tutorial_data$smoker_indicator == 0),
    sum(tutorial_data$cvd_indicator == 1),
    round(100 * mean(tutorial_data$cvd_indicator[tutorial_data$smoker_indicator == 1]), 2),
    round(100 * mean(tutorial_data$cvd_indicator[tutorial_data$smoker_indicator == 0]), 2)
  )
)

knitr::kable(sample_summary)
Quantity Value
Complete-case sample size 6299.00
Number of smokers 2559.00
Number of non-smokers 3740.00
Number with CVD 646.00
Observed smoker risk (%) 15.12
Observed non-smoker risk (%) 6.93

6 A One-Predictor Logistic Regression: Age and CVD

Before using the full lecture model, it helps to see one simple logistic regression:

\[ \operatorname{logit}\{P(\text{CVD}=1\mid \text{age})\} = \beta_0+\beta_1\text{age}. \]

age_model <- glm(
  cvd_indicator ~ age_yr,
  data = tutorial_data,
  family = binomial()
)

summary(age_model)
## 
## Call:
## glm(formula = cvd_indicator ~ age_yr, family = binomial(), data = tutorial_data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.172083   0.220240  -28.02   <2e-16 ***
## age_yr       0.072471   0.003561   20.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4165.7  on 6298  degrees of freedom
## Residual deviance: 3605.6  on 6297  degrees of freedom
## AIC: 3609.6
## 
## Number of Fisher Scoring iterations: 6

The gray points are individual binary outcomes. Each participant contributes a 0 or a 1, so the raw data live on two horizontal bands. The red points compress those binary outcomes into local empirical risks within age bins. The blue curve is the logistic model’s estimate of \[ P(\text{CVD}=1\mid \text{age}), \] so it should be read as a conditional mean curve, not as a deterministic prediction of who will have CVD. At age \(a\), the curve answers: among participants with similar age \(a\), what fraction would the fitted model expect to have CVD?

How to Read the Age-Risk Plot

The plot moves through three levels of aggregation. First, the gray points show individual binary outcomes. Second, the red points show empirical averages within age ranges. Third, the blue logistic curve smooths those empirical averages by imposing a linear relation on the log-odds scale. The curve is useful because causal estimators need predicted probabilities, but the curve is still a working statistical model.

7 Interpreting Coefficients

The age coefficient is a log-odds coefficient. Exponentiating it gives an odds ratio.

age_beta <- coef(age_model)["age_yr"]
age_or <- exp(age_beta)
age_or_10 <- exp(10 * age_beta)

age_interpretation <- data.frame(
  Quantity = c(
    "Age coefficient",
    "Odds ratio for one additional year",
    "Odds ratio for ten additional years"
  ),
  Value = round(c(age_beta, age_or, age_or_10), 3)
)

knitr::kable(age_interpretation)
Quantity Value
Age coefficient 0.072
Odds ratio for one additional year 1.075
Odds ratio for ten additional years 2.064

For a one-year increase in age, the model multiplies the conditional odds of CVD by 1.075. For a ten-year increase, it multiplies the conditional odds by 2.064.

Odds Ratio Versus Risk Ratio

An odds ratio is not the same as a risk ratio. When the outcome is rare, they can be close. When the outcome is not rare, they can differ substantially. The causal lectures therefore report marginal risks, risk differences, and risk ratios after standardization or weighting, rather than treating the logistic coefficient itself as the final causal estimand.

8 The Outcome Regression Model Used in the Lecture

Lecture 3 used logistic regression to estimate the conditional CVD risk:

\[ m(t,X)=P(Y=1\mid T=t,X). \]

The working model was:

\[ \operatorname{logit}\{P(Y=1\mid T,X)\} = \beta_0+\beta_1T+\beta_2\text{age}+\cdots. \]

outcome_model <- glm(
  cvd_indicator ~ smoker_indicator + age_yr + gender_label +
    race_label + educ_label + inc_to_pov_ratio + bmi,
  data = tutorial_data,
  family = binomial()
)

outcome_coef <- summary(outcome_model)$coefficients

main_terms <- c("smoker_indicator", "age_yr", "bmi", "inc_to_pov_ratio")

outcome_or_table <- data.frame(
  Term = main_terms,
  Log_odds_coefficient = round(outcome_coef[main_terms, "Estimate"], 3),
  Conditional_odds_ratio = round(exp(outcome_coef[main_terms, "Estimate"]), 3),
  P_value = signif(outcome_coef[main_terms, "Pr(>|z|)"], 3)
)

knitr::kable(outcome_or_table)
Term Log_odds_coefficient Conditional_odds_ratio P_value
smoker_indicator smoker_indicator 0.431 1.539 5.5e-06
age_yr age_yr 0.072 1.075 0.0e+00
bmi bmi 0.030 1.031 1.0e-07
inc_to_pov_ratio inc_to_pov_ratio -0.183 0.833 0.0e+00

The coefficient for smoker_indicator compares smokers and non-smokers conditional on the model covariates. The standardized risk comparison below is a different object: it averages predicted risks over the same empirical covariate distribution.

The curves hold the other covariates fixed at reference values. This picture is useful for understanding the fitted model, but the lecture’s causal estimand is marginal over the whole classroom covariate distribution.

9 From Logistic Model to Standardized Risks

Outcome-regression standardization uses the fitted logistic model in two counterfactual prediction steps:

  1. Copy the analytic sample and set everyone to \(T=0\).
  2. Copy the analytic sample and set everyone to \(T=1\).
  3. Predict each person’s risk in both copied datasets.
  4. Average the predicted risks over the same empirical covariate distribution.

data_if_non_smoker <- tutorial_data
data_if_non_smoker$smoker_indicator <- 0

data_if_smoker <- tutorial_data
data_if_smoker$smoker_indicator <- 1

predicted_risk_if_non_smoker <- predict(
  outcome_model,
  newdata = data_if_non_smoker,
  type = "response"
)

predicted_risk_if_smoker <- predict(
  outcome_model,
  newdata = data_if_smoker,
  type = "response"
)

or_risk_non_smoker <- mean(predicted_risk_if_non_smoker)
or_risk_smoker <- mean(predicted_risk_if_smoker)
or_risk_difference <- or_risk_smoker - or_risk_non_smoker
or_risk_ratio <- or_risk_smoker / or_risk_non_smoker

standardized_table <- data.frame(
  Quantity = c(
    "Standardized risk if non-smoker (%)",
    "Standardized risk if smoker (%)",
    "Risk difference (percentage points)",
    "Risk ratio"
  ),
  Value = c(
    round(100 * or_risk_non_smoker, 2),
    round(100 * or_risk_smoker, 2),
    round(100 * or_risk_difference, 2),
    round(or_risk_ratio, 3)
  )
)

knitr::kable(standardized_table)
Quantity Value
Standardized risk if non-smoker (%) 8.530
Standardized risk if smoker (%) 11.980
Risk difference (percentage points) 3.450
Risk ratio 1.405

The density plot is not showing the observed CVD outcome. It is showing two columns of fitted counterfactual risks. Every participant is predicted twice using the same covariates \(X_i\): once as if \(T_i=0\), giving \(\hat m(0,X_i)\), and once as if \(T_i=1\), giving \(\hat m(1,X_i)\). The solid curves show the distribution of those person-specific predicted risks across the analytic sample.

The dashed vertical lines are not modes, medians, or cutpoints. They are the averages of the two prediction columns: \[ \hat\psi_0^{OR}=\frac{1}{n}\sum_{i=1}^n \hat m(0,X_i), \qquad \hat\psi_1^{OR}=\frac{1}{n}\sum_{i=1}^n \hat m(1,X_i). \] Here the sum is over all analytic-sample participants because the target is the analytic-sample marginal risk, \[ \psi_t = E_X\{m(t,X)\}, \] where the outer expectation is over the target covariate distribution \(P_X\), not over the smoker-specific distribution \(P_{X\mid T=1}\) or the non-smoker-specific distribution \(P_{X\mid T=0}\). Its empirical plug-in version replaces \(P_X\) by the empirical distribution of all observed \(X_i\)’s.

The blue dashed line is therefore the average predicted risk if everyone in the analytic sample were set to non-smoking; the orange dashed line is the average predicted risk if everyone in the analytic sample were set to smoking. Their horizontal distance is the standardized risk difference, \[ \hat\Delta^{OR} = \hat\psi_1^{OR}-\hat\psi_0^{OR}. \] Thus, standardization turns conditional probabilities into marginal risks by averaging both columns over the same target empirical covariate distribution.

If instead one averaged only within the observed group, \[ \frac{1}{n_t}\sum_{i:T_i=t}\hat m(t,X_i), \] then the averaging distribution would be \(P_{X\mid T=t}\). That is a different object: it describes the fitted mean risk within the covariate distribution of the observed group \(T=t\). It is not the analytic-sample marginal counterfactual risk unless the group covariate distribution happens to equal the full analytic-sample covariate distribution. For an ATT target, for example, one would intentionally average over the treated covariate distribution, such as \[ \hat\psi_0^{ATT,OR} = \frac{1}{n_1}\sum_{i:T_i=1}\hat m(0,X_i), \] but that is a different estimand from the marginal risk comparison used in this tutorial.

10 Logistic Regression for Propensity Scores

Lecture 2 used logistic regression for a different binary outcome: smoking status.

\[ e(X)=P(T=1\mid X). \]

The fitted values from this model are propensity scores. They are probabilities, so they must lie between 0 and 1.

propensity_model <- glm(
  smoker_indicator ~ age_yr + gender_label + race_label +
    educ_label + inc_to_pov_ratio + bmi,
  data = tutorial_data,
  family = binomial()
)

tutorial_data$propensity_score <- predict(propensity_model, type = "response")

summary(tutorial_data$propensity_score)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05416 0.26826 0.39915 0.40625 0.53934 0.85971

The propensity score plot compares the fitted distribution of \(\hat e(X)\) among observed smokers and non-smokers. The middle region, where both colors appear, is the empirical overlap region: people with similar covariates appear in both treatment groups. The tails matter because IPW uses \[ \frac{T_i}{\hat e(X_i)} \qquad\text{and}\qquad \frac{1-T_i}{1-\hat e(X_i)}. \] Therefore, smokers with very small \(\hat e(X_i)\) receive large treated weights, and non-smokers with very large \(\hat e(X_i)\) receive large untreated weights. A plot with little overlap would signal a practical positivity problem.

11 Model Fit: Calibration on the Risk Scale

A logistic model can be useful even when individual predictions are noisy. For teaching purposes, one simple diagnostic is calibration: among people assigned similar predicted risk, do the observed risks roughly match the predicted risks?

tutorial_data$predicted_cvd_risk <- predict(outcome_model, type = "response")

cal_breaks <- unique(quantile(
  tutorial_data$predicted_cvd_risk,
  probs = seq(0, 1, by = 0.1),
  na.rm = TRUE
))

tutorial_data$risk_decile <- cut(
  tutorial_data$predicted_cvd_risk,
  breaks = cal_breaks,
  include.lowest = TRUE
)

calibration_table <- aggregate(
  cbind(predicted_cvd_risk, cvd_indicator) ~ risk_decile,
  data = tutorial_data,
  FUN = mean
)

calibration_table$n <- as.integer(table(tutorial_data$risk_decile)[
  as.character(calibration_table$risk_decile)
])

names(calibration_table) <- c(
  "Risk_decile",
  "Mean_predicted_risk",
  "Observed_risk",
  "N"
)

knitr::kable(
  data.frame(
    Risk_decile = calibration_table$Risk_decile,
    Mean_predicted_percent = round(100 * calibration_table$Mean_predicted_risk, 2),
    Observed_percent = round(100 * calibration_table$Observed_risk, 2),
    N = calibration_table$N
  )
)
Risk_decile Mean_predicted_percent Observed_percent N
[0.00205,0.00924] 0.62 0.16 630
(0.00924,0.0162] 1.26 1.59 630
(0.0162,0.0257] 2.07 1.27 630
(0.0257,0.0382] 3.14 1.75 630
(0.0382,0.0569] 4.72 5.24 630
(0.0569,0.0848] 7.07 7.63 629
(0.0848,0.122] 10.21 12.06 630
(0.122,0.181] 14.90 15.56 630
(0.181,0.271] 22.20 23.33 630
(0.271,0.674] 36.36 33.97 630

The calibration plot groups participants by predicted-risk decile. For each decile, the \(x\)-coordinate is the average predicted CVD risk and the \(y\)-coordinate is the observed fraction with CVD. Points close to the diagonal line mean that the model’s average prediction in that risk stratum is close to the empirical event rate. Points above the line indicate underprediction; points below the line indicate overprediction.

This diagnostic does not prove the model is correct, and it does not prove the causal assumptions. It only checks whether predicted risks and observed risks are broadly aligned across risk strata.

12 What Logistic Regression Does Not Solve

Logistic regression is a statistical model. It does not by itself guarantee causal identification.

The lectures use logistic regression in a deliberately simple way: main effects for measured baseline covariates. That makes the tutorial transparent, but it is still a working model. If the true response surface has nonlinearities or interactions, the fitted probabilities may be misspecified. Doubly robust methods help with some nuisance-model misspecification, but they do not repair violations of the causal identification assumptions.

13 Summary

Logistic regression is useful in these lectures because it estimates probabilities for binary variables.

  1. On the log-odds scale, the model is linear.
  2. On the probability scale, predictions stay between 0 and 1.
  3. For the propensity score, logistic regression estimates \(e(X)=P(T=1\mid X)\).
  4. For outcome regression, logistic regression estimates \(m(t,X)=P(Y=1\mid T=t,X)\).
  5. For causal risk comparisons, the coefficient of \(T\) is not the final target; standardized or weighted marginal risks are the target.

The practical lesson is:

\[ \text{logistic regression estimates conditional probabilities;} \qquad \text{causal estimators transform them into marginal causal contrasts.} \]