Example 9.7: New stuff in SAS 9.3-- Frailty models

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.