This session contains a complete, self-contained set of R commands that
simulate a data set on lung cancer status and six explanatory variables.
assemble the variables into a data frame with properly labelled factors, and
fit a binary logistic regression model whose parameters are reported as odds ratios with 95% confidence intervals.
Because the data are simulated from a model with known parameters, the estimated odds ratios can be compared directly with the true values used to generate the data.
This makes the example useful for teaching, for testing analysis code before real data arrive, and for sample-size or power exploration.
All output shown in this document was produced by running the code in R version 4.6.1.
With set.seed(12345) the results are exactly
reproducible on any machine.
The outcome is binary (lung cancer: No / Yes), so a logistic
regression model is used. On the logit scale the model is:
logit[ P(lung cancer = Yes) ] = b0 + b1(smoking) + b2(sex) + b3(residence) + b4(married) + b5(mining) + b6(age)
Exponentiating a coefficient gives the odds ratio for that variable: exp(b1) is the odds of lung cancer among smokers relative to non-smokers, holding the other variables constant.
For the numeric variable age, exp(b6) is the odds ratio for a one-year increase in age.
Table 1: Description of the variables
| Variable | Type | Categories | Reference | True OR used |
|---|---|---|---|---|
| lung_cancer | Factor | No, Yes | No | - |
| smoking | Factor | No, Yes | No | 4.0 |
| sex | Factor | Female, Male | Female | 1.6 |
| residence | Factor | Rural, Urban | Rural | 1.3 |
| married | Factor | No, Yes | No | 0.8 |
| mining | Factor | No, Yes | No | 2.5 |
| age | Numeric | Years (18–85) | - | 1.05 per year |
In R, the first level of a factor is automatically the reference category, so the levels are deliberately ordered No, Female and Rural first.
The seed makes the simulation reproducible; n is the number of individuals.
set.seed(12345) # makes the simulation reproducible
n <- 3000 # number of individuals to simulateAge is drawn from a normal distribution and trimmed to a plausible
adult range. The five binary variables are drawn from Bernoulli
distributions using rbinom() with size = 1; the probability
supplied is the prevalence of the “Yes” category.
# Age: numeric, mean 45 years, SD 12, restricted to 18-85
age <- round(rnorm(n, mean = 45, sd = 12))
age[age < 18] <- 18
age[age > 85] <- 85
# Binary regressors, coded 0/1 for now (labels are added in Step 4)
smoking <- rbinom(n, 1, 0.45) # 45% smokers
sex <- rbinom(n, 1, 0.50) # 50% male
residence <- rbinom(n, 1, 0.40) # 40% urban
married <- rbinom(n, 1, 0.60) # 60% married
mining <- rbinom(n, 1, 0.15) # 15% work in miningThe linear predictor is built from the true log-odds (the logarithms of the odds ratios in Table 1), converted to a probability with the inverse-logit function, and the outcome is then drawn from a Bernoulli distribution with that probability.
# True regression coefficients = log(odds ratio)
b0 <- -5.5 # intercept, controls overall prevalence
b_smoking <- log(4.0) # OR = 4.0 for smokers
b_sex <- log(1.6) # OR = 1.6 for males
b_resid <- log(1.3) # OR = 1.3 for urban residents
b_married <- log(0.8) # OR = 0.8 for married
b_mining <- log(2.5) # OR = 2.5 for mine workers
b_age <- log(1.05) # OR = 1.05 per additional year
# Linear predictor
eta <- b0 + b_smoking * smoking + b_sex * sex + b_resid * residence +
b_married * married + b_mining * mining + b_age * age
# Convert to a probability and draw the outcome
p <- 1 / (1 + exp(-eta)) # same as plogis(eta)
lung_cancer <- rbinom(n, 1, p)factor() converts the 0/1 codes into labelled
categorical variables. The order given in levels = decides
the reference category, which is what the odds ratios will be compared
against. Age is left as numeric.
lung <- data.frame(lung_cancer = factor(lung_cancer, levels = c(0, 1), labels = c("No", "Yes")),
smoking = factor(smoking, levels = c(0, 1), labels = c("No", "Yes")),
sex = factor(sex, levels = c(0, 1), labels = c("Female", "Male")),
residence = factor(residence, levels = c(0, 1), labels = c("Rural", "Urban")),
married = factor(married, levels = c(0, 1), labels = c("No", "Yes")),
mining = factor(mining, levels = c(0, 1), labels = c("No", "Yes")),
age = age
)str(lung)## 'data.frame': 3000 obs. of 7 variables:
## $ lung_cancer: Factor w/ 2 levels "No","Yes": 2 1 1 1 1 1 1 1 1 1 ...
## $ smoking : Factor w/ 2 levels "No","Yes": 2 1 1 1 1 2 1 2 2 1 ...
## $ sex : Factor w/ 2 levels "Female","Male": 2 2 1 1 2 1 1 2 2 1 ...
## $ residence : Factor w/ 2 levels "Rural","Urban": 1 2 1 2 1 1 2 1 1 1 ...
## $ married : Factor w/ 2 levels "No","Yes": 2 2 1 2 1 2 2 1 1 1 ...
## $ mining : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
## $ age : num 52 54 44 40 52 23 53 42 42 34 ...head(lung, n = 5)summary(lung)## lung_cancer smoking sex residence married mining
## No :2670 No :1669 Female:1490 Rural:1820 No :1171 No :2550
## Yes: 330 Yes:1331 Male :1510 Urban:1180 Yes:1829 Yes: 450
##
##
##
##
## age
## Min. :18.00
## 1st Qu.:37.00
## Median :45.00
## Mean :45.02
## 3rd Qu.:53.00
## Max. :85.00table(lung$lung_cancer, lung$smoking, dnn = c("Lung cancer", "Smoking"))## Smoking
## Lung cancer No Yes
## No 1570 1100
## Yes 99 231table(lung$lung_cancer, lung$mining, dnn = c("Lung cancer", "Mining"))## Mining
## Lung cancer No Yes
## No 2308 362
## Yes 242 88glm() with
family = binomial(link = "logit") fits the model. Because
lung_cancer is a factor with levels No then Yes, R models the
probability of the second level, Yes, which is what is wanted.
model <- glm(lung_cancer ~ smoking + sex + residence + married + mining + age,
data = lung, family = binomial(link = "logit"))
summary(model)##
## Call:
## glm(formula = lung_cancer ~ smoking + sex + residence + married +
## mining + age, family = binomial(link = "logit"), data = lung)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.908473 0.328692 -17.976 < 2e-16 ***
## smokingYes 1.238108 0.130794 9.466 < 2e-16 ***
## sexMale 0.423950 0.125319 3.383 0.000717 ***
## residenceUrban 0.221638 0.124875 1.775 0.075918 .
## marriedYes -0.274972 0.124859 -2.202 0.027648 *
## miningYes 0.898748 0.147135 6.108 1.01e-09 ***
## age 0.058129 0.005588 10.402 < 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: 2079.1 on 2999 degrees of freedom
## Residual deviance: 1806.9 on 2993 degrees of freedom
## AIC: 1820.9
##
## Number of Fisher Scoring iterations: 5The Estimate column contains log-odds, not odds ratios. The next step converts them.
Exponentiating the coefficients and their confidence limits gives the odds ratios.
confint.default() returns Wald intervals;
confint() returns profile-likelihood intervals,
which are preferable in small samples but slower.
# Quick version: odds ratios with 95% Wald confidence intervals
exp(cbind(OR = coef(model), confint.default(model)))## OR 2.5 % 97.5 %
## (Intercept) 0.002716333 0.001426264 0.005173278
## smokingYes 3.449080933 2.669141768 4.456922979
## sexMale 1.527984959 1.195218266 1.953398890
## residenceUrban 1.248120090 0.977152480 1.594227913
## marriedYes 0.759593703 0.594704071 0.970201183
## miningYes 2.456526698 1.841109998 3.277655016
## age 1.059851851 1.048307386 1.071523449# Tidy table with p-values
OR <- exp(coef(model))
CI <- exp(confint.default(model))
pval <- summary(model)$coefficients[, 4]
or_table <- data.frame(
Term = names(OR),
OR = round(OR, 3),
Lower_95 = round(CI[, 1], 3),
Upper_95 = round(CI[, 2], 3),
p_value = round(pval, 4),
row.names = NULL
)
print(or_table)## Term OR Lower_95 Upper_95 p_value
## 1 (Intercept) 0.003 0.001 0.005 0.0000
## 2 smokingYes 3.449 2.669 4.457 0.0000
## 3 sexMale 1.528 1.195 1.953 0.0007
## 4 residenceUrban 1.248 0.977 1.594 0.0759
## 5 marriedYes 0.760 0.595 0.970 0.0276
## 6 miningYes 2.457 1.841 3.278 0.0000
## 7 age 1.060 1.048 1.072 0.0000Smokers have about 3.4 times the odds of lung cancer compared with non-smokers, after adjusting for the other variables (OR = 3.45, 95% CI 2.67–4.46).
Working in mining is associated with about 2.5 times the odds of lung cancer relative to not working in mining (OR = 2.46, 95% CI 1.84–3.28).
Males have roughly 53% higher odds than females (OR = 1.53, 95% CI 1.20–1.95).
Each additional year of age raises the odds by about 6% (OR = 1.06, 95% CI 1.05–1.07). Multiplying over 10 years gives 1.06^10 ≈ 1.78, i.e. 78% higher odds per decade.
Being married is associated with lower odds (OR = 0.76, 95% CI 0.60–0.97); the confidence interval excludes 1, so the association is statistically significant at the 5% level.
Urban residence shows a raised but non-significant odds ratio (OR = 1.25, 95% CI 0.98–1.59, p = 0.076). The interval includes 1, so with this sample size the effect cannot be distinguished from no effect even though a true effect of 1.30 was built into the simulation — a useful reminder about power.
We perform overall model fit to see whether the fitted model is significant (i.e. the regressors included in the model explains the response).
lr <- model$null.deviance - model$deviance # likelihood-ratio statistic
pchisq(lr, df = 6, lower.tail = FALSE) # p-value for the whole model## [1] 7.5042e-56AIC(model)## [1] 1820.937The likelihood-ratio chi-square is significant (p < 0.001), so the six predictors jointly explain a significant amount of variation in lung cancer status. AIC = 1820.94.
Changing the seed: remove or change set.seed(12345)
to get a different simulated sample. Re-running with the same seed
always reproduces these exact numbers.
Changing the sample size or prevalence: adjust n, or adjust the intercept b0. A more negative intercept gives a rarer outcome.
Changing the reference category: use relevel(), for
example
lung$residence <- relevel(lung$residence, ref = "Urban").
Profile-likelihood intervals: replace
confint.default(model) with confint(model).
These are generally preferred, especially with few events.
Tidy output with an add-on package:
install.packages("broom"), then
broom::tidy(model, exponentiate = TRUE, conf.int = TRUE)
produces the odds-ratio table in one line.
Age as a categorical variable:
cut(age, breaks = c(17, 39, 59, 85), labels = c("18-39", "40-59", "60+"))
creates age groups if categorical age is preferred over a linear
term.
Checking the fitted model: predicted probabilities come from
predict(model, type = "response"); classification accuracy
and the ROC curve can be examined with the pROC package.