Event History Analysis - Example 5 Cox Proportional Hazards Model Part 2 - Model Checking
coreysparks
February 18, 2015
This example will illustrate how to fit parametric the Cox Proportional hazards model to a discrete-time (longitudinal) data set and examine various model diagnostics to evaluate the overall model fit. The data example uses data from the ECLS-K. Specifically, we will examine the transition into poverty between kindergarten and third grade.
#Load required libraries
library(foreign)
library(survival)## Loading required package: splineslibrary(car)
library(survey)## Loading required package: grid
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchartlibrary(muhaz)Using Longitudinal Data
As in the other examples, I illustrate fitting these models to data that are longitudinal, instead of person-duration. In this example, we will examine how to fit the Cox model to a longitudinally collected data set.
First we load our data
load("~/Google Drive/dem7903_App_Hier/data/eclsk.Rdata")
names(eclsk)<-tolower(names(eclsk))
#get out only the variables I'm going to use for this example
myvars<-c( "childid","gender", "race", "r1_kage","r4age", "r5age", "r6age", "r7age","c1r4mtsc", "c4r4mtsc", "c5r4mtsc", "c6r4mtsc", "c7r4mtsc", "w1povrty","w1povrty","w3povrty", "w5povrty", "w8povrty","wkmomed", "s2_id", "c1_5fp0", "c15fpstr", "c15fppsu")
eclsk<-eclsk[,myvars]
eclsk$age1<-ifelse(eclsk$r1_kage==-9, NA, eclsk$r1_kage/12)
eclsk$age2<-ifelse(eclsk$r4age==-9, NA, eclsk$r4age/12)
#for the later waves, the NCES group the ages into ranges of months, so 1= <105 months, 2=105 to 108 months. So, I fix the age at the midpoint of the interval they give, and make it into years by dividing by 12
eclsk$age3<-recode(eclsk$r5age,recodes="1=105; 2=107; 3=109; 4=112; 5=115; 6=117; -9=NA")/12
eclsk$pov1<-ifelse(eclsk$w1povrty==1,1,0)
eclsk$pov2<-ifelse(eclsk$w3povrty==1,1,0)
eclsk$pov3<-ifelse(eclsk$w5povrty==1,1,0)
#Recode race with white, non Hispanic as reference using dummy vars
eclsk$hisp<-recode (eclsk$race, recodes="3:4=1;-9=NA; else=0")
eclsk$black<-recode (eclsk$race, recodes="2=1;-9=NA; else=0")
eclsk$asian<-recode (eclsk$race, recodes="5=1;-9=NA; else=0")
eclsk$nahn<-recode (eclsk$race, recodes="6:7=1;-9=NA; else=0")
eclsk$other<-recode (eclsk$race, recodes="8=1;-9=NA; else=0")
eclsk$race_gr<-recode(eclsk$race, recodes="3:4='hisp'; 2='nh black'; 5='nh asian'; 6:7='nahn'; 8='other'; 1='nh white'; else=NA", as.factor.result = T)
eclsk$race_gr<-relevel(eclsk$race_gr, ref = 'nh white')
eclsk$male<-recode(eclsk$gender, recodes="1=1; 2=0; -9=NA")
eclsk$mlths<-recode(eclsk$wkmomed, recodes = "1:2=1; 3:9=0; else = NA")
eclsk$mgths<-recode(eclsk$wkmomed, recodes = "1:3=0; 4:9=1; else =NA") Now, I need to form the transition variable, this is my event variable, and in this case it will be 1 if a child enters poverty between the first wave of the data and the third grade wave, and 0 otherwise. NOTE I need to remove any children who are already in poverty age wave 1, because they are not at risk of experiencing this particular transition.
eclsk<-subset(eclsk, is.na(pov1)==F&is.na(pov2)==F&is.na(pov3)==F&is.na(age1)==F&is.na(age2)==F&is.na(age3)==F&pov1!=1&is.na(eclsk$c15fpstr)==F)
eclsk$povtran1<-ifelse(eclsk$pov1==0&eclsk$pov2==0, 0,1)
eclsk$povtran2<-ifelse(eclsk$povtran1==1, NA,ifelse(eclsk$pov2==0&eclsk$pov3==0,0,1))Now we do the entire data set. To analyze data longitudinally, we need to reshape the data from the current “wide” format (repeated measures in columns) to a “long” format (repeated observations in rows). The reshape() function allows us to do this easily. It allows us to specify our repeated measures, time varying covariates as well as time-constant covariates.
e.long<-reshape(eclsk, idvar="childid", varying=list(age=c("age1","age2"), age2=c("age2", "age3"), povtran=c("povtran1", "povtran2")), times=1:2, direction="long" , drop = names(eclsk)[4:20])
e.long<-e.long[order(e.long$childid, e.long$time),]
#find which kids failed in the first time period and remove them from the second risk period risk set
failed1<-which(is.na(e.long$povtran1)==T)
e.long<-e.long[-failed1,]
e.long$age1r<-round(e.long$age1, 0)
e.long$age2r<-round(e.long$age2, 0)
head(e.long, n=10)## childid gender race c1_5fp0 c15fpstr c15fppsu pov1 pov2 pov3
## 0001002C.1 0001002C 2 1 352.35 63 1 0 0 0
## 0001002C.2 0001002C 2 1 352.35 63 1 0 0 0
## 0001007C.1 0001007C 1 1 358.99 63 1 0 0 0
## 0001007C.2 0001007C 1 1 358.99 63 1 0 0 0
## 0001010C.1 0001010C 2 1 352.35 63 1 0 0 0
## 0001010C.2 0001010C 2 1 352.35 63 1 0 0 0
## 0002004C.1 0002004C 2 1 228.05 63 1 0 0 0
## 0002004C.2 0002004C 2 1 228.05 63 1 0 0 0
## 0002006C.1 0002006C 2 1 287.25 63 1 0 1 1
## 0002008C.1 0002008C 2 1 259.19 63 1 0 0 0
## hisp black asian nahn other race_gr male mlths mgths time
## 0001002C.1 0 0 0 0 0 nh white 0 0 0 1
## 0001002C.2 0 0 0 0 0 nh white 0 0 0 2
## 0001007C.1 0 0 0 0 0 nh white 1 0 1 1
## 0001007C.2 0 0 0 0 0 nh white 1 0 1 2
## 0001010C.1 0 0 0 0 0 nh white 0 0 1 1
## 0001010C.2 0 0 0 0 0 nh white 0 0 1 2
## 0002004C.1 0 0 0 0 0 nh white 0 0 1 1
## 0002004C.2 0 0 0 0 0 nh white 0 0 1 2
## 0002006C.1 0 0 0 0 0 nh white 0 NA NA 1
## 0002008C.1 0 0 0 0 0 nh white 0 0 1 1
## age1 age2 povtran1 age1r age2r
## 0001002C.1 6.433333 7.894167 0 6 8
## 0001002C.2 7.894167 9.750000 0 8 10
## 0001007C.1 5.300000 6.825000 0 5 7
## 0001007C.2 6.825000 8.916667 0 7 9
## 0001010C.1 5.885833 7.385833 0 6 7
## 0001010C.2 7.385833 9.333333 0 7 9
## 0002004C.1 5.450000 6.977500 0 5 7
## 0002004C.2 6.977500 9.083333 0 7 9
## 0002006C.1 5.158333 6.685833 1 5 7
## 0002008C.1 5.905833 7.433333 0 6 7Now we fit the Cox model using full survey design. In the ECLS-K, I use the longitudinal weight for waves 1-5, as well as the associated psu and strata id’s for the longitudinal data from these waves from the parents of the child, since no data from the child themselves are used in the outcome.
des2<-svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights=~c1_5fp0, data=e.long[complete.cases(e.long),], nest=T)
#Fit the model
fitl1<-svycoxph(Surv(time = time, event = povtran1)~mlths+mgths+race_gr, design=des2)
summary(fitl1) ## Stratified 1 - level Cluster Sampling design (with replacement)
## With (501) clusters.
## svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights = ~c1_5fp0,
## data = e.long[complete.cases(e.long), ], nest = T)## Call:
## svycoxph(formula = Surv(time = time, event = povtran1) ~ mlths +
## mgths + race_gr, design = des2)
##
## n= 12980, number of events= 617
##
## coef exp(coef) se(coef) z Pr(>|z|)
## mlths 0.8410 2.3186 0.1301 6.466 1.00e-10 ***
## mgths -1.0004 0.3677 0.1267 -7.898 2.78e-15 ***
## race_grhisp 0.8431 2.3236 0.1207 6.986 2.82e-12 ***
## race_grnahn 1.1961 3.3073 0.3511 3.407 0.000658 ***
## race_grnh asian 0.2013 1.2229 0.3302 0.610 0.542138
## race_grnh black 1.3878 4.0062 0.1640 8.463 < 2e-16 ***
## race_grother 0.9351 2.5474 0.3226 2.899 0.003744 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## mlths 2.3186 0.4313 1.7969 2.9918
## mgths 0.3677 2.7194 0.2869 0.4713
## race_grhisp 2.3236 0.4304 1.8342 2.9437
## race_grnahn 3.3073 0.3024 1.6618 6.5818
## race_grnh asian 1.2229 0.8177 0.6403 2.3358
## race_grnh black 4.0062 0.2496 2.9049 5.5248
## race_grother 2.5474 0.3926 1.3537 4.7935
##
## Concordance= 0.757 (se = 0.01 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 7 df, p=NA
## Wald test = 452.1 on 7 df, p=0
## Score (logrank) test = NA on 7 df, p=NAModel Residuals
First, we will extract the Shoenfeld residuals, which are useful for examining non-proportional hazards with repect to time. This means that the covariate effect could exhibit time-dependency. First we extract the residuals from the model, then we fit a linear model to the residual and the observed (uncensored) failure times, we DO NOT WANT TO SEE A SIGNIFICANT MODEL!!!!!
schoenresid<-resid(fitl1, type="schoenfeld")
fit.sr<-lm(schoenresid~des2$variables$time[des2$variables$povtran==1])
summary(fit.sr)## Response mlths :
##
## Call:
## lm(formula = mlths ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2353 -0.2317 -0.2111 -0.2075 0.7925
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.051550 0.050684
## des2$variables$time[des2$variables$povtran == 1] 0.003608 0.034392
## t value Pr(>|t|)
## (Intercept) -1.017 0.310
## des2$variables$time[des2$variables$povtran == 1] 0.105 0.916
##
## Residual standard error: 0.4168 on 615 degrees of freedom
## Multiple R-squared: 1.789e-05, Adjusted R-squared: -0.001608
## F-statistic: 0.011 on 1 and 615 DF, p-value: 0.9165
##
##
## Response mgths :
##
## Call:
## lm(formula = mgths ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3890 -0.3677 -0.3571 0.6217 0.6429
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 0.03779 0.05877
## des2$variables$time[des2$variables$povtran == 1] 0.01065 0.03988
## t value Pr(>|t|)
## (Intercept) 0.643 0.52
## des2$variables$time[des2$variables$povtran == 1] 0.267 0.79
##
## Residual standard error: 0.4833 on 615 degrees of freedom
## Multiple R-squared: 0.000116, Adjusted R-squared: -0.00151
## F-statistic: 0.07132 on 1 and 615 DF, p-value: 0.7895
##
##
## Response race_grhisp :
##
## Call:
## lm(formula = race_grhisp ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3158 -0.3059 -0.2720 0.6842 0.7379
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.08854 0.05493
## des2$variables$time[des2$variables$povtran == 1] 0.04380 0.03727
## t value Pr(>|t|)
## (Intercept) -1.612 0.108
## des2$variables$time[des2$variables$povtran == 1] 1.175 0.240
##
## Residual standard error: 0.4517 on 615 degrees of freedom
## Multiple R-squared: 0.00224, Adjusted R-squared: 0.0006173
## F-statistic: 1.38 on 1 and 615 DF, p-value: 0.2405
##
##
## Response race_grnahn :
##
## Call:
## lm(formula = race_grnahn ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06177 -0.06177 -0.06002 -0.05060 0.95114
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 0.03237 0.02817
## des2$variables$time[des2$variables$povtran == 1] -0.01117 0.01912
## t value Pr(>|t|)
## (Intercept) 1.149 0.251
## des2$variables$time[des2$variables$povtran == 1] -0.584 0.559
##
## Residual standard error: 0.2317 on 615 degrees of freedom
## Multiple R-squared: 0.0005543, Adjusted R-squared: -0.001071
## F-statistic: 0.3411 on 1 and 615 DF, p-value: 0.5594
##
##
## Response race_grnh asian :
##
## Call:
## lm(formula = `race_grnh asian` ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.05472 -0.05326 -0.05149 -0.05003 0.94997
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 0.019694 0.027014
## des2$variables$time[des2$variables$povtran == 1] 0.003232 0.018331
## t value Pr(>|t|)
## (Intercept) 0.729 0.466
## des2$variables$time[des2$variables$povtran == 1] 0.176 0.860
##
## Residual standard error: 0.2221 on 615 degrees of freedom
## Multiple R-squared: 5.056e-05, Adjusted R-squared: -0.001575
## F-statistic: 0.0311 on 1 and 615 DF, p-value: 0.8601
##
##
## Response race_grnh black :
##
## Call:
## lm(formula = `race_grnh black` ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1827 -0.1735 -0.1547 -0.1456 0.8544
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.08127 0.04484
## des2$variables$time[des2$variables$povtran == 1] 0.02795 0.03043
## t value Pr(>|t|)
## (Intercept) -1.813 0.0704 .
## des2$variables$time[des2$variables$povtran == 1] 0.918 0.3587
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3687 on 615 degrees of freedom
## Multiple R-squared: 0.00137, Adjusted R-squared: -0.000254
## F-statistic: 0.8436 on 1 and 615 DF, p-value: 0.3587
##
##
## Response race_grother :
##
## Call:
## lm(formula = race_grother ~ des2$variables$time[des2$variables$povtran ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.02757 -0.02757 -0.02473 -0.02207 0.98077
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -0.001473 0.018765
## des2$variables$time[des2$variables$povtran == 1] -0.005504 0.012733
## t value Pr(>|t|)
## (Intercept) -0.078 0.937
## des2$variables$time[des2$variables$povtran == 1] -0.432 0.666
##
## Residual standard error: 0.1543 on 615 degrees of freedom
## Multiple R-squared: 0.0003037, Adjusted R-squared: -0.001322
## F-statistic: 0.1868 on 1 and 615 DF, p-value: 0.6657And we jump up and down! None of our predictors appear to be associated with the residuals, which suggests the effects are constant over time
Not so soon
We can also get a formal test using weighted residuals in a nice pre-rolled form with a plot, a la Grambsch and Therneau (1994) :
fit.test<-cox.zph(fitl1)
fit.test## rho chisq p
## mlths 0.00056 3.32e-04 0.985473
## mgths 0.03938 2.71e+00 0.099900
## race_grhisp 0.03464 1.46e+00 0.227101
## race_grnahn 0.06499 1.13e+01 0.000757
## race_grnh asian -0.05807 7.22e+00 0.007198
## race_grnh black -0.04165 2.56e+00 0.109776
## race_grother 0.07559 7.13e+00 0.007594
## GLOBAL NA 2.41e+01 0.001100par(mfrow=c(3,3))
plot(fit.test, df=2)
par(mfrow=c(1,1))
Here, we see that the nahn, asian and other race/ethnic groups have significant tests in the formal test, which suggests a non-proportionality according to race/ethnicity.
Next we examine Martingale residuals. Martingale residuals are also useful for assessing the functional form of a covariate. A plot of the martingale residuals against the values of the covariate will indicate if the covariate is being modeled correctly, i.e. linearly in the Cox model. If a line fit to these residuals is a straight line, then the covariate has been modeled effectively, if it is curvilinear, you may need to enter the covariate as a quadratic, although this is not commonly a problem for dummy variables.
#extract Martingale residuals
res.mar<-resid(fitl1, type="martingale")
#plot vs maternal education
scatter.smooth(des2$variables$mlths, res.mar,degree = 2,span = 1, ylab="Martingale Residual",col=1, cex=.25, lpars=list(col = "red", lwd = 3))
title(main="Martingale residuals for Mother < High School")
See, this doesn’t really make sense, why would a dummy variable have a quadratic effect? The scatter plot smoother lies!
Stratification
Above, we observed evidence of non-proportional effects by race/ethnicity. There are a few standard ways of dealing with this in practice. The first is stratification of the model by the offending predictor. If one of the covariates exhibits non-proportionality we can re-specify the model so that each group will have its own baseline hazard rate. This is direct enough to do by using the strata() function within a model. This is of best use when a covariate is categorical, and not of direct importance for our model (i.e. a control variable).
fitl2<-svycoxph(Surv(time = time, event = povtran1)~mlths+mgths+strata(race_gr), design=des2, subset=race!=-8)
summary(fitl2) ## Stratified 1 - level Cluster Sampling design (with replacement)
## With (501) clusters.
## svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights = ~c1_5fp0,
## data = e.long[complete.cases(e.long), ], nest = T)## Call:
## svycoxph(formula = Surv(time = time, event = povtran1) ~ mlths +
## mgths + strata(race_gr), design = des2, subset = race !=
## -8)
##
## n= 12980, number of events= 617
##
## coef exp(coef) se(coef) z Pr(>|z|)
## mlths 0.8434 2.3242 0.1307 6.452 1.11e-10 ***
## mgths -1.0000 0.3679 0.1267 -7.893 3.00e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## mlths 2.3242 0.4303 1.799 3.0030
## mgths 0.3679 2.7183 0.287 0.4716
##
## Concordance= 0.688 (se = 0.016 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 2 df, p=NA
## Wald test = 159.6 on 2 df, p=0
## Score (logrank) test = NA on 2 df, p=NA#plot the two models
survfit(fitl2)## Call: survfit(formula = fitl2)
##
## records n.max n.start events median 0.95LCL 0.95UCL
## race_gr=nh white 9223 0.6781 0.6781 0.02397 NA NA NA
## race_gr=hisp 1667 0.1589 0.1589 0.01972 NA 1 NA
## race_gr=nahn 291 0.0186 0.0186 0.00251 NA 1 NA
## race_gr=nh asian 651 0.0336 0.0336 0.00174 NA 1 NA
## race_gr=nh black 799 0.0892 0.0892 0.01289 NA 1 NA
## race_gr=other 349 0.0216 0.0216 0.00211 NA 1 NApar(mfrow=c(1,2))
plot(survfit(fitl1, newdata = data.frame(mlths=0, mgths=0, race_gr=levels(e.long$race_gr))), col = 1:6, ylim=c(.7,1))
legend("bottomleft", legend=c("nh white", "hisp", "nanh", "nh asian", "nh black", "other"), col=1:6, lty=1)
title(main="Un-stratified Cox Model")
plot(survfit(fitl2,newdata = data.frame(mlths=0, mgths=0)), col = 1:6, ylim=c(.7,1))
legend("bottomleft", legend=c("nh white", "hisp", "nanh", "nh asian", "nh black", "other"), col=1:6, lty=1)
title(main="Stratified Cox Model")
Non-proportional effects with time
We can also include a timecovariate interaction term to model directly any time-dependence in the covariate effect. Different people say to do different things, some advocate for simply interacting time covariate, others say use a nonlinear function of time, e.g. log(time) * covariate, others say use time-1 * covariate. Below, I simply use time * covariate.
fit3<-svycoxph(Surv(time = time, event = povtran1)~mlths+mgths+time*race_gr, design=des2, subset=race!=-8)## Warning in fitter(X, Y, strats, offset, init, control, weights = weights,
## : Loglik converged before variable 3 ; beta may be infinite.summary(fit3)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (501) clusters.
## svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights = ~c1_5fp0,
## data = e.long[complete.cases(e.long), ], nest = T)## Call:
## svycoxph(formula = Surv(time = time, event = povtran1) ~ mlths +
## mgths + time * race_gr, design = des2, subset = race != -8)
##
## n= 12980, number of events= 617
##
## coef exp(coef) se(coef) z Pr(>|z|)
## mlths 8.127e-01 2.254e+00 1.264e-01 6.429 1.28e-10
## mgths -9.872e-01 3.726e-01 1.245e-01 -7.926 2.22e-15
## time -1.948e+01 3.469e-09 1.826e-01 -106.703 < 2e-16
## race_grhisp 7.952e-01 2.215e+00 3.862e-01 2.059 0.039496
## race_grnahn 2.846e-01 1.329e+00 1.107e+00 0.257 0.797223
## race_grnh asian 1.594e+00 4.925e+00 6.556e-01 2.432 0.015027
## race_grnh black 1.358e+00 3.889e+00 3.647e-01 3.725 0.000196
## race_grother -2.786e-01 7.568e-01 1.307e+00 -0.213 0.831131
## time:race_grhisp 2.447e-02 1.025e+00 2.623e-01 0.093 0.925676
## time:race_grnahn 6.432e-01 1.903e+00 6.047e-01 1.064 0.287468
## time:race_grnh asian -1.118e+00 3.268e-01 4.650e-01 -2.405 0.016172
## time:race_grnh black -2.410e-04 9.998e-01 2.766e-01 -0.001 0.999305
## time:race_grother 8.382e-01 2.312e+00 7.667e-01 1.093 0.274322
##
## mlths ***
## mgths ***
## time ***
## race_grhisp *
## race_grnahn
## race_grnh asian *
## race_grnh black ***
## race_grother
## time:race_grhisp
## time:race_grnahn
## time:race_grnh asian *
## time:race_grnh black
## time:race_grother
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## mlths 2.254e+00 4.436e-01 1.759e+00 2.888e+00
## mgths 3.726e-01 2.684e+00 2.919e-01 4.756e-01
## time 3.469e-09 2.883e+08 2.426e-09 4.961e-09
## race_grhisp 2.215e+00 4.515e-01 1.039e+00 4.722e+00
## race_grnahn 1.329e+00 7.523e-01 1.517e-01 1.165e+01
## race_grnh asian 4.925e+00 2.031e-01 1.362e+00 1.780e+01
## race_grnh black 3.889e+00 2.571e-01 1.903e+00 7.948e+00
## race_grother 7.568e-01 1.321e+00 5.847e-02 9.797e+00
## time:race_grhisp 1.025e+00 9.758e-01 6.128e-01 1.714e+00
## time:race_grnahn 1.903e+00 5.256e-01 5.816e-01 6.224e+00
## time:race_grnh asian 3.268e-01 3.060e+00 1.314e-01 8.131e-01
## time:race_grnh black 9.998e-01 1.000e+00 5.814e-01 1.719e+00
## time:race_grother 2.312e+00 4.325e-01 5.145e-01 1.039e+01
##
## Concordance= 0.851 (se = 0.011 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 13 df, p=NA
## Wald test = 71636 on 13 df, p=0
## Score (logrank) test = NA on 13 df, p=NA