The number of events

table(df$cut_hi)

## 
##    0    1 
## 1336   40

We see that there were 40 challengers entering the race. In other words, it is a fairly rare event.

We can chech whether there are several events per incumbent candidate.

tab <- 
df %>%
  group_by(caseid) %>%
  summarize(challengers = sum(cut_hi))

tab$challengers %>% table

## .
##   0   1 
## 357  40

There isn’t. The study stops when an incumbent is challenged, meaning this is not a repeated event study.

A combined outcome variable: segments for durations

We are interested in the combiation of these two outcomes. How can we explore that? One way is to plot the durations and flag the ones that finishes in a challenge.

The base-line code is this.

df0 %>% 
ggplot +
  #Plot a segment per incumbent candidate
  geom_segment(aes(
    #The segment starts at week 0
    x = 0,
    #... and ends at the end of the study
    xend = te,
    #y axis reports the candidate id; here I order them according to their duration
    y = caseid %>% as.factor,
    yend = caseid %>% as.factor,
    color = cut_hi %>% as.factor),
    position = "identity") 

We see that most candidates last (“survive”) for a long time, and that the number of challenges is small relative to the number of candidates. The plot is too crammed, however, so let’s do some partitionning and esthetical tweaks.

I split the data into Democrats and Republicans. Here, I only do the democrats.

df0 %>% 
  #There are so many candidates, that I split them by party
  filter(dem == 1) %>%
  #I give the challenger variable some more human values
  mutate(challenger = if_else(cut_hi == 1, 
                              "Challenger enters",
                              "No challenge")) %>%
ggplot +
  #Plot a segment per incumbent candidate
  geom_segment(aes(
    #The segment starts at week 0
    x = 0,
    #... and ends at the end of the study
    xend = te,
    #y axis reports the candidate id; here I order them according to their duration
    y = reorder(caseid, te),
    yend = reorder(caseid, te),
    color = challenger %>% as.factor),
    position = "identity", 
    #Transparent segments with a size
    alpha = 0.5,
    size = 0.8) +
  #I parse out the events in a different pane
  facet_grid(cols = vars(challenger)) +
  scale_color_manual(values = c("blue", "black"),
                     name = "Challenger enters") +
  #Smaller text on the y axis
  theme(axis.text = element_text(size = 5),
        #supress legend
        legend.position = "none") +
  ggtitle("Campaign duration and challenger entry",
          subtitle = "Among DEMOCRATS") +
  #no y-axis label
  ylab("") +
  xlab("Campaign duration (weeks)")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

A combined variable for censored data

We are now moving into the domain of event history functions. Let’s explore how R combines the two outcomes when the estimation is done. Note that I am now working in the main data frame again (df)

The Surv function in the survival package creates a new, combined dependent variable for us. (If you havent installed the package yet, you will have to do it before you fetch it in the library).

We report the duration with start and stop times. (time= and time2=) and whether or not an event has taken place at the end of this period (event=).

In this data, there is no starting time. The te variable only reports the duration at the end of the period for the observation, so we need to create a starting point. This starting point is either 0 (the first risk set) or the end-point of the previous risk set. In other words, we need to do some recoding where we calculate the lag of the end-time.

df <- 
  df %>%
  #Sort data so that first event is first etc
  arrange(caseid, te) %>%
  #within each candidate
  group_by(caseid) %>%
  #find the previous duration at the end of the period
  mutate(start = lag(te))

#First period starts at 0
df$start[is.na(df$start)] = 0

Now, we are good to go. R can now label the survival time/duration with whether or not this duration is censored (+) or not.

Surv(time = df$start, time2 = df$te, event = df$cut_hi) %>% head

## [1] ( 0,26+] (26,50]  ( 0,26+] (26,53+] (53,65+] (65,75+]

This will constitute the dependent variable in all of our analyses. It’s a combination of three pieces of information.

We can explore the base-line hazard/survival in a model with no covariates by creating a Kaplan-Meyer curve. It is somewhat of a homologue to the geom_smooth() statistics in ggplot2. It estimates hazard rates and survival curves for us from the empirical distribution. That is, it estimates the survival rate at each event over time. It gives us a sneak peak of the shape of the base-line hazard. Of course, our predictors may change that shape, but it gives us an indication of the distribution, just as histogram would do. Incidentially, it also happens to function the same way as the Cox proportional hazard model.

KaplanMeier <- survfit(Surv(start, te, cut_hi) ~ 1, data = df)
KaplanMeier

