SAS and R

Shared frailty models are a way of allowing correlated observations into proportional hazards models. Briefly, instead of l_i(t) = l_0(t)e^(x_iB), we allow l_ij(t) = l_0(t)e^(x_ijB + g_i), where observations j are in clusters i, g_i is typically normal with mean 0, and g_i is uncorrelated with g_i'. The nomenclature frailty comes from examining the logs of the g_i and rewriting the model as l_ij(t) = l_0(t)u_i*e^(xB) where the u_i are now lognormal with median 0. Observations j within cluster i share the frailty u_i, and fail faster (are frailer) than average if u_i> 1.

In SAS 9.2, this model could not be fit, though it is included in the survival package in R. (Section 4.3.2) With SAS 9.3, it can now be fit. We explore here through simulation, extending the approach shown in example 7.30.

SAS
To include frailties in the model, we loop across the clusters to first generate the frailties, then insert the loop from example 7.30, which now represents the observations within cluster, adding the frailty to the survival time model. There's no need to adjust the censoring time.

data simfrail;
 beta1 = 2;
 beta2 = -1;
 lambdat = 0.002; *baseline hazard;
 lambdac = 0.004; *censoring hazard;
 do i = 1 to 250; *new frailty loop;
 frailty = normal(1999) * sqrt(.5);
 do j = 1 to 5; *original loop;
 x1 = normal(0);
 x2 = normal(0);
 * new model of event time, with frailty added;
 linpred = exp(-beta1*x1 - beta2*x2 + frailty);
 t = rand("WEIBULL", 1, lambdaT * linpred);
 * time of event;
 c = rand("WEIBULL", 1, lambdaC);
 * time of censoring;
 time = min(t, c); * which came first?;
 censored = (c lt t);
 output;
 end;
 end;
run;

For comparison's sake, we replicate the naive model assuming independence:

proc phreg data=simfrail;
 model time*censored(1) = x1 x2;
run;

 Parameter Standard Hazard
 Parameter DF Estimate Error Chi-Square Pr> ChiSq Ratio

 x1 1 1.68211 0.05859 824.1463 <.0001 5.377
 x2 1 -0.88585 0.04388 407.4942 <.0001 0.412

The parameter estimates are rather biased. In contrast, here is the correct frailty model.

proc phreg data=simfrail;
 class i;
 model time*censored(1) = x1 x2;
 random i / noclprint;
run;
 Cov REML Standard
 Parm Estimate Error

 i 0.5329 0.07995

 Parameter Standard Hazard
 Parameter DF Estimate Error Chi-Square Pr> ChiSq Ratio

 x1 1 2.03324 0.06965 852.2179 <.0001 7.639
 x2 1 -1.00966 0.05071 396.4935 <.0001 0.364

This returns estimates gratifyingly close to the truth. The syntax of the random statement is fairly straightforward-- the noclprint option prevents printing all the values of i. The clustering variable must be specified in the class statement. The output shows the estimated variance of the g_i.

R
In our book (section 4.16.14) we show an example of fitting the uncorrelated data model, but we don't display a frailty model. Here, we use the data generated in SAS, so we omit the data simulation in R. As in SAS, it would be a trivial extension of the code presented in example 7.30. For parallelism, we show the results of ignoring the correlation, first.

> library(survival)
> with(simfrail, coxph(formula = Surv(time, 1-censored) ~ x1 + x2))

 coef exp(coef) se(coef) z p
x1 1.682 5.378 0.0586 28.7 0
x2 -0.886 0.412 0.0439 -20.2 0

with identical results to above. Note that the Surv function expects an indicator of the event, vs. SAS expecting a censoring indicator.

As with SAS, the syntax for incorporating the frailty is simple.

> with(simfrail, coxph(formula = Surv(time, 1-censored) ~ x1 + x2 
 + frailty(i)))

 coef se(coef) se2 Chisq DF p
x1 2.02 0.0692 0.0662 850 1 0
x2 -1.00 0.0506 0.0484 393 1 0
frailty(i) 332 141 0

Variance of random effect= 0.436

Here, the results differ slightly from the SAS model, but still return parameter estimates that are quite similar. We're not familiar enough with the computational methods to diagnose the differences.

In addition to the non-parametric tools discussed in recent entries, it's common to use proportional hazards regression, (section 4.3.1) also called Cox regression, in evaluating survival data.

It's important in such models to test the proportionality assumption. Below, we demonstrate doing this for a simple model from the HELP data, available at the book web site. Our results below draw on this web page, part of the excellent resources provided by the UCLA ATS site.

SAS

First, we run a proportional hazards regression to assess the effects of treatment on the time to linkage with primary care. (Data were read in and observations with missing values removed in example 7.40.) Only a portion of the results are shown.

proc phreg data=h2; 
 model dayslink*linkstatus(0)=treat;
 output out= propcheck ressch = schres;
run; 

 Parameter Standard
 Parameter DF Estimate Error Chi-Square Pr> ChiSq

 treat 1 1.59103 0.19142 69.0837 <.0001

Being in the intervention group would appear to increase the probability of being linked to primary care. But do the model assumptions hold? The output statement above makes a new data set that contains the Schoenfeld residuals. (Schoenfeld D. Residuals for the proportional hazards regresssion model. Biometrika, 1982, 69(1):239-241.) One assessment of proportional hazards is based on these residuals, which ought to show no association with time if proportionality holds. We can plot them against the time to linkage using proc sgplot (section 5.1.1) and adding a loess curve to assess the relationship.

proc sgplot data=propcheck; 
 loess x = dayslink y = schres / clm;
run;

From the resulting plot, shown above, there is an indication of a possible problem. One way to assess this is to include a time-varying covariate, an interaction between the suspect predictor(s) and the event time. This can be done easily within proc phreg. If the interaction term is significant, the null hypothesis of proportionality has been rejected.

proc phreg data=h2; 
 model dayslink*linkstatus(0)=treat treat_time; 
 treat_time = treat*log(dayslink);
run; 

 Parameter Standard
 Parameter DF Estimate Error Chi-Square Pr> ChiSq

 treat 1 4.13105 0.97873 17.8153 <.0001
 treat_time 1 -0.61846 0.22569 7.5093 0.0061

The proportional hazards assumption doesn't hold in this case.

R

In R, the time-varying covariate approach is harder to implement. However, the survival library includes a formal test based on the Schoenfeld residuals. (See Therneau and Grambsch, 2000, pages 127-142.) This is accessed by applying the cox.zph() function to the output of the coxph() function. The smallds object used below was created in the previous example. (Some output omitted.)

> library(survival)
> propcheck = coxph(Surv(dayslink, linkstatus) ~ treat, method="breslow")

> summary(propcheck) 

 coef exp(coef) se(coef) z Pr(>|z|) 
treat 1.5910 4.9088 0.1914 8.312 <2e-16 ***

> cox.zph(propcheck) 
 rho chisq p
treat -0.233 8.47 0.00361

This test agrees with the SAS approach that the assumptions that proportionality holds can be rejected. To understand why, you can plot() the test. The transform='identity" option is used to make results more similar to those obtained with SAS.

plot(cox.zph(propcheck, transform='identity'))