Module 1 · Homework 2

Worked answers with reasoning. References point to the lecture notes PDF and the notes pages. Questions 5–8 use hw2_braincancer2.csv (88 patients, 53 events; the BrainCancer covariates sex, diagnosis, loc, ki, gtv, stereo, with time in months and status).

← Back to notes index

Q1 — When to use the log-rank test

Answer: When we have one (or more) groups defined by a categorical variable and we want to compare their survival curves over time.

The log-rank test compares $S_1(t)$ vs. $S_0(t)$ (or $J$ curves at once) across the whole observed time frame, pooling the observed-minus-expected event counts $U_k = d_{1k} - n_{1k}d_k/n_k$ from the $2 \times 2$ table at every event time. Notes §1.6: “When $Z_i$ is a single variable that is either binary or categorical, we can use the log-rank test. But when $Z_i$ is continuous or represents multiple variables, we must instead consider an approach to model the distribution of $T \mid Z$.”

Why the others are wrong:

Estimating a hazard ratio under proportional hazards — that is the Cox model. The log-rank test yields a $p$-value, not an effect estimate.
Comparing nested models — that is the partial likelihood ratio test, $2(l_a - l_b) \sim \chi^2_{a-b}$ (notes p. 33).
Multiple continuous predictors, adjusted — needs a regression model (Cox); the log-rank test handles neither continuous covariates nor adjustment.

Q2 — Two treatments, controlling for trial site

Answer: Stratified log-rank test.

Two treatment groups, a censored time-to-event outcome, and a single categorical variable (hospital, $M = 12$ strata) to adjust for. Notes p. 22: “The stratified log-rank test is one way to adjust for a single categorical variable when testing for a difference between the two treatment groups.” It forms the $2 \times 2$ table within each stratum at each event time and pools across strata and event times,

\[T_S = \frac{\big[\sum_{m=1}^{M}\sum_{k=1}^{K} U_{k,m}\big]^2} {\sum_{m=1}^{M}\sum_{k=1}^{K} V_{k,m}} \sim \chi^2_1\]

testing $H_0 : S_{1,m}(t) = S_{0,m}(t)$ for all strata $m$ and all $t$. In R: survdiff(Surv(time, relapse) ~ drug + strata(site)). This is also the design the notes describe — randomization stratified by a prognostic factor, then the matching stratified analysis.

Why the others are wrong:

Likelihood ratio test — compares nested models, not two treatment groups.
Paired t-test — paired continuous outcomes with no censoring; patients here are neither paired nor uncensored.
Rank-sum test — uncensored continuous data, no stratification, cannot handle censoring.

Q3 — Correct statement about the Cox model

Answer: Under the proportional hazards assumption, the hazard ratio between two covariate profiles $x_1$ vs. $x_2$ is constant over time and equals $\exp{\beta^T(x_1 - x_2)}$.

From $\lambda(t \mid Z) = \lambda_0(t)\exp{\beta^T Z}$,

\[\frac{\lambda(t \mid x_1)}{\lambda(t \mid x_2)} = \frac{\lambda_0(t)\, e^{\beta^T x_1}}{\lambda_0(t)\, e^{\beta^T x_2}} = \exp\{\beta^T(x_1 - x_2)\}\]

The baseline hazard cancels, so the ratio has no $t$ in it — that is the proportional hazards property (notes p. 27–28, generalized from a one-unit change in a single $Z_j$, where it reduces to $e^{\beta_j}$).

Why the others are wrong:

“$h_0(t)$ must be specified” — backwards. $\lambda_0(t)$ is left unspecified; that is what makes the Cox model semiparametric. Specifying it as exponential or Weibull is parametric (AFT) regression.
“Violations can always be fixed by adding baseline covariates” — false. The remedies in the notes are stratifying the baseline hazard, interaction with $t$ or $\log t$ (a time-varying term), piecewise constant hazards, splines. Time-dependent terms are one of the fixes, not something to avoid.
“Partial likelihood estimates $\beta$ and $\Lambda_0$ simultaneously” — false. $\lambda_0(t_j)$ cancels out of every factor $\exp{\beta^T Z_j} / \sum_{l \in R_j}\exp{\beta^T Z_l}$, so the partial likelihood involves $\beta$ only. $\hat{\Lambda}_0(t)$ comes afterwards from Breslow’s estimator with $\hat\beta$ plugged in.