## Call: survfit(formula = Surv(start, te, cut_hi) ~ 1, data = df)
## 
##      records   n events median 0.95LCL 0.95UCL
## [1,]    1376 397     40     NA      NA      NA

The r object tells us how many observations we have and how many events. The object is a list containing a lot of information. It reports the cumulative hazard (the increasing hazard of being challenged), the time and the survival rates.

Among other things, it reports the risk-sets. Risk sets are the set of people/observations that are “alive” at time \(t\), running the risk of being challenged (i.e. experiencing the event). These are the “mini-datasets” that R uses when it calculates its “mini-logits”. For each risk set, a model is estimated. The \(\betas\) are then aggregated/averaged over to report a single effect of \(x\) with a single slope coefficient (\(\beta\)). In other words, it is all the unique duration periods where an event happened or – in the end of the study – did not happen. We could have calculated that from our df0.

KaplanMeier$n.risk %>% length

## [1] 80

df0$te%>% unique %>% length

## [1] 80

We can plot the descriptive statistics in many ways. We can do it by hand simply by extracting the information ourselves. The robject contains all of these variables: upper (conf.int), lower, time, surv (mean predicted survival rate).

  ggplot() +
  geom_ribbon(aes(ymin = KaplanMeier$lower,
                  ymax = KaplanMeier$upper,
                  x = KaplanMeier$time),
              alpha = 0.2) +
  geom_line(aes(x = KaplanMeier$time,
                y = KaplanMeier$surv)) +
  ggtitle("Kaplan-Meyer: a base-line survival curve")

There are ready-mades that can do this for you in no-time. Here, in baseR

KaplanMeier %>%
  plot

Or in ggplot, using the ggsurvfit package.

library(ggsurvfit)

## Warning: package 'ggsurvfit' was built under R version 4.2.3

ggsurvfit(KaplanMeier) +
  add_confidence_interval() 

It is clear that most incumbent candidates stand the race all the way through to election. We see that from the survival probability at the end of the study.

Segments

Using the segments code from the previous section, I could split or color the data according to their values on the main predictors. We’re moving towards bivariate exploration.

For fun, I recode the warchest into rich/not rich incumbents.

df0$ec %>%hist

df0$rich <- ifelse(df0$ec > 0.25, "rich", "not rich")

… and color code accordingly. The code is mostly a copy-paste from prevously.

The only new thing here, is the color specification and me moving the legend to the bottom of the plot.

df0 %>% 
  #There are so many candidates, that I split them by party
  filter(dem == 1) %>%
  #I give the challenger variable some more human values
  mutate(challenger = if_else(cut_hi == 1, 
                              "Challenger enters",
                              "No challenge")) %>%
ggplot +
  #Plot a segment per incumbent candidate
  geom_segment(aes(
    #The segment starts at week 0
    x = 0,
    #... and ends at the end of the study
    xend = te,
    #y axis reports the candidate id; here I order them according to their duration
    y = reorder(caseid, te),
    yend = reorder(caseid, te),
    color = rich),
    position = "identity", 
    #Transparent segments with a size
    alpha = 0.5,
    size = 0.8) +
  #I parse out the events in a different pane
  facet_grid(cols = vars(challenger)) +
  scale_color_manual(values = c("blue", "black"),
                     name = "Challenger enters") +
  #Smaller text on the y axis
  theme(axis.text = element_text(size = 5),
        #move legend
        legend.position = "bottom") +
  ggtitle("Campaign duration and challenger entry",
          subtitle = "Among DEMOCRATS") +
  #no y-axis label
  ylab("") +
  xlab("Campaign duration (weeks)")

Now, we can clearly see that he density of “rich” candidates is higher in the “No challenge” group, and high up in the plot where the long durations are. In other words the idea that rich candidates deter newcomers is visible both in the duration (somewhat) and the event/no event data.

Your turn

Can you do this for republicans? And give the bars a continous color range according to the size of their warchest?

Kaplan Meyer

We can also partition the data into groups and create different Kaplan-Meyer curves for subgroups of the data. This is only useful for categorical predictors, so I recode my data again. Note that I am now working in the main data frame (df).

#Create a bin
df$rich <- ifelse(df$ec > 0.25, "rich", "not rich")

KaplanMeier <- survfit(Surv(start, te, cut_hi) ~ rich, data = df)

#Plot using the ggsurvfit package
ggsurvfit(KaplanMeier) +
  add_confidence_interval() +
  ggtitle("Base-line survival rate using ggsurvfit")

