Week 4: More Likelihood Models
PS 818 - Data Analysis with Statistical Models
September 28, 2026
More Likelihood Models
\[ \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator*{\argmax}{arg\,max} \]
Previously
Generalized linear models
\[g(\mathbb{E}[Y_i|X_i]) = X_i^{\prime}\beta\]
Binary outcome models
- Careful with interpreting the raw model coefficients
- For a logit, \(\beta_k\) is a change in log-odds, and \(\frac{\partial \pi_i}{\partial X_{ik}} = \beta_k\pi_i(1-\pi_i)\)
Inference on transformations:
- Invariance + delta method
This week
- Latent variables and discrete choice
- The latent-variable representation of the logit and probit
- Ordinal logit and the proportional odds assumption
- Multinomial logit and independence of irrelevant alternatives (IIA)
- Random utility and conditional logit
- Choice data and alternative-specific covariates
- The model based approach to conjoint experiments
- Duration models
- Survival functions, the hazard, and censoring
- Parametric, semi-parametric (Cox), and non-parametric (Kaplan-Meier) approaches
“Latent variables” and logit
Latent variables
Another common way of formulating discrete outcome regressions is in terms of an unobserved continuous “latent” variable and an observation mechanism.
- Instead of putting a distribution on \(Y_i\), we put one on an unobserved \(Y_i^*\)
- \(Y_i\) is some function of \(Y_i^*\) (usually an indicator function of some sort)
For the logistic regression
\[Y_i = \mathbb{1}(Y_i^* \ge 0)\] \[Y_i^* = X_i^{\prime}\beta + \epsilon_i\] \[\epsilon_i \sim \text{Logistic}(0, 1)\]
The Logistic distribution is parameterized in terms of a “location” (mean) \(\mu\) and a “scale” parameter \(s>0\).
- For \(\mu = 0\), \(s= 1\) we have the “standard” logistic
\[\Pr(\epsilon_i \le x) = \frac{1}{1 + \exp(-x)}\]
Latent variables
So what’s \(\mathbb{E}[Y_i|X_i] = \Pr(Y_i = 1|X_i)\)?
\[\Pr(Y_i = 1|X_i) = \Pr(Y_i^* \ge 0 | X_i) = \Pr(\epsilon_i \ge -X_i^{\prime}\beta) = \Pr(\epsilon_i \le X_i^{\prime}\beta) = \frac{1}{1 + \exp(-X_i^{\prime}\beta)}\]
- The last-but-one step uses the symmetry of the logistic distribution around zero.
We can write the probit similarly, defining the error distribution as
\[\epsilon_i \sim \text{Normal}(0, 1)\]
And by extension, the probability \(\Pr(Y_i = 1|X_i)\)
\[\Pr(Y_i = 1|X_i) = \Pr(Y_i^* \ge 0|X_i) = \Pr(\epsilon_i \le X_i^{\prime}\beta) = \Phi(X_i^{\prime}\beta)\]
Note the scale of \(\epsilon_i\) is fixed by assumption, not estimated
- \(Y_i^*\) is never observed, so only the ratio \(\beta/s\) is identified – so we normalize to \(s =1\)
Latent variables
- We can visualize this as an integral over the density of the random “error” \(\epsilon_i\)
Latent variables
- Same thing for a probit
Ordinal logit
- Using the latent variable formulation lets us extend binary outcome models to outcomes with multiple discrete levels.
- For example, consider an “ordinal” response variable
- (e.g.) a Likert scale (1 = “strongly disagree”, 3 = “neutral”, 5 = “strongly agree”)
- The responses are ordered but the intervals between them do not have any meaning.
- A popular modeling technique with an ordinal \(Y\) is to treat it as an expression of a latent variable.
Ordinal logit
In an ordered logit model, an ordinal outcome \(Y\) with \(L\) unique, ordered values \(Y \in \{1, 2, 3, \dotsc L\}\) is a function of a latent variable and a set of \(L-1\) cutpoints \(\kappa\)
\[Y_i = \begin{cases}1 &\text{if } Y_i^* \le \kappa_1 \\ 2 &\text{if } \kappa_1 < Y_i^* \le \kappa_2\\ 3 &\text{if } \kappa_2 < Y_i^* \le \kappa_3\\ \vdots &\\ L &\text{if } Y_i^* > \kappa_{L-1}\end{cases}\]
We still assume \(Y_i^*\) has a logistic distribution
\[Y_i^* = X_i^{\prime}\beta + \epsilon_i, \qquad \epsilon_i \sim \text{Logistic}(0, 1)\]
Ordinal logit
To estimate the model, we now have to estimate both the \(\beta\) parameters and the \(\kappa\) parameters
- Note we omit the intercept from \(X_i\) – the \(\kappa\) parameters are essentially the “intercepts” for each boundary between the choices.
How does the model map changes in the latent variable to changes in probability? Work with the cumulative probabilities:
\[\Pr(Y_i \le l \,|\, X_i) = \Pr(Y_i^* \le \kappa_l \,|\, X_i) = \Pr(\epsilon_i \le \kappa_l - X_i^{\prime}\beta)\]
And then difference adjacent cumulatives to get each cell probability
\[\Pr(Y_i = l \,|\, X_i) = \Pr(Y_i \le l \,|\, X_i) - \Pr(Y_i \le l-1 \,|\, X_i)\]
with the conventions \(\Pr(Y_i \le 0) = 0\) and \(\Pr(Y_i \le L) = 1\).
Ordinal logit
- Let’s think in terms of the distribution of \(Y_i^*\) and the cut-points for a case with \(L=3\). For \(X_i^{\prime}\beta = 0\):
- The green lines are the cutpoints \(\kappa_1, \kappa_2\); the shaded areas are \(\Pr(Y_i = 1), \Pr(Y_i = 2), \Pr(Y_i = 3)\).
Ordinal logit
- Shifting the mean of \(Y^*\) up by 3:
- The cutpoints stay put but the distribution slides. All three probabilities move together as a result.
Ordinal logit
This structural mapping from \(Y^*\) to \(\Pr(Y_i = l)\) encodes a particular assumption about how these probabilities relate to one another.
Referred to as the proportional odds assumption
\[\log\left(\frac{\Pr(Y_i \le l)}{1 - \Pr(Y_i \le l)}\right) = \kappa_l - X_i^{\prime}\beta\]
- The same \(\beta\) appears for every \(l\); only the intercept \(\kappa_l\) changes.
Intuitively: the \(\beta\) parameters are shared across all the possible outcome levels.
- Information about whether someone responds “strongly agree” also tells us something about how likely they would be to respond “agree”, “neutral”, “disagree”, or “strongly disagree”.
- Not innocuous if there’s non-monotonicity or heterogeneous effects.
- Be careful with these models!
Multinomial logit
In a setting with multiple unordered outcomes, the natural extension of the logit regression is the multinomial logit
Suppose we have an outcome with \(K\) unordered levels (e.g. “choices”). Define a latent “utility” \(Y_{ik}^*\) for each of the \(K\) categories
- For identification, we set one category to be the “baseline” by anchoring its coefficients to zero: \(\beta_1 = 0\)
\[\begin{align*}Y_{i1}^* &= X_i^{\prime}\beta_1 + \epsilon_{i1} = \epsilon_{i1}\\ Y_{i2}^* &= X_i^{\prime}\beta_2 + \epsilon_{i2}\\ &\vdots \\ Y_{iK}^* &= X_i^{\prime}\beta_K + \epsilon_{iK}\end{align*}\]
The probability of observing choice \(Y_{i} = k\) is the probability that the \(k\)th latent variable is greater than all of the others
\[\Pr(Y_i = k \,|\, X_i) = \Pr\big(Y^*_{ik} = \max\{Y^*_{i1}, Y^*_{i2}, \dotsc Y^*_{iK}\}\big)\]
Multinomial logit
It turns out that if the error terms \(\epsilon_{ik}\) have i.i.d. “Type I extreme value” distributions, we can write the choice probabilities in closed form
\[\Pr(Y_i = k \,|\, X_i) = \frac{\exp(X_i^{\prime}\beta_{k})}{\sum_{k^\prime=1}^K \exp(X_i^{\prime}\beta_{k^\prime})}\]
This is sometimes referred to as the “softmax” function in machine learning.
Note: with the normalized baseline category \(\beta_1 = 0\), we can write this equivalently as
\[\Pr(Y_i = k \,|\, X_i) = \frac{\exp(X_i^{\prime}\beta_{k})}{1 + \sum_{k^\prime=2}^K \exp(X_i^{\prime}\beta_{k^\prime})}\]
The i.i.d. error assumption implies independence of irrelevant alternatives (IIA).
- The relative odds of choice \(k\) vs. \(l\) don’t depend on what other options are available.
- Famously implausible in some choice settings (the “red bus/blue bus” problem).
Example: Appointer Affiliation Bias
- Puig and Strezhnev (2017) ran a survey experiment on a group of international arbitration experts to see how manipulating the appointing party affected their decisionmaking in a vignette
- Respondents were given a hypothetical dispute and asked to determine who should pay the legal costs
- Three options: Loser pays all, Loser pays some, Both sides split
Example: Appointer Affiliation Bias
- Read in the responses –
-99is a missing-value code
- The outcome is stored as an integer, so recode it into an ordered factor
study1 <- study1 %>%
mutate(isReimbursementFolded_Factor = case_when(
isReimbursementFolded == 1 ~ "Winner pays all",
isReimbursementFolded == 2 ~ "Winner pays some",
isReimbursementFolded == 3 ~ "Split costs",
isReimbursementFolded == 4 ~ "Loser pays some",
isReimbursementFolded == 5 ~ "Loser pays all"))
study1$isReimbursementFolded_Factor <- factor(
study1$isReimbursementFolded_Factor,
levels = c("Winner pays all", "Winner pays some", "Split costs",
"Loser pays some", "Loser pays all"), ordered = TRUE)Example: Appointer Affiliation Bias
- Keep the completed responses, and drop the two “winner pays” categories
- Winner pays is basically an error, so we’re left with the three options in the vignette
Winner pays all Winner pays some Split costs Loser pays some
0 0 76 59
Loser pays all
109 Example: Appointer Affiliation Bias
Example: Appointer Affiliation Bias
- It would seem like an ordinal logit regression would be appropriate – a three-category discrete outcome with a clear ordering
- Implemented in the
ordinalpackage
- Implemented in the
formula: isReimbursementFolded_Factor ~ isAppointerSide
data:
study1_complete %>% filter(isAppointerSide == "Appointed by Winner" | isAppointerSide == "Appointed by Loser")
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 123 -128.28 262.57 5(0) 2.98e-13 1.5e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
isAppointerSideAppointed by Winner 0.411 0.340 1.21 0.23
Threshold coefficients:
Estimate Std. Error z value
Split costs|Loser pays some -0.672 0.262 -2.56
Loser pays some|Loser pays all 0.302 0.256 1.18Example: Appointer Affiliation Bias
- The model implies that arbitrators appointed by the winning party shift mostly from splitting the costs to making the loser pay for everything
- But is this real or an artifact of the proportional odds assumption? Is there a “constant effect” on the latent variable scale?
Example: Appointer Affiliation Bias
- Compare with multinomial logit (which in this case is equivalent to OLS for each category since we’re using a fully-saturated regression)
- Running the ordinal logit (as a reviewer asked us to do!) would have given us the wrong result
- Remember, don’t just fixate on the outcome when choosing a modeling strategy - think about what your entire bundle of assumptions implies!
Conditional logit
Same latent utility setup as the multinomial logit, but now the covariates are attribute-specific but their types are shared across alternatives
- Price of each mode of transportation, the attributes of each profile in a conjoint
\(X_{ik}\) varies across categories, but all categories share a common \(\beta\) vector
\[Y_{ik}^* = X_{ik}^{\prime}\beta + \epsilon_{ik}\]
With the same i.i.d. Type I extreme value errors, we get the same “softmax” expression for the choice probability.
\[\Pr(Y_i = k \,|\, X_{i1}, \dotsc, X_{iK}) = \frac{\exp(X_{ik}^{\prime}\beta)}{\sum_{k^\prime=1}^K \exp(X_{ik^\prime}^{\prime}\beta)}\]
IIA
We still have the independence of irrelevant alternatives assumption: the relative odds of two choices only depend on those two choices’ attributes.
\[\frac{\Pr(Y_i = j)}{\Pr(Y_i = k)} = \exp\big((X_{ij} - X_{ik})^{\prime}\beta\big)\]
Extensions involve imposing additional hierarchical structure:
- Nested logit - group similar alternatives and allow correlation within groupings
- A random effect common to every alternative in the same group \(g(k)\): \(Y_{ik}^* = X_{ik}^{\prime}\beta + \zeta_{i,g(k)} + \epsilon_{ik}\)
- Mixed logit / random coefficients – let \(\beta\) vary across individuals (Train 2009)
- Unit-specific coefficients drawn from a common distribution: \(Y_{ik}^* = X_{ik}^{\prime}\beta_i + \epsilon_{ik}\)
- Increasingly estimated via Bayesian methods (for regularization + partial pooling)
- Nested logit - group similar alternatives and allow correlation within groupings
Conjoint experiments
- A conjoint experiment shows respondents two (or more) profiles with randomly assigned attributes and asks them to pick one.
- e.g. choosing between two candidates with different backgrounds, partisan affiliation, policy positions, etc…
- Often designed as a forced choice between alternatives described by alternative-specific covariates
- In political science: analyzed as a factorial experiment
- Estimand is the AMCE (Hainmueller et al. 2014) – the average change in the probability of selection from changing one attribute, marginalizing over the others.
- Estimated by OLS, no distributional assumptions, “causal” interpretation.
- In literally every other field that does conjoints: analyzed with a random utility model
- Model latent utilities and simulate choice probabilities.
Example: Workplace democracy
- Mazumder and Yan (2024) ask whether Americans want to work at democratically governed firms
- YouGov Blue, June 2020: 1,002 respondents, four tasks each, pairwise comparison of two hypothetical employers
- Outcome: at which firm would you rather work? (forced choice)
- Each firm profile randomizes 19 attributes: the one emphasized in the paper is corporate governance
- Private, non-worker shareholders (baseline); Private, one owner; Public, shareholders; Workers are shareholders; Workers on the board; Workers elect managers
- Plus pay and benefits, hours, location, owners’ race and gender, unionization, political donations, workplace conditions, …
- The paper reports AMCEs from OLS. We’ll fit the conditional logit to the same choices.
Example: Workplace democracy
- One row per respondent-task-profile: exactly the long format the conditional logit wants
- We drop the 122 tasks where the respondent picked neither firm, leaving 3,886 tasks from 993 respondents
# A tibble: 4 × 7
rid task profile chosen governance income hours
<int> <dbl> <int> <int> <fct> <fct> <fct>
1 1 1 1 1 Workers are shareholders $200k 80
2 1 1 2 0 Private, one owner $260k 80
3 1 2 1 1 Private, one owner $90k 40
4 1 2 2 0 Workers on the board $200k 40 dfidx()tellsmlogitwhich rows form a choice situation (task_id, nested in respondentrid) and which alternative each row is (profile)
Example: Workplace democracy
- All 19 attributes are alternative-specific with a generic coefficient: the first part of the
mlogitformula| 0drops the alternative-specific constant (a preference for whichever profile is listed second)
[1] 73- Respondents complete four tasks each, so cluster by respondent
estfun()has one row per task (one likelihood contribution each), so the cluster vector needs one entry per task
Example: Workplace democracy
- The coefficients are part-worths: each level’s contribution to the systematic component of utility, relative to its attribute’s baseline
| level | part_worth | se | odds_ratio | pr_vs_baseline |
|---|---|---|---|---|
| Private, one owner | 0.061 | 0.088 | 1.06 | 0.515 |
| Public, shareholders | 0.283 | 0.088 | 1.33 | 0.570 |
| Workers are shareholders | 0.420 | 0.091 | 1.52 | 0.604 |
| Workers on the board | 0.359 | 0.090 | 1.43 | 0.589 |
| Workers elect managers | 0.192 | 0.089 | 1.21 | 0.548 |
- Codetermination (Workers on the board) beats an otherwise identical firm that is privately owned by non-worker shareholders (the baseline) about 59% of the time
Part-worths: the firm
Conditional logit with all 19 attributes; hollow points are baselines. 95% CIs clustered by respondent. Data from Mazumder and Yan (2024).
Part-worths: the job
Conditional logit with all 19 attributes; hollow points are baselines. 95% CIs clustered by respondent. Data from Mazumder and Yan (2024).
Part-worths: salary
- 28 salary levels, so each dummy is noisy. A single log-salary term seems to fit well:
Red: salary dummies with 95% CIs. Blue: log-salary fit with 95% CI, centered on the dummies’ average level. Both clustered by respondent.
Predicted choices
Conditional logit with salary as log(salary); respondent-clustered covariance. Data from Mazumder and Yan (2024). CIs: delta method on the log-odds scale.
Survival models
Survival models
Often we have outcome data that is measured as the “time-to-event”
- From the medical literature, this is often the number of days before a patient suffers some adverse event (e.g. death).
- But lots of other applications (e.g. time to treaty ratification)
Assume we observe \(n\) i.i.d. event times \(T_1, T_2, \dotsc, T_n\)
Define the “survival function”
\[S(t) = Pr(T > t) = 1 - F(t)\]
We’ll also often work with the “hazard” or the “instantaneous” probability of an event happening given that it has not yet occurred
\[h(t) = \underset{\Delta t \to 0}{\lim} \frac{P(t < T_i \le t + \Delta t | T_i > t)}{\Delta t} = \frac{f(t)}{S(t)} = \frac{f(t)}{1-F(t)}\]
Intuitively: Given that I know that an event has not occurred until \(t\), what is the probability that it will immediately occur next!
Survival models
We’ll also define the cumulative hazard
\[H(t) = \int_{0}^{t} h(u)du = \int_0^{t} \frac{f(u)}{S(u)} du = -\log(S(t))\]
Non-parametric approach
- The easiest way to estimate \(S(t) = 1 - F(t)\) is to just plug in \(\hat{F}(t)\) - our empirical CDF in the sample
- At any given value \(t\), what share of events occur after \(t\)
- Problem: This doesn’t address censoring of observations – sometimes we stop observing a unit at time \(t^*\) rather than observing an actual event.
- We know that \(T_i > t^*\) but not what \(T_i\) actually equals
- Let \(C = \{C_1, C_2, \dotsc, C_n\}\) denote the censoring indicator
- If \(C_i = 1\), then the observation is censored, the observed value \(t\) only gives us information on \(T_i > t\) as opposed to \(T_i = t\).
- Crucially, we assume that censoring is “uninformative”
- One approach to constructing a non-parametric estimator of the survival function under censoring is the Kaplan-Meier estimator.
Kaplan-Meier Estimator
Start by arranging the possible event times from \(0 < t_{(1)} < t_{(2)} < t_{(3)} < t_{(4)}, \dotsc, < t_{(n)}\)
- Suppose \(t\) is after \(t_{(j)}\) and before \(t_{(j+1)}\). We can write the survival function in terms of conditional probabilities:
\[S(t) = Pr(T_i > t | T_i > t_{(j)}) \times Pr(T_i > t_{(j)} | T_i > t_{(j-1)}) \times \dotsc \times Pr(T_i > t_{(2)} | T_i > t_{(1)}) \times Pr(T_i > t_{(1)})\]
The Kaplan-Meier estimator \(\hat{S}(t)\) is a “plug-in” estimator for each of these probabilities.
\[\hat{S}(t) = \prod_{j: t_{(j)} \le t} \bigg(1 - \frac{d_j}{r_j} \bigg)\]
\(d_j\) denotes the number of units that experience an event at time \(t_{(j)}\)
\(r_j\) is the number of units that are still in the risk set just before \(t_{(j)}\)
- Units drop out of the risk set \(r\) over time by either experiencing an event or being censored
Kaplan-Meier Estimator
Intuition: Consider one of the terms in the product. Using the definition of conditional probability:
\[Pr(T_i > t_{(j)} | T_i > t_{(j-1)}) = 1 - Pr(T_i \le t_{(j)}|T_i > t_{(j-1)}) = 1 - \frac{Pr(T_i \le t_{(j)}, T_i > t_{(j-1)})}{Pr(T_i > t_{(j-1)})}\]
In the numerator, our “plug-in” estimator is just the share of units that have events occurring at \(t_{(j)}\).
In the denominator we have the share of units that we know for sure are still “active” in the interval after \(t_{(j-1)}\)
Hence our overall “plug-in” estimator:
\[\widehat{Pr}(T_i > t_{(j)} | T_i > t_{(j-1)}) = 1 - \frac{d_j}{r_j}\]
Taking the product across all of the \(j\) where \(t_{(j)} \le t\) yields the Kaplan-Meier estimator.
Illustration: Lung Cancer
- Let’s look at a sample dataset from the
survivalR package on survival times for patients with advanced lung cancer.
# A tibble: 228 × 10
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
7 7 310 2 68 2 2 70 60 384 10
8 11 361 2 71 2 2 60 80 538 1
9 1 218 2 53 1 1 70 80 825 16
10 7 166 2 61 1 2 70 70 271 34
# ℹ 218 more rows- We want to estimate \(S(t | \text{Sex} = \text{Male})\) and \(S(t | \text{Sex} = \text{Female})\)
Illustration: Lung Cancer
- In the
survivalpackage, we need to define a “survival” outcome in terms of the event time and the censoring indicator
Illustration: Lung Cancer
- Plot the curve for the sample as a whole
Illustration: Lung Cancer
- Now let’s do this for gender
Illustration: Lung Cancer
- Plot separate curves for men and women
Parametric inference
When we want to model the survival function using covariates, we run into a clear constraint using purely non-parametric methods
- Splitting the sample and constructing separate K-M curves for each unique covariate profile quickly becomes infeasible.
We may want to assume a model for the event time. We can write that density in terms of the hazard and the survival function
\[f(t) = h(t)S(t)\]
Note that by the definition of the cumulative hazard \(H(t) = - \log(S(t))\) and \(S(t) = \exp(-H(t))\), so
\[f(t) = h(t)\exp(-H(t))\]
Parametric inference
Depending on the assumptions we make about the nature of the hazard, we can obtain different distributions for \(f(t)\).
- In the simplest case, let’s assume that the hazard does not vary over time: \(h(t) = \lambda\), with \(\lambda > 0\)
- And the cumulative hazard is \(H(t) = \int_{0}^{t} \lambda du = \lambda t\)
This yields the exponential distribution
\[f(t) = \lambda \exp(-\lambda t)\]
Exponential Model
The mean of \(T_i \sim \text{Exponential}(\lambda)\) is \(E[T_i] = \frac{1}{\lambda}\). We’ll often reparameterize the mean as \(\theta\), yielding a density of
\[f(t) = \frac{1}{\theta} \exp\bigg(-\frac{1}{\theta} t\bigg)\]
\(\theta\) is strictly greater than \(0\), so we’ll use a log-link to incorporate the covariates:
\[\log(E[T_i | X_i]) = \log(\theta_i) = X_i^{\prime}\beta\]
Alternatively, we can write
\[E[T_i | X_i] = \exp(X_i^{\prime}\beta)\]
Note that the constant hazard assumption is strong
- The exponential distribution is “memoryless” in that knowing that a unit has survived up until time \(t\) tells you nothing about its instantaneous probability of failure.
- More flexible distributions like the Weibull allow for a time-varying hazard.
Exponential Model
- We can interpret changes in the linear predictor in terms of changes in the hazard since
\[h(t) = \frac{1}{\theta_i} = \frac{1}{\exp(X_i^{\prime}\beta)}\]
Increases in the linear predictor reduce the hazard.
This also lets us write down the survival function
\[S(t) = \exp\bigg(- \frac{1}{\exp(X_i^{\prime}\beta)}t \bigg)\]
Application
- Let’s compare our parametric fit to the non-parametric Kaplan-Meier estimates
- Get the survival curves for men and women
Application
- Overlay the exponential survival curves on the Kaplan-Meier estimates
Application: Weibull
The fit seems to be off.
- This is because the hazard appears to be time-varying.
- You can see a kind of “s-curve” which the exponential cannot capture.
The Weibull adds a shape parameter \(k\) that lets the hazard change over time
\[h(t) = \frac{k}{\lambda}\left(\frac{t}{\lambda}\right)^{k-1}\qquad S(t) = \exp\left(-\left(\frac{t}{\lambda}\right)^{k}\right)\]
- \(k = 1\) is the exponential; \(k > 1\) is a rising hazard, \(k < 1\) a falling one
Application: Weibull
- Get the Weibull survival curves (\(\lambda_i = \exp(X_i^{\prime}\beta)\)) and overlay them on Kaplan-Meier
Cox Proportional Hazards
- The Cox proportional hazards model assumes that covariates enter the hazard log-linearly, but the baseline hazard is unspecified
\[h_i(t) = h_0(t)\exp(X_i^{\prime}\beta)\]
\(h_0(t)\) is the unknown baseline hazard function.
Note that while the baseline hazard is not explicitly identified, the ratio of the hazards of two observations \(i\) and \(j\) at time \(t\) does not depend on the baseline
\[\frac{h_i(t)}{h_j(t)} = \frac{\exp(X_i^{\prime}\beta)}{\exp(X_j^{\prime}\beta)}\]
The model is semi-parametric in that some components are left entirely unspecified…
- …but others are not – we assume that the covariates act proportionally on the hazard irrespective of time.
Cox Proportional Hazards
The model is estimated by maximizing a partial likelihood
- We observe event time \(Y_i\) - assume no ties
The individual partial likelihood is
\[L_i(\beta) = \frac{\exp(X_i^{\prime}\beta)}{\sum_{k; Y_k \ge Y_i} \exp(X_k^{\prime}\beta)}\]
Intuitively: The denominator contains the hazards of all observations in the risk set at the event time \(Y_i\)
Then the full likelihood is just the product over all the individual partial likelihoods
\[L(\beta) = \prod_{i: C_i = 0} \frac{\exp(X_i^{\prime}\beta)}{\sum_{k; Y_k \ge Y_i} \exp(X_k^{\prime}\beta)}\]
Cox and conditional logit
The individual partial likelihood should look familiar - it has exactly the form of the conditional logit
\[L_i(\beta) = \frac{\exp(X_i^{\prime}\beta)}{\sum_{k; Y_k \ge Y_i} \exp(X_k^{\prime}\beta)}\]
At each event time we are asking: given that one unit in the risk set failed, which one was it?
- The “choice set” is the risk set \(\{k: Y_k \ge Y_i\}\), and the “chosen alternative” is the unit that failed.
Application
- What does the Cox fit look like
Call:
coxph(formula = Surv(time, status) ~ I(sex == 2), data = lung)
coef exp(coef) se(coef) z p
I(sex == 2)TRUE -0.5 0.6 0.2 -3 0.001
Likelihood ratio test=11 on 1 df, p=0.001
n= 228, number of events= 165 - Note that we don’t get an estimate of the hazard because the baseline is left unspecified
- If we want to plot survival curves, we’ll need some estimate of \(h_0\) using another method (like Kaplan-Meier)
Application
- Estimate the baseline cumulative hazard non-parametrically (the Breslow estimator) and get the Cox-predicted survival curves for men and women
Next week
- Introduction to Bayesian inference
- What happens if we assume parameters are not fixed but carry uncertainty?
- The posterior distribution \(p(\theta | \mathbf{Y})\) as a target of inference
- The prior distribution and how it shapes the posterior.
- Midterm Exam
References
PS 818 - University of Wisconsin-Madison