Week 4: More Likelihood Models

Previously

• Generalized linear models: a link function \(g()\) relating the CEF to the linear predictor

  \[g(\mathbb{E}[Y_i|X_i]) = X_i^{\prime}\beta\]

• Binary outcome models: the logit and probit, estimated by maximum likelihood

  • Interpretation is harder than OLS – \(\beta_k\) is a change in log-odds, and \(\frac{\partial \pi_i}{\partial X_{ik}} = \beta_k\pi_i(1-\pi_i)\) depends on every covariate.

• Inference on transformations: invariance gives us the point estimate, the delta method gives us the standard error.

• We used all of this to estimate propensity scores and construct an IPTW estimator.

This week

• Latent variables and discrete choice (finishing up from last week)
  • 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, alternative-specific covariates, and the RUM foundation
  • Conjoint experiments and model-based approaches to analyzing them
• Duration models
  • The survival function, 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_{l=1}^K \exp(X_i^{\prime}\beta_{l})}\]

• 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_{l=2}^K \exp(X_i^{\prime}\beta_{l})}\]

• 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 – -99 is the survey platform’s missing-value code

study1 <- read_csv("data/experiment1_cleaned.csv", na = c("", "NA", "-99"))

• The outcome is stored as an integer, so recode it into an ordered factor
  • The ordering is what makes this an ordinal outcome rather than a nominal one

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
  • Essentially no one chose them, so we’re left with the three options in the vignette

study1_complete <- study1 %>%
  filter(!is.na(isReimbursementFolded) & Finished == 1) %>%
  filter(isReimbursementFolded_Factor != "Winner pays some" &
         isReimbursementFolded_Factor != "Winner pays all")

table(study1_complete$isReimbursementFolded_Factor)

 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 here – a three-category discrete outcome with a clear ordering
  • Implemented in the ordinal package

ordered_logit <- ordinal::clm(isReimbursementFolded_Factor ~ isAppointerSide, data=study1_complete %>% 
                       filter(isAppointerSide == "Appointed by Winner"|isAppointerSide == "Appointed by Loser"))

summary(ordered_logit)

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.18

Example: 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?

           Group Estimate Std. Error     z Pr(>|z|)   S   2.5 % 97.5 %
 Loser pays all    0.1022     0.0839  1.22    0.223 2.2 -0.0623 0.2667
 Loser pays some  -0.0170     0.0156 -1.09    0.275 1.9 -0.0476 0.0135
 Split costs      -0.0851     0.0706 -1.21    0.228 2.1 -0.2235 0.0532

Term: isAppointerSide
Type: prob
Comparison: Appointed by Winner - Appointed by Loser

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)

multinom_logit <- nnet::multinom(isReimbursementFolded_Factor ~ isAppointerSide, data=study1_complete %>% 
                       filter(isAppointerSide == "Appointed by Winner"|isAppointerSide == "Appointed by Loser"))

marginaleffects::avg_comparisons(multinom_logit)

           Group Estimate Std. Error      z Pr(>|z|)   S    2.5 %  97.5 %
 Loser pays all     0.180     0.0892  2.015   0.0439 4.5  0.00494  0.3546
 Loser pays some   -0.197     0.0773 -2.544   0.0110 6.5 -0.34834 -0.0452
 Split costs        0.017     0.0826  0.205   0.8372 0.3 -0.14498  0.1789

Term: isAppointerSide
Type: probs
Comparison: Appointed by Winner - Appointed by Loser

• 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!

Random utility and conditional logit

From latent variables to utility

• The latent-variable story generalizes naturally: let \(U_{ij}\) be the utility individual \(i\) gets from alternative \(j\)

  \[U_{ij} = V_{ij} + \epsilon_{ij} = X_{ij}^{\prime}\beta + \epsilon_{ij}\]

• Individual \(i\) chooses the alternative with the highest utility

  \[Y_i = \argmax_{j \in \{1, \dotsc, J\}} U_{ij}\]

• The systematic component \(V_{ij}\) is what we model; \(\epsilon_{ij}\) is everything we don’t observe.

• This is the random utility model (RUM) – the microeconomic foundation for discrete choice (McFadden 1974).

The conditional logit

• McFadden (1974) showed that if the \(\epsilon_{ij}\) are i.i.d. Type-I Extreme Value (Gumbel), the choice probabilities have a closed form

  \[\Pr(Y_i = j \,|\, X_i) = \frac{\exp(X_{ij}^{\prime}\beta)}{\sum_{k=1}^J \exp(X_{ik}^{\prime}\beta)}\]