Q4 — Reading a Cox model table (pharmacoSmoking)

Answer: Holding age, treatment group, employment, and race constant, the estimated hazard of relapse for those with a light level of smoking was 1.17 times the hazard of relapse for those with a heavy smoking level.

Reference levels (R takes the first factor level): grp = combination, employment = ft, race = black, levelSmoking = heavy. Every row in the table is a comparison against those.

Statement	Check	Verdict
A patch-only decreased the probability of relapse by 0.65	$e^{0.65163} = 1.92$ — hazard is higher; and Cox output describes the hazard, never a probability	wrong
B age, patch-only, pt/other employment all decreased the hazard	only age has $\beta < 0$; the other three are $\beta > 0$ (HR $> 1$, harm)	wrong
C light vs. heavy smoking HR $= 1.17$	$e^{0.15678} = 1.170$, reference = heavy, hazard language, “holding the rest constant”	correct
D other vs. part-time employment HR $= 2.27$	$e^{0.81910 - 0.59316} = 1.25$; $2.27$ is other vs. full time (the reference)	wrong

The general rule from the notes (p. 29–30): $e^{\beta_j}$ is the relative change in the hazard for a one-unit change in $Z_j$ (here, one factor level vs. its reference) holding everything else constant. Do not interpret in terms of survival probability or risk. Comparing two non-reference levels needs the difference of their $\beta$s, not one of them alone.

Note C’s HR is not significant ($p = 0.55$), but the statement only describes the point estimate, and that description is accurate. Also, the table’s employmentpt row prints $\exp(\beta) = 2.00003$, yet $e^{0.59316} = 1.81$ — a typo in the question that does not affect the answer.

Script: hw2_q4.R — recomputes every $\exp(\beta)$ and contrast from the table’s values, and refits the model on asaur::pharmacoSmoking to confirm the reference levels (the packaged data give slightly different coefficients, so the question’s table is from another data version).

Q5 — Log-rank test: SRS vs. SRT

Answer: 0.10

bc = read.csv("hw2_braincancer2.csv")
survdiff(Surv(time, status) ~ stereo, data = bc)

            N Observed Expected (O-E)^2/E (O-E)^2/V
stereo=SRS 23       18     12.9     2.007      2.75
stereo=SRT 65       35     40.1     0.646      2.75

 Chisq= 2.7  on 1 degrees of freedom, p= 0.0974

$T_{LR} = \big[\sum_k U_k\big]^2 / \sum_k V_k = 2.75 \sim \chi^2_1$, so $p = P(\chi^2_1 > 2.75) = 0.0974 \to \mathbf{0.10}$.

Reading the table: SRS patients had more deaths than expected under $H_0: S_{SRS}(t) = S_{SRT}(t)$ (18 observed vs. 12.9 expected) and SRT fewer (35 vs. 40.1), but the pooled evidence does not reach significance at 0.05. Note survdiff prints the $p$-value rounded to one digit (p= 0.1); the $pvalue element gives 0.0974.

Script: hw2_q5_8.R — Q5–Q8 share the dataset, so they live in one script.

Q6 — Kaplan–Meier $\hat{S}(24)$ in the SRT group

Answer: 0.61

km.srt = survfit(Surv(time, status) ~ 1, data = bc, subset = (stereo == "SRT"))
summary(km.srt, times = 24)$surv        # 0.6104

KM is a step function that changes only at event times, so $\hat{S}(24)$ is the value at the last SRT death at or before 24 months. The steps around $t = 24$:

time	$n_k$	$d_k$	$\hat{S}_{KM}$
19.555	35	1	0.6498
22.430	33	1	0.6301
23.611	32	1	0.6104
28.854	30	1	0.5901