Here, we see that the “not rich” candidates have a higher probability of being challenged, and that these challenges seem to arrive after a year (50-75 weeks). The richest candidates have a really high probability of survival (no challenge) throughout the campaign, although – there too – there’s a bit of rush in challenges towards the end.

Cox proportional hazard

The Cox proportional hazard model is semi-parametric in the sense that it simply uses the empirical base-line hazard (survival rate; number of observations still “alive” and those that experience the event at time \(t\)), once our predictors have accounted for their variation.

mod.cox <- coxph(Surv(start, te, cut_hi) ~ 
                   ec
                 + cluster(caseid),
                 df)

summary(mod.cox)

## Call:
## coxph(formula = Surv(start, te, cut_hi) ~ ec, data = df, cluster = caseid)
## 
##   n= 1376, number of events= 40 
## 
##        coef exp(coef) se(coef) robust se      z Pr(>|z|)   
## ec -3.68704   0.02505  1.35286   1.38564 -2.661  0.00779 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##    exp(coef) exp(-coef) lower .95 upper .95
## ec   0.02505      39.93  0.001657    0.3786
## 
## Concordance= 0.665  (se = 0.046 )
## Likelihood ratio test= 10.8  on 1 df,   p=0.001
## Wald test            = 7.08  on 1 df,   p=0.008
## Score (logrank) test = 6.57  on 1 df,   p=0.01,   Robust = 12.71  p=4e-04
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

The model output gives us quite some information. It reports the raw coefficient, i.e. the log hazard ratio, but also the hazard ratios, the standard error, z-stats and p-value. In many studies, the results are reported as hazard ratios (the exponentiated coefficient), but the standard error is not transformed. This is extremely useless: How can I compare the size of the coefficient and its estimated variation (s.e.)? Usually, it is better to report either the non-transformed or both in combination with the standard error. Alternatively, you might report the confidence interval.

1-exp(mod.cox$coefficients[1]*0.1)

##        ec 
## 0.3083697

Clustering and frailties?

When the data is split into a panel, as is the case here, each candidate is observed several times. However, in event history models we don’t necessarily account for the additional observations. This is obvious in single-event studies with continious time where there is only one duration (and line in the data) per candidate. But even when the data is split into a panel but with only one event, as long as the candidate is present only once per risk-set, we don’t have to account for unobserved heterogeneity.

If had been in the presence of a multiple-event study where each candidate can go in and out of spells (being challenged several times), we would account for it by adding a “frailty” term. This is equivalent to a random-intercept model. We would have done that by the function frailty(caseid). Here, this would overidentify the model, so it is not possible. Another possibility in all models is to account for correlations in the observations by clustering the standard errors, but leave the regression coefficients intact. We can do this directly with the cluster function. I could, for example, cluster per candidate: cluster(caseid).

In this study, I do neither.

Weibull model

The Weibull survival model is a commonly used parametric model in survival analysis. Just like the Cox model, it allows for the estimation of the hazard function of a survival time. It assumes that the hazard function follows a Weibull distribution. It has two parameters: a scale parameter (often denoted as \(\lambda\)) and a shape parameter (often denoted as \(k\)).

The Weibull distribution has a monotonic hazard function, meaning that the hazard either increases or decreases over time, depending on the value of the shape parameter.

mod.weibull <- phreg(Surv(start, te, cut_hi) ~ 
                        ec
                      + iv
                      + dem
                      + south,
                      dist = "weibull",
                        df) 

These parameters are reported in the output.

The exponential model

The exponential model is just a sub-category of the Weibull model. While the weibull usually estimates its two parameters in view of the data, we can also force the shape parameter to 1. The result is an assumption that the hazard rate is constant over time. This means that the probability of an event occurring in a given time interval is the same regardless of how long the candidate has been in the study. The exponential model is often used as a baseline model for comparison with more complex models.

mod.exponential <- phreg(Surv(start, te, cut_hi) ~ 
                        ec
                      + iv
                      + dem
                      + south,
                      dist = "weibull",
                      shape = 1,
                        df) 

stargazer(mod.full.cox, mod.weibull, mod.exponential,
          # mod.full.curv.cox, mod.full.tt.cox,
          colnames = T,
          title = "Survival models",
          type = "html")

How do I choose the right model? If I have a clear theory about the shape of the hazard function, then I have a fairly easy time choosing a model (based on theory). Often this is not the case, the variation in hazard over the duration of the study is often noise/a disturbance in relation to our research question. This is often why people opt for the Cox model. However, if the cox model underperforms (i.e. it is overidentified, has many ties etc.) or we are looking to make predictions, we might opt for the parametric models (remember, the cox model is close to the data in terms of the base-line hazard).