• Compare this to the multinomial logit we just derived:

  • Multinomial logit: covariates vary over individuals \((X_i)\), and we estimate a separate \(\beta_j\) per alternative.
  • Conditional logit: covariates vary over alternatives \((X_{ij})\), and we estimate a single \(\beta\).

• Conditional logit lets us model the attributes of the choices themselves – price, distance, ideological position – not just the attributes of the chooser.

Multinomial vs. conditional logit

• The practical difference is the shape of the data
  • Multinomial logit wants one row per respondent
  • Conditional logit wants one row per respondent-alternative pair
• The two can be combined – alternative-specific covariates and individual-specific covariates in one model (sometimes called the “mixed” or “universal” logit).

IIA, revisited

• The i.i.d. extreme value assumption buys us the closed form, but it also imposes IIA

  \[\frac{\Pr(Y_i = j)}{\Pr(Y_i = k)} = \exp\big((X_{ij} - X_{ik})^{\prime}\beta\big)\]

  • The relative odds of \(j\) vs. \(k\) depend only on \(j\) and \(k\) – adding or removing other alternatives changes nothing.

• The red bus/blue bus problem again, now with a name for the culprit: the errors are independent across alternatives.

• Ways out:

  • Nested logit – group similar alternatives and allow correlation within nests
  • Mixed logit / random coefficients – let \(\beta\) vary across individuals (Train 2009)
  • Multinomial probit – a full covariance matrix for \(\epsilon\), at real computational cost

Conjoint experiments

• A conjoint experiment shows respondents two (or more) profiles with randomly assigned attributes and asks them to pick one.
  • Immigrants differing in education, country of origin, job; candidates differing in party, gender, experience.
• This is a forced choice between alternatives described by alternative-specific covariates – exactly the conditional logit data structure.
• Randomization of the attributes is what makes the design so attractive: the attributes are independent of potential outcomes by construction.

Design-based vs. model-based

• The standard estimand is the AMCE (Hainmueller, Hopkins, and Yamamoto 2014) – the average change in the probability of selection from changing one attribute, marginalizing over the others.
  • Estimated by OLS of the choice indicator on attribute dummies. No distributional assumptions at all.
• The RUM gives us a model-based alternative: a conditional logit over the profiles in each task.
  • Buys us a utility interpretation and lets us model interactions and heterogeneity more cleanly.
  • Costs us the extreme value assumption on \(\epsilon_{ij}\) – and IIA comes along with it.
• Worth asking every time: what does the extra structure buy, and what would go wrong if it’s false?

Example: (TBD)

# AMCE, design-based
# lm_robust(chosen ~ educ + origin + job, data = conjoint, clusters = respondent_id)

# conditional logit, model-based
# survival::clogit(chosen ~ educ + origin + job + strata(task_id), data = conjoint)

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\), how many 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)}\) can be understood as 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 up to time \(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 < 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 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)} < t\) yields the Kaplan-Meier estimator.

Illustration: Lung Cancer

• Let’s look at a sample dataset from the survival R 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 survival package, 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 – with continuous covariates, we may only have one observation that has a unique covariate value

  • Splitting the sample and constructing separate K-M curves 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, lets assume that the hazard does not vary over time: \(h(t) = \lambda\). \(\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[Y_i | X_i]) = \log(\theta_i) = X_i^{\prime}\beta\]

• Alternatively, we can write

  \[E[Y_i | X_i] = \exp(X_i^{\prime}\beta)\]

• Note that the constant hazard assumption is a strong one – the exponential distribution is “memoryless” in that knowing that a unit has survived up until time \(t\) tells you nothing about its instantenous 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

Semi-parametric

• You’ll notice that the fit is okay, but there is some discrepancy
  • This is because the hazard appears to be time-varying.
  • You can see a kind of “s-curve” which the exponential by construction cannot capture.
• Other parametric models allow for more flexibility in the hazard rate’s evolution over time
  • For example, the Weibull model, which adds an additional shape parameter to the hazard to allow it to vary over time actually fits the data quite well!
• But suppose that we want to make very minimal assumptions about the structure of the hazard function – can we still estimate the relationship between covariates and the hazard?

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

• Look again at the individual partial likelihood – it has exactly the form of the conditional logit from earlier today

  \[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.

• This isn’t an analogy – it’s the same likelihood. survival::clogit() rewrites its call and hands it to coxph().

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

• Calculate the baseline hazard via Kaplan-Meier and overlay the Cox predicted hazards for men and women

Next week

• Bayesian inference
  • Priors, posteriors, and how this differs from what we’ve been doing
  • Conjugate models and the Beta-Binomial
  • Credible intervals vs. confidence intervals
• Readings: TBD