The last death before 24 is at 23.611 months, so $\hat{S}_{SRT}(24) = 0.6104 \to \mathbf{0.61}$. (The curve does not move again until 28.854.) Both approx(km$time, km$surv, xout = 24, method = "constant") and summary(km, times = 24) return the same value.

Script: hw2_q5_8.R.

Q7 — Cox model: adjusted HR for SRT vs. SRS

Answer: 0.64

bc$stereo = factor(bc$stereo)          # levels SRS, SRT -> SRS is the reference
cox7 = coxph(Surv(time, status) ~ ki + gtv + stereo, data = bc)

               coef exp(coef) se(coef)      z Pr(>|z|)
ki         0.005322    1.0053  0.01269  0.419    0.675
gtv        0.012901    1.0130  0.01699  0.759    0.448
stereoSRT -0.450318    0.6374  0.30258 -1.488    0.137

R takes the first factor level, SRS, as the reference, so the stereoSRT row is already the contrast asked for — SRT relative to SRS:

\[\widehat{HR}_{SRT\,vs\,SRS} = e^{\hat\beta_{SRT}} = e^{-0.4503} = 0.6374 \to \mathbf{0.64}\]

Interpretation: holding Karnofsky index and tumour volume constant, the hazard of death for SRT patients is estimated to be 0.64 times that of SRS patients. The 95% CI $(0.35,\ 1.15)$ includes 1 ($p = 0.14$), in line with the non-significant log-rank test in Q5. Had the question asked for SRS vs. SRT instead, the answer would be $1/0.6374 = 1.57$.

Script: hw2_q5_8.R.

Q8 — Test $H_0: S_{SRS}(30) = S_{SRT}(30)$ with Kaplan–Meier

Answer: 0.66

The single-time-point comparison from the notes (p. 17–18), done exactly as in the Lecture Video 9 exercise: difference of the two KM estimates at $t^* = 30$, variance from Greenwood’s formula in each group (independent groups, so the variances add), then a $z$-test.

km.g = survfit(Surv(time, status) ~ stereo, data = bc)
summary(km.g, times = 30)

Group	$n$ at risk at 30	$\hat{S}(30)$	Greenwood se
SRS	10	0.5349	0.1085
SRT	29	0.5901	0.0666

\[\hat{S}_{SRT}(30) - \hat{S}_{SRS}(30) = 0.5901 - 0.5349 = 0.0552\] \[\text{se} = \sqrt{0.1085^2 + 0.0666^2} = 0.1273 \qquad z = \frac{0.0552}{0.1273} = 0.434\] \[p = 2\,\{1 - \Phi(0.434)\} = 0.6646 \to \mathbf{0.66}\]

Equivalently, the 95% CI $0.0552 \pm 1.96 \times 0.1273 = (-0.19,\ 0.30)$ includes 0, so we fail to reject $H_0$. The sign of the difference does not matter for the two-sided $p$-value.

Contrast with Q5: the log-rank test ($p = 0.10$) uses every event time and is more sensitive than this one-time-point test ($p = 0.66$) — the “inefficient use of the available information” limitation the notes list for single-$t^*$ tests.

Script: hw2_q5_8.R.

Q9 — Hazard ratio between two patients

Answer: 2.39

For any two covariate profiles the Cox model gives (see Q3)

\[\frac{\lambda(t \mid z_A)}{\lambda(t \mid z_B)} = \frac{\lambda_0(t)\,e^{\beta^T z_A}}{\lambda_0(t)\,e^{\beta^T z_B}} = \exp\{\beta^T (z_A - z_B)\}\]

so only the differences in covariates matter and the baseline hazard cancels.

Covariate	$\beta$	$z_A$	$z_B$	$\beta\,(z_A - z_B)$
Age	0.05	35	21	$0.05 \times 14 = 0.70$
Treatment	$-0.10$	1	0	$-0.10$
Sex = Female	0.15	1	0	$0.15$
Ulcer = Present	0.12	1	0	$0.12$
Sum				0.87

\[HR_{A\,vs\,B} = e^{0.87} = 2.387 \to \mathbf{2.39}\]

