The PHREG Procedure

The Penalized Partial Likelihood Approach for Fitting Frailty Models

Let bold-italic gamma equals left-parenthesis gamma 1 comma ellipsis comma gamma Subscript s Baseline right-parenthesis prime be the vector of random components for the s clusters.

Gamma Frailty Model

Assume each normal e Superscript gamma Super Subscript i has an independent and identically distributed gamma distribution with mean 1 and a common unknown variance theta; that is, normal e Superscript gamma Super Subscript i is iid upper G left-parenthesis StartFraction 1 Over theta EndFraction comma StartFraction 1 Over theta EndFraction right-parenthesis. The penalty function is

minus StartFraction 1 Over theta EndFraction sigma-summation Underscript i equals 1 Overscript s Endscripts left-parenthesis gamma Subscript i Baseline minus normal e Superscript gamma Super Subscript i Superscript Baseline right-parenthesis

plus a function of theta. The penalized partial log likelihood is given by

l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis equals l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis plus StartFraction 1 Over theta EndFraction sigma-summation Underscript i equals 1 Overscript s Endscripts left-parenthesis gamma Subscript i Baseline minus normal e Superscript gamma Super Subscript i Superscript Baseline right-parenthesis

where l Subscript partial Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis is the log of any of the partial likelihoods in the sections Partial Likelihood Function for the Cox Model and The Multiplicative Hazards Model.

The profile marginal log-likelihood of this shared frailty model (Therneau and Grambsch 2000, pp. 257–258) is

l Subscript m Baseline left-parenthesis theta right-parenthesis equals l Subscript p Baseline left-parenthesis ModifyingAbove bold-italic beta With caret left-parenthesis theta right-parenthesis comma ModifyingAbove bold-italic gamma With caret left-parenthesis theta right-parenthesis comma theta right-parenthesis plus sigma-summation Underscript i equals 1 Overscript s Endscripts StartSet theta Superscript negative 1 Baseline minus left-parenthesis theta Superscript negative 1 Baseline plus d Subscript i Baseline right-parenthesis log left-bracket theta Superscript negative 1 Baseline plus d Subscript i Baseline right-bracket plus theta Superscript negative 1 Baseline log left-parenthesis theta Superscript negative 1 Baseline right-parenthesis plus log left-bracket StartFraction normal upper Gamma left-parenthesis theta Superscript negative 1 Baseline plus d Subscript i Baseline right-parenthesis Over normal upper Gamma left-parenthesis theta Superscript negative 1 Baseline right-parenthesis EndFraction right-bracket EndSet

where d Subscript i is the number of events in the ith cluster.

The maximization of this approximate likelihood is a doubly iterative process that alternates between the following two steps:

  • For a provisional value of theta, the best linear unbiased predictors (BLUP) of bold-italic beta and bold-italic gamma are computed by maximizing the penalized partial log-likelihood l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis. The marginal likelihood is evaluated. This constitutes the inner loop.

  • A new value of theta is obtained by the golden section search based on the marginal likelihood of all the previous iterations. This constitutes the outer loop.

The outer loop is iterated until the bracketing interval of theta is small.

Lognormal Frailty Model

With each gamma Subscript i having a zero-mean normal distribution and a common variance theta, the penalty function is

StartFraction 1 Over 2 theta EndFraction bold-italic gamma prime bold-italic gamma

plus a function of theta. The penalized partial log likelihood is given by

l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis equals l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis minus StartFraction 1 Over 2 theta EndFraction bold-italic gamma prime bold-italic gamma

where l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis is the log of any of the partial likelihoods in the sections Partial Likelihood Function for the Cox Model and The Multiplicative Hazards Model.

For a given theta, let bold upper H be the negative Hessian of the penalized partial log likelihood l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis; that is,

bold upper H equals bold upper H left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis equals Start 2 By 2 Matrix 1st Row 1st Column bold upper H 11 2nd Column bold upper H 12 2nd Row 1st Column bold upper H 21 2nd Column bold upper H 22 EndMatrix

where bold upper H 11 equals minus StartFraction partial-differential squared l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis Over partial-differential bold-italic beta squared EndFraction comma bold upper H 12 equals bold upper H prime 21 equals minus StartFraction partial-differential squared l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis Over partial-differential bold-italic beta partial-differential bold-italic gamma EndFraction, and bold upper H 22 equals minus StartFraction partial-differential squared l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis Over partial-differential bold-italic gamma squared EndFraction.

The marginal log likelihood of this shared frailty model is

l Subscript m Baseline left-parenthesis bold-italic beta comma theta right-parenthesis equals minus one-half log left-parenthesis theta Superscript s Baseline right-parenthesis plus log left-bracket integral normal e Superscript l Super Subscript p Superscript left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis Baseline d bold-italic gamma right-bracket

Using a Laplace approximation to the integral as in Breslow and Clayton (1993), an approximate marginal log likelihood (Ripatti and Palmgren 2000; Therneau and Grambsch 2000) is given by

l Subscript m Baseline left-parenthesis bold-italic beta comma theta right-parenthesis almost-equals minus one-half log left-parenthesis theta Superscript s Baseline right-parenthesis minus one-half log left-parenthesis StartAbsoluteValue bold upper H 22 left-parenthesis bold-italic beta comma bold-italic gamma overTilde comma theta right-parenthesis EndAbsoluteValue right-parenthesis minus l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma overTilde comma theta right-parenthesis