The choice between different parametric models often ends up being a comparison of which one performs the best. How are their predictions.

Lastly, robust research results are not modl-dependent. If we find support for our claims in different models, that’s an asset, while if we only find support in one type of model, we need to ask ourselves why.

Fixing the problem

If the proportional hazards assumption is not met for one or more covariates in the Cox model, there are several strategies that can be used to address the violation. These include stratifying on the covariate, including interaction terms between the covariate and time, or using time-varying coefficients.

Strata

Stratification in the event history models functions as fixed-effects. We can ask the model to make comparisons strictly within each stratum. The difference is that the model output doesn’t report the effect of the fixed-effect; it just makes the comparisons without much ado.

The other difference is more conceptual: To fix the time-problem, we basically create dummies (strata) for different slices of the duration/spell. That is, we manually ask the model to change the level of the base-line hazard for slices of the data.

The following is a dirty hack insofar as it is completely informed by the data (i.e. I am looking to describe the data appropriately, assuming it is noise for my theory), but it can work miracles. Look at the Schoenfield residuals plot again. It seems that something seems to happen around 75 weeks into the campaign. We saw that already in the descriptive statistics thanks to the Kaplan Meier curves. The challenging of candidates picks up speed/the challenger rate increases. Here, we can see a slight bend upwards for both ev and ec.

To do this, we also need to split the data into new appropriate period/observations according to this new definition. The survSplit() function allows me to make many manual turning points and reshape the data accordingly. I store the output in a new dataframe.

df1 <- survSplit(Surv(start, te, cut_hi)~.,
          data = df,
          cut = c(73),
          #The name of the new "dummy"
          episode = "stratum")

#compare
dim(df);dim(df1)

## [1] 1376    9

## [1] 1519   10

I now have a new dataset with new observations, and we can estimate our model anew, but tell R that the baseline hazard is different towards the end; i.e. in our two strata (before and after week 73).

mod.full.strata.cox <- coxph(Surv(start, te, cut_hi) ~ 
                        #warchest
                        ec
                      #prior vote share
                      + iv
                      #democrats
                      + dem
                      #from southern US states
                      + south
                      #Stratum
                      + strata(stratum)
                      ,
                      ties = "breslow",
                      tt = function(x,t, ...){
                        x*log(t)},
                        id = caseid,
                      x = T,
                      na.action = "na.omit",
                        df1)

summary(mod.full.strata.cox)

## Call:
## coxph(formula = Surv(start, te, cut_hi) ~ ec + iv + dem + south + 
##     strata(stratum), data = df1, na.action = "na.omit", ties = "breslow", 
##     x = T, tt = function(x, t, ...) {
##         x * log(t)
##     }, id = caseid)
## 
##   n= 1519, number of events= 40 
## 
##             coef  exp(coef)   se(coef)      z Pr(>|z|)    
## ec    -3.0088853  0.0493467  1.3864091 -2.170    0.030 *  
## iv    -6.9736554  0.0009362  1.6574337 -4.208 2.58e-05 ***
## dem    0.2340996  1.2637703  0.3239421  0.723    0.470    
## south -0.4380054  0.6453223  0.4239028 -1.033    0.301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## ec    0.0493467    20.2648 3.259e-03   0.74709
## iv    0.0009362  1068.1200 3.636e-05   0.02411
## dem   1.2637703     0.7913 6.698e-01   2.38456
## south 0.6453223     1.5496 2.812e-01   1.48116
## 
## Concordance= 0.776  (se = 0.042 )
## Likelihood ratio test= 34.11  on 4 df,   p=7e-07
## Wald test            = 25.01  on 4 df,   p=5e-05
## Score (logrank) test = 26.32  on 4 df,   p=3e-05

None of the coefficients changed; now remember we didn’t really have a problem to begin with, so to us, this is good news. The Wald test tells us that the explanatory power of our new model is better. That’s unsurprising, since we just added some fixed-effects into the mix. The Log likelihood is the same, while \(R^2\) is lower.

Linear fixes

The stratification provides manual and descrete fixes to the problem that the baseline hazard and some of the covariates correlate by changing the (ability of) the base-line hazard (to vary). We could do this with linear predictors instead. That is, we do as we’d do in an OLS: We recode the predictors so that they are no longer linear. I can do this by creating an interaction between x and time or an interaction between x and x for example (i.e. a curvilinear effect). The argument tt= allows me to define an R-function that is applied to all the variables of my choice in the model. I apply the function by wrapping them into tt().