Patient A’s hazard is about 2.4 times patient B’s, mostly from the 14-year age gap; the treatment (protective, $\beta < 0$) offsets only a little of it. Sanity check on direction: A is older, female, has an ulcer — all $\beta > 0$ — and received treatment ($\beta < 0$), so the net should exceed 1, and it does.

Script: hw2_q9.R — computes the term-by-term contributions and cross-checks $e^{\beta^T z_A} / e^{\beta^T z_B}$.

Q10 — Best model by AIC

Answer: III

\[\mathrm{AIC} = -2\,l_a + 2a\]

(notes p. 33): $-2l_a$ rewards fit, $2a$ penalizes the number of parameters, smaller is better. The sample size $n = 575$ is not needed for AIC — it only enters BIC’s penalty $a \log n$.

Model	$a$	$l_a$	$-2l_a$	$2a$	AIC	BIC
I	1	$-585$	1170	2	1172	1176.4
II	2	$-572$	1144	4	1148	1156.7
III	4	$-553$	1106	8	1114	1131.4
IV	6	$-552$	1104	12	1116	1142.1
V	8	$-550$	1100	16	1116	1150.8

Model III has the smallest AIC. Beyond it the fit improves too slowly to pay for the extra parameters: III $\to$ IV gains 1 log-likelihood unit (worth 2 on the $-2l$ scale) for a penalty of 4; III $\to$ V gains 3 units (worth 6) for a penalty of 8. BIC, with the heavier $\log(575) \approx 6.35$ per parameter, also selects III — and, as the notes say, BIC would never pick a larger model than AIC does.

Script: hw2_q10.R.

Q11 — Best model by BIC

Answer: III

\[\mathrm{BIC} = -2\,l_a + a \log(n), \qquad \log(750) = 6.62\]

(notes p. 33). Same structure as AIC but each parameter now costs $\log n = 6.62$ instead of $2$, so BIC leans harder toward small models. Smaller is better.

Model	$a$	$l_a$	$-2l_a$	$a\log 750$	BIC	AIC
I	1	$-1536$	3072	6.6	3078.6	3074
II	2	$-1531$	3062	13.2	3075.2	3066
III	4	$-1506$	3012	26.5	3038.5	3020
IV	8	$-1506$	3012	53.0	3065.0	3028
V	16	$-1499$	2998	105.9	3103.9	3030

Model III has the smallest BIC. Two things make this clear without even finishing the arithmetic:

IV vs. III — identical log-likelihood, twice the parameters, so IV can only be worse under any penalized criterion.
V vs. III — fit improves by 7 units ($14$ on the $-2l$ scale) but costs $12 \times 6.62 = 79$ extra penalty. Under AIC that same trade costs only $12 \times 2 = 24$, which is why V is much closer to III on AIC (3030 vs. 3020) than on BIC (3104 vs. 3038).

AIC also picks III, but the gap illustrates the notes’ point that BIC “will tend to result in models with fewer parameters as compared to AIC.”

Script: hw2_q11.R.

Q12 — Mean absolute error ($L_1$ loss)

Answer: 2.23

\[\text{MAE} = \frac{1}{n}\sum_{i=1}^{n} |Y_i - \hat{Y}_i|\]

(notes §1.8.1). Term by term:

$Y_i$	$\hat{Y}_i$	$\lvert Y_i - \hat{Y}_i \rvert$
14	12.5	1.5
5	9.3	4.3
9	13.8	4.8
8	6.2	1.8
14	13.1	0.9
6	5.9	0.1
	sum	13.4

\[\text{MAE} = \frac{13.4}{6} = 2.2333 \to \mathbf{2.23}\]

(For contrast, the $L_2$ loss / MSE is $\frac{1}{6}\sum (Y_i - \hat{Y}_i)^2 = 7.97$ — it is much larger because the two errors near 4–5 get squared.)

Script: hw2_q12.R.

Q13 — PPV at two thresholds

Answer: 0.20

Lower scores suggest disease, so the notes’ cognitive-test convention applies: $R \le \theta \Rightarrow T^{+}$, $R > \theta \Rightarrow T^{-}$. Then