The maximization of this approximate likelihood is a doubly iterative process that alternates between the following two steps:

  • For a provisional value of theta, PROC PHREG computes the best linear unbiased predictors (BLUP) of bold-italic beta and bold-italic gamma by maximizing the penalized partial log likelihood l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis. This constitutes the inner loop.

  • For bold-italic beta and bold-italic gamma fixed at the BLUP values, PROC PHREG estimates theta by maximizing the approximate marginal likelihood l Subscript m Baseline left-parenthesis bold-italic beta comma theta right-parenthesis. This constitutes the outer loop.

The outer loop is iterated until the difference between two successive estimates of theta is small.

The ML estimate of theta is

ModifyingAbove theta With caret equals StartFraction ModifyingAbove bold-italic gamma With caret prime ModifyingAbove bold-italic gamma With caret plus normal t normal r normal a normal c normal e left-parenthesis bold upper H 22 Superscript negative 1 Baseline right-parenthesis Over s EndFraction

The variance for ModifyingAbove theta With caret is

normal v normal a normal r left-parenthesis ModifyingAbove theta With caret right-parenthesis equals 2 ModifyingAbove theta With caret left-bracket s plus StartFraction 1 Over ModifyingAbove theta With caret squared EndFraction normal t normal r normal a normal c normal e left-parenthesis bold upper H 22 Superscript negative 1 Baseline bold upper H 22 Superscript negative 1 Baseline right-parenthesis minus StartFraction 2 Over ModifyingAbove theta With caret EndFraction normal t normal r normal a normal c normal e left-parenthesis bold upper H 22 Superscript negative 1 Baseline right-parenthesis right-bracket Superscript negative 1

The REML estimation of theta is obtained by replacing left-parenthesis bold upper H 22 right-parenthesis Superscript negative 1 by left-parenthesis bold upper H Superscript negative 1 Baseline right-parenthesis Subscript 22.

The inverse of the final bold upper H matrix is used as the variance estimate of left-parenthesis ModifyingAbove bold-italic beta With caret comma ModifyingAbove bold-italic gamma With caret right-parenthesis prime.

The final BLUP estimates of the random components gamma 1 comma ellipsis comma gamma Subscript s Baseline can be displayed by using the SOLUTION option in the RANDOM statement. Also displayed are estimates of the lognormal frailties, which are the exponentiated estimates of the BLUP estimates.

Wald-Type Tests for Penalized Models

Let bold upper I be the negative Hessian of the partial log likelihood l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis:

bold upper I equals Start 2 By 2 Matrix 1st Row 1st Column bold upper I 11 2nd Column bold upper I 12 2nd Row 1st Column bold upper I 21 2nd Column bold upper I 22 EndMatrix

where bold upper I 11 equals minus StartFraction partial-differential squared l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis Over partial-differential bold-italic beta squared EndFraction comma bold upper I 12 equals bold upper I prime 21 equals minus StartFraction partial-differential squared l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis Over partial-differential bold-italic beta partial-differential bold-italic gamma EndFraction, and bold upper I 22 equals minus StartFraction partial-differential squared l Subscript normal p normal a normal r normal t normal i normal a normal l Baseline left-parenthesis bold-italic beta comma bold-italic gamma right-parenthesis Over partial-differential bold-italic gamma squared EndFraction. Write bold-italic tau prime equals left-parenthesis bold-italic beta prime comma bold-italic gamma Superscript prime Baseline right-parenthesis prime. The Wald-type chi-square statistic for testing upper H 0 colon bold upper C bold-italic tau equals bold 0 is

left-parenthesis bold upper C ModifyingAbove bold-italic tau With caret right-parenthesis prime left-parenthesis bold upper C bold upper H Superscript negative 1 Baseline bold upper C prime right-parenthesis Superscript negative 1 Baseline left-parenthesis bold upper C ModifyingAbove bold-italic tau With caret right-parenthesis

Let bold upper H be the negative Hessian of the penalized partial log likelihood l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma theta right-parenthesis at the ML estimate theta; that is, bold upper H equals StartFraction partial-differential squared Over partial-differential bold-italic beta partial-differential bold-italic gamma EndFraction l Subscript p Baseline left-parenthesis bold-italic beta comma bold-italic gamma comma ModifyingAbove theta With caret right-parenthesis. Let bold upper V equals bold upper H Superscript negative 1 Baseline bold upper I bold upper H Superscript negative 1. Gray (1992) recommends the following generalized degrees of freedom for the Wald test:

normal upper D normal upper F equals normal t normal r normal a normal c normal e left-bracket left-parenthesis bold upper C bold upper H Superscript negative 1 Baseline bold upper C Superscript prime Baseline right-parenthesis Superscript negative 1 Baseline bold upper C bold upper V bold upper C Superscript prime Baseline right-parenthesis right-bracket

See Therneau and Grambsch (2000, Section 5.8) for a discussion of this Wald-type test.

PROC PHREG uses the label "Adjusted DF" to represent this generalized degrees of freedom in the output.

Last updated: March 08, 2022