Take 1: Here, I recode the relationship the candidate’s past vote share and itself (iv) such that I get a curvlinear effect.

mod.full.curv.cox <- coxph(Surv(start, te, cut_hi) ~ 
                           #warchest
                           ec
                         #prior vote share
                         + iv
                         + tt(iv)
                         #democrats
                         + dem
                         #from southern US states
                         + south
                         ,
                         ties = "breslow",
                         tt = function(x, ...){x*x},
                         id = caseid,
                         df)

summary(mod.full.curv.cox)

## Call:
## coxph(formula = Surv(start, te, cut_hi) ~ ec + iv + tt(iv) + 
##     dem + south, data = df, ties = "breslow", tt = function(x, 
##     ...) {
##     x * x
## }, id = caseid)
## 
##   n= 1376, number of events= 40 
## 
##              coef  exp(coef)   se(coef)      z Pr(>|z|)  
## ec     -2.906e+00  5.471e-02  1.380e+00 -2.105   0.0353 *
## iv     -2.645e+01  3.262e-12  1.192e+01 -2.219   0.0265 *
## tt(iv)  1.377e+01  9.529e+05  8.224e+00  1.674   0.0941 .
## dem     1.817e-01  1.199e+00  3.263e-01  0.557   0.5777  
## south  -5.170e-01  5.963e-01  4.270e-01 -1.211   0.2260  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## ec     5.471e-02  1.828e+01 3.657e-03 8.185e-01
## iv     3.262e-12  3.066e+11 2.343e-22 4.540e-02
## tt(iv) 9.529e+05  1.049e-06 9.525e-02 9.534e+12
## dem    1.199e+00  8.339e-01 6.326e-01 2.273e+00
## south  5.963e-01  1.677e+00 2.582e-01 1.377e+00
## 
## Concordance= 0.773  (se = 0.041 )
## Likelihood ratio test= 36.6  on 5 df,   p=7e-07
## Wald test            = 32.78  on 5 df,   p=4e-06
## Score (logrank) test = 40.56  on 5 df,   p=1e-07

Take 2: Here, I create an interaction term. That is I recode the relationship between time and the candidate’s past vote share (iv) such that the effect of x is proportional/exponential with time.

mod.full.tt.cox <- coxph(Surv(start, te, cut_hi) ~ 
                           #warchest
                           ec
                         #prior vote share
                         + tt(iv)
                         #democrats
                         + dem
                         #from southern US states
                         + south
                         ,
                         ties = "breslow",
                         tt = function(x,t, ...){
                           x*log(t)},
                         id = caseid,
                         df)

summary(mod.full.tt.cox)

## Call:
## coxph(formula = Surv(start, te, cut_hi) ~ ec + tt(iv) + dem + 
##     south, data = df, ties = "breslow", tt = function(x, t, ...) {
##     x * log(t)
## }, id = caseid)
## 
##   n= 1376, number of events= 40 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## ec     -3.12756   0.04382  1.38792 -2.253   0.0242 *  
## tt(iv) -1.67940   0.18649  0.41627 -4.034 5.47e-05 ***
## dem     0.22727   1.25517  0.32392  0.702   0.4829    
## south  -0.42788   0.65189  0.42360 -1.010   0.3124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## ec       0.04382    22.8182  0.002886    0.6654
## tt(iv)   0.18649     5.3623  0.082474    0.4217
## dem      1.25517     0.7967  0.665246    2.3682
## south    0.65189     1.5340  0.284188    1.4954
## 
## Concordance= 0.771  (se = 0.044 )
## Likelihood ratio test= 32.05  on 4 df,   p=2e-06
## Wald test            = 23.74  on 4 df,   p=9e-05
## Score (logrank) test = 24.71  on 4 df,   p=6e-05

This latter solution is the classic way of fixing violations of the proprtional hazard assumption. You can define the shape of the relationship yourself; often by inferring it from the Shoenfield residuals plot.

stargazer(mod.full.cox, mod.full.strata.cox,
          mod.full.curv.cox, mod.full.tt.cox,
          colnames = T,
          title = "Different fixes to violations of the proportional hazard assumption",
          type = "html")

For each attempt at fixing the problem, I would obviously go back and check whether the proportional hazard assumption is now respected. Unfortunately, the cox.zph() function doesn’t handle tt definititions.