\[\text{PPV} = P(D^{+} \mid T^{+}) = \frac{a}{a+b}\]

with $a$ = diseased and test-positive, $b$ = not diseased and test-positive. No score equals exactly 30 or 50, so $\le$ vs. $<$ is moot.

$\theta = 30$ — test-positive means score $\le 30$:

	$D^{+}$	$D^{-}$
$T^{+}$ ($\le 30$)	2, 5, 23 → $a = 3$	7 → $b = 1$
$T^{-}$ ($> 30$)	$c = 7$	$d = 7$

\[\text{PPV}_{30} = \frac{3}{3+1} = 0.7500\]

$\theta = 50$ — test-positive means score $\le 50$:

	$D^{+}$	$D^{-}$
$T^{+}$ ($\le 50$)	2, 5, 23, 36, 41, 46 → $a = 6$	7, 34, 38, 46, 46 → $b = 5$
$T^{-}$ ($> 50$)	$c = 4$	$d = 3$

\[\text{PPV}_{50} = \frac{6}{6+5} = 0.5455\] \[\text{PPV}_{30} - \text{PPV}_{50} = 0.7500 - 0.5455 = 0.2045 \to \mathbf{0.20}\]

The stricter threshold flags fewer people, so it misses more cases (sensitivity $0.30$ vs. $0.60$), but the ones it does flag are more often truly diseased — the usual sensitivity/PPV trade-off as $\theta$ moves.

Script: hw2_q13.R — builds the full $2 \times 2$ table at each threshold and prints sensitivity and specificity too.

Q14 — Correct statement about the ROC curve

Answer: Each point on the ROC curve corresponds to a (sensitivity, 1 − specificity) pair at a particular threshold.

Notes §1.8.2: the ROC curve “is constructed by plotting 1-specificity vs. sensitivity for each potential value of $\theta$.” Sweep $\theta$ from one extreme to the other and each threshold contributes one point, $\big(1 - \text{Spec}(\theta),\ \text{Sens}(\theta)\big)$, from $(0, 0)$ to $(1, 1)$.

Why the others are wrong:

“Plots PPV on the y-axis versus FPR” — the $y$-axis is sensitivity, the $x$-axis is $1 - $ specificity $=$ FPR. A curve with PPV on the $y$-axis is a precision–recall curve, not an ROC.
“Depends on disease prevalence, so the AUC changes too” — sensitivity $P(T^{+} \mid D^{+})$ and specificity $P(T^{-} \mid D^{-})$ each condition on disease status, so they are “characteristics of the test” (notes p. 36) and do not move with prevalence; neither does the curve or its AUC. It is PPV and NPV that depend on prevalence.
“AUC can decrease under an increasing rescaling (e.g. logit)” — the AUC is $P(X > Y)$, the Mann–Whitney probability that a random $D^{+}$ score exceeds a random $D^{-}$ score. A strictly increasing transform preserves every pairwise ordering, so every $I(X_i > Y_j)$, every (sens, spec) pair at the matched threshold, and the AUC are all unchanged.

Q15 — Likelihood ratio test for nested Cox models

Answer: 0.04

Model 1 (10 predictors, $l_b = -722$) is nested in Model 2 (14 predictors, $l_a = -717$). Notes p. 33:

\[2\,(l_a - l_b) \sim \chi^2_{a-b}\]

$\text{LRT} = 2\,{-717 - (-722)} = 2 \times 5 = 10$
$\text{df} = 14 - 10 = 4$ — one per added predictor
$p = P(\chi^2_4 > 10) = 0.0404 \to \mathbf{0.04}$

pchisq(10, df = 4, lower.tail = FALSE)   # 0.0404

Since $p < 0.05$, Model 2 — the larger one — is the superior model: the four extra predictors improve the fit by more than chance would explain. Common slips to avoid: forgetting the factor of 2 (giving $\chi^2 = 5$, $p = 0.29$), or using df $= 14$ instead of the difference in parameters.

Script: hw2_q15.R.