Firth’s penalized likelihood for proportional hazards regressions for complex surveys
Section 1. Introduction

• Previous
• Next

The Cox proportional hazards regression model (Cox, 1972) is widely used to analyze survival data. It is a semiparametric model that explains the effect of explanatory variables on hazard rates. The model assumes a linear form for the effect of the explanatory variables but allows an unspecified form for the underlying survivor function. The parameters of the model are estimated by maximizing a partial likelihood (Cox, 1972, 1975).

For estimating canonical parameters in the exponential family distributions, Firth (1993) suggested multiplying the likelihood by the Jeffreys prior to obtain a maximum likelihood estimate that is first-order unbiased. The penalized likelihood is of the form

$$L_p\left(\beta \right)\mathrm{<mo>=</mo>}L\left(\beta \right){\left|I\left(\beta \right)\right|}^{0.5}$$

where $L\left(\beta \right)$ is the unpenalized likelihood, $I$ is the information matrix, and $\beta$ is a vector of regression parameters. Firth’s penalized likelihood is a very useful technique in practice, not only to reduce bias but also to correct for monotone likelihoods.

Proportional hazards regression models often suffer from monotone likelihoods, in which the likelihood converges to a finite value but at least one parameter diverges (Heinze, 1999). Firth’s penalized likelihood is also used to correct monotone likelihoods and to obtain parameter estimates that converge (Heinze, 1999; Heinze and Schemper, 2001; Heinzel, Rüdiger and Schilling, 2002).

Although Firth’s penalized likelihood is useful for reducing biases and for obtaining estimates from monotone likelihoods, the penalized likelihood is not studied for complex surveys involving unequal weights. It is reasonable to use a weighted likelihood for complex surveys to compensate for unequal weighting (Fuller, 1975; Binder and Patak, 1994). Survey data sets commonly include design weights or analysis weights for which the sum of the weights is an estimator of the population size. However, these unscaled weights will not appropriately scale the information matrix that is used in the penalty term. It is desirable for proportional hazards regression parameters for survey data to have the following two properties:

• Invariance: The point estimates and standard errors for the regression parameters should be invariant to the scale of the weights.
• Closeness: The Taylor linearized variance for the estimated regression parameters should be close to the delete-one jackknife variance.

In this article, we first show that if the Firth correction is not used, then both the invariance and closeness are satisfied; but if the Firth correction is used with the unscaled weights, then the point estimates and the standard errors are not invariant to the scale of the weights. That is, if the weights are multiplied by a constant and the Firth correction is used, then the point estimates and standard errors will be different. We then propose a commonsense weight scaling method and demonstrate that the Firth correction using the scaled weights has both properties. The only difference between the scaled and unscaled weights is that the sum of the scaled weights is equal to the sample size, but the sum of the unscaled weights is an estimator of the population size.

1.1  Example that uses unscaled weights

 We used a data set from a study of 65 myeloma patients who were treated with alkylating agents (Lee, Wei and Amato, 1992) to demonstrate the properties of Firth’s penalized likelihood that uses unscaled weights. Survival times in months were recorded for each patient. Patients who were alive after the study period were considered to be censored. The following variables were available for each patient:

• Time: Survival time in months,
• Vstatus: Patient status, zero or one, indicating whether the patient was alive or dead, respectively,
• LogBUN: Log of blood urea nitrogen level,
• HGB: Blood hemoglobin level.

To create a monotone likelihood, we added a new explanatory variable, Contrived, such that its value at all event times is the largest of all values in the risk set (see the example “Firth’s Correction for Monotone Likelihood” in “The PHREG Procedure” in SAS Institute Inc. (2018)). The variable Contrived has the value 1 if the observed survival time is less than or equal to 65; otherwise it has the value 0.

To demonstrate the effect of weights in Firth’s penalized likelihood, we created three weight variables, w1, w3, and w5, with the values of 1, 1,000, and 100,000 for each observation, respectively. Proportional hazards regression parameters are estimated by maximizing a weighted likelihood as described in Section 1.2. Because $w1$ has the value 1 for all observations, using $w1$ in the analysis is equivalent to performing the unweighted analysis.

We fitted the following two proportional hazards models using the PHREG procedure in SAS/STAT^® (see “The PHREG Procedure” in SAS Institute Inc. (2018)):

$$\lambda \left(t\mathrm{,}Z\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left({\beta }_1\text{LogBUN}+{\beta }_2\text{HGB}\right)$$

$$\lambda \left(t\mathrm{,}Z\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left({\beta }_1\text{LogBUN}+{\beta }_2\text{HGB}+{\beta }_3\text{Contrived}\right)$$

where $\lambda \left(t\right)$ and ${\lambda }_0\left(t\right)$ are the hazard function and the baseline hazard function, respectively. Firth’s penalized likelihood is not required in order to fit the first model without Contrived (the likelihood converged in three iteration steps), but the second model containing the variable Contrived does not converge without the Firth penalty in the likelihood. Table 1.1 displays the value of the likelihood and the three regression coefficients for 14 iterations. Although the objective function and the coefficients for LogBun and HGB converge to a finite value after the fourth iteration, the coefficients for Contrived diverges. This is an example of a monotone likelihood for the variable Contrived. Because of this monotonicity, Firth’s penalized likelihood must be used to fit the second model containing Contrived.

If Contrived is not used as an explanatory variable, then all three sets of weights produce the same point estimates and Taylor linearized variance estimates (Table 1.2). The delete-one jackknife variance estimates are also the same for all three sets of weights. Thus, the point estimates and the standard errors are invariant to the scale of the weights when the Firth correction is not used.

However, if the unscaled weights are used, then the point estimates for Contrived are not invariant to the scale of the weights. Table 1.3 displays the parameter estimates for three sets of weights when Contrived is used as an explanatory variable (and Firth’s penalized likelihood is applied). Because the likelihood is not monotone (Table 1.1) for LogBun and HGB, the point estimates for these two coefficients are not affected by the scale of the weights.

If Contrived is not used as an explanatory variable, then the ratio of jackknife standard errors to Taylor linearized standard errors is 1.13 and 1.10 for all three sets of weights for the variables LogBun and HGB, respectively. Thus the ratio of the jackknife variance to the Taylor linearized variance for the unpenalized likelihood is invariant to the scale of weights, and it is reasonable to expect the ratio to be invariant when the penalized likelihood is used.

1.2  A brief review of point and variance estimates for regression parameters for finite populations

Before we discuss the weight scaling method, we briefly review point and variance estimates for regression parameters for proportional hazards regression of complex surveys involving unequal weights. Lin and Wei (1989); Binder (1990, 1992); Lin (2000); and Boudreau and Lawless (2006) discussed pseudo-maximum likelihood estimation of proportional hazard regression parameters for survey data. For a more general description for estimating regression parameters for complex surveys, see Kish and Frankel (1974); Godambe and Thompson (1986); Pfeffermann (1993), Korn and Graubard (1999, Chapter 3), Chambers and Skinner (2003, Chapter 2), and Fuller (2009, Section 6.5). Wolter (2007) described different variance estimation techniques for survey data.

Let $U_N\mathrm{<mo>=</mo>}\left\{\mathrm{1,}\mathrm{2,}\dots \mathrm{,}N\right\}$ be the set of indices and let $F_N$ be the set of values for a finite population of size $N.$ The survival time of each member of the finite population is assumed to follow its own hazard function, ${\lambda }_i\left(t\right),$ expressed as

$${\lambda }_i\left(t\right)\mathrm{<mo>=</mo>}\lambda \left(t\mathrm{<mo>;</mo>}{Z}_i\left(t\right)\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left(Z_i^{\prime }\left(t\right)\beta \right)$$

where ${\lambda }_0\left(t\right)$ is an arbitrary and unspecified baseline hazard function, $Z_i\left(t\right)$ is a vector of size $P$ of explanatory variables for the $i^{\text{th}}$ unit at time $t,$ and $\beta$ is a vector of unknown regression parameters.

The partial likelihood function introduced by Cox (1972, 1975) eliminates the unknown baseline hazard ${\lambda }_0\left(t\right)$ and accounts for censored survival times. If the entire population is observed, then this partial likelihood function can be used to estimate $\beta .$ Let ${\beta }_N$ be the desired estimator.

Assuming a working model with uncorrelated responses, ${\beta }_N$ is obtained by maximizing the partial log likelihood,

$$l_N\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in U_N}}\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}$$

with respect to $\beta ,$ where $\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)$ is Cox’s partial likelihood function.

Assume that a probability sample $A_N$ is selected from the finite population $U_N.$ Let ${\pi }_i$ be the selection probability and $w_i\left(={\pi }_i^{-1}\right)$ be the sampling weight for unit $i.$ Further assume that explanatory variables $Z_i\left(t\right)$ and survival time $t_i$ are available for every unit in sample $A_N.$ A design unbiased estimator for the finite population log likelihood is

$$l\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in A_N}}{\pi }_i^{-1}\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in A_N}}{w}_i\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}.$$

A sample-based estimator ${\widehat{\beta }}_N$ for the finite population quantity ${\beta }_N$ can be obtained by maximizing the partial pseudo-log likelihood $l\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)$ with respect to $\beta .$ The design-based variance for ${\widehat{\beta }}_N$ is obtained by assuming that the set of finite population values $F_N$ is fixed.

The weighted Breslow likelihood can be expressed as

$$L\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \prod _{k\mathrm{<mo>=</mo>\; 1}}^K}\frac{\exp \left({\beta }^{\prime }{\displaystyle {\sum }_{D_k}{w}_i}{Z}_i\left(t\right)\right)}{{\left\{{\displaystyle {\sum }_{R_k}{w}_i}\exp \left({\beta }^{\prime }{Z}_i\left(t\right)\right)\right\}}^{{\displaystyle {\sum }_{D_k}{w}_i}}}$$

where $R_k$ is the risk set just before the $k^{\text{th}}$ ordered event time $t_{\left(k\right)},$ $D_k$ is the set of individuals who fail at the $t_{\left(k\right)},$ and $K$ is the number of distinct event times.

The point estimates for $\beta$ are obtained by maximizing $l\left(\beta \right)\mathrm{<mo>=</mo>}\log \left[L\left(\beta \right)\right].$

Although the weights are sufficient for estimating regression coefficients for the finite population, stratification and clustering information must also be used to estimate sampling variability. In order to estimate sampling variability, you can use either the Taylor series linearization method or a replication method.

1.2.1  Analytic variance estimator using the Taylor series linearization method

The Taylor series linearization method uses a sum of squares of the weighted score residuals to estimate the sampling variability.

Define $\overline{Z}\left(\beta \mathrm{,}t\right)\mathrm{<mo>=</mo>}\frac{S^{\left(1\right)}\left(\beta \mathrm{,}t\right)}{S^{\left(0\right)}\left(\beta \mathrm{,}t\right)},$ where

$$S^{\left(0\right)}\left(\beta ,t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{A_N}}{w}_{i}I\left(t_i\ge t\right)\exp \left({\beta }^{\prime }{Z}_i\left(t\right)\right)$$

and

$$S^{\left(1\right)}\left(\beta ,t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{A_N}}{w}_{i}I\left(t_i\ge t\right)\exp \left({\beta }^{\prime }{Z}_i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\right)Z_i\left(t\right).$$

The score residual for the $i^{\text{th}}$ subject is

$$\begin{array}{ll}{u}_i\left(\beta \right)=\hfill & {\Delta }_i\left\{Z_i\left(t_i\right)-\overline{Z}\left(\beta ,t_i\right)\right\}\hfill \\ \hfill & -{\displaystyle \sum _{j\in A_N}}\left[{\Delta }_j\frac{w_{j}I\left(t_i\ge t_j\right)\exp \left({\beta }^{\prime }{Z}_i\left(t_j\right)\right)}{S^{\left(0\right)}\left(\beta ,t_j\right)}\left\{Z_i\left(t_j\right)-\overline{Z}\left(\beta \mathrm{,}{t}_j\right)\right\}\right]\hfill \end{array}$$

where ${\Delta }_i$ is the event indicator.

Then the Taylor linearized variance estimator is

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{I}^{-1}\left(\widehat{\beta }\right)\text{G}{I}^{-1}\left(\widehat{\beta }\right)$$

where $I\left(\widehat{\beta }\right)$ is the observed information matrix and the $p\times p$ matrix $G$ is defined as

$$G\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\in A_N\mathrm{<mo>:</mo>}i\mathrm{<mo><</mo>}j}}\frac{{\pi }_i{\pi }_j-{\pi }_{ij}}{{\pi }_{ij}}{\left(\frac{{\widehat{u}}_i}{{\pi }_i}-\frac{{\widehat{u}}_j}{{\pi }_j}\right)}^{\prime }\left(\frac{{\widehat{u}}_i}{{\pi }_i}-\frac{{\widehat{u}}_j}{{\pi }_j}\right)$$

where ${\pi }_{ij}$ are the joint inclusion probabilities for units $i$ and $j.$

In particular, for stratified cluster designs in which the PSUs are selected by using a simple random sample without replacement, the $p\times p$ matrix $G$ reduces to

$$G\mathrm{<mo>=</mo>}{\displaystyle \sum _{h\mathrm{<mo>=</mo>\; 1}}^H}\frac{n_h\left(1-f_h\right)}{n_h-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^{n_h}}{\left(e_{hi+}-{\overline{e}}_{h\cdot \cdot }\right)}^{\prime }\left(e_{hi+}-{\overline{e}}_{h\cdot \cdot }\right)$$

where $e_{hi+}$ is the weighted sum of the score residuals, ${\widehat{u}}_{hij},$ in stratum $h$ and PSU $i;$ ${\overline{e}}_{h\mathrm{..}}$ is the mean of $e_{hi+};$ $n_h$ is the number of PSUs; and $f_h$ is the sampling fraction in stratum $h.$

These estimators are well studied in the sample survey literature. For example, Binder (1992) and Lin (2000) provide conditions under which $\widehat{\beta }$ and $\widehat{V}\left(\widehat{\beta }\right)$ are consistent. Chambless and Boyle (1985) derived the design-based variance and asymptotic normality for discrete proportional hazards models.

1.2.2  Replication variance estimator using the delete-one jackknife method

The jackknife method is a commonly used replication variance estimation method for complex surveys. To create replicates, it deletes (assigns a zero weight to) one PSU at a time from the full sample. In each replicate, the sampling weights of the remaining PSUs are modified by the jackknife coefficient ${\alpha }_r.$ The modified weights are called replicate weights.

Let PSU $i_r$ in stratum $h_r$ be omitted from the $r^{\text{th}}$ replicate; then the replicate weights and jackknife coefficients are given by

$$w_{hij}^{\left(r\right)}\mathrm{<mo>=</mo>}\{\begin{array}{ll}0\hfill & i\mathrm{<mo>=</mo>}{i}_r\text{and}h\mathrm{<mo>=</mo>}{h}_r\hfill \\ w_{hij}/{\alpha }_r\hfill & i\ne i_r\text{and}h\mathrm{<mo>=</mo>}{h}_r\hfill \\ w_{hij}\hfill & h\ne h_r\hfill \end{array}$$

and ${\alpha }_r\mathrm{<mo>=</mo>}\frac{n_{h_r}-1}{n_{h_r}},$ respectively, for all observation units $j$ in stratum $h$ and PSU $i.$ The number of PSUs in stratum $h_r$ is $n_{h_r}.$

The jackknife method can be applied to estimate variances for the estimated regression parameters for Cox’s model because the model parameters are solutions of a set of estimating equations that are smooth functions of totals (the corresponding score functions are given in Section 2). Properties of jackknife variance estimators for proportional hazard regression models are discussed in Shao and Tu (1995, Section 8.3).

To apply the jackknife method, model parameters are estimated by using the full sample and by using every replicate sample. Let $\widehat{\beta }$ be the estimated proportional hazards regression coefficients from the full sample, and let ${\widehat{\beta }}_r$ be the estimated regression coefficients from the $r^{\mathrm{th}}$ replicate. Then the covariance matrix of $\widehat{\beta }$ is estimated by

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R}{\alpha }_r\left({\widehat{\beta }}_r-\widehat{\beta }\right){\left({\widehat{\beta }}_r-\widehat{\beta }\right)}^{\prime }.$$

If the sampling fractions are not ignorable, then the covariance matrix of $\widehat{\beta }$ is estimated by

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R}{\alpha }_r\left(1-f_r\right)\left({\widehat{\beta }}_r-\widehat{\beta }\right){\left({\widehat{\beta }}_r-\widehat{\beta }\right)}^{\prime }$$

where $f_r\mathrm{<mo>=</mo>}\frac{n_{h_r}}{N_{h_r}}$ is the sampling fraction in stratum $h_r.$

In practice, both Taylor linearized variance and jackknife variance estimates are used to construct Wald $t$ confidence intervals with $R-H$ degrees of freedom, where $R$ is the number of PSUs (or the number of replicates) and $H$ is the number of strata.

It is straightforward to show that the jackknife variance estimator is algebraically equivalent to the Taylor linearized estimator for design linear estimators. But for design nonlinear estimators, such as the regression coefficients for proportional hazards regression models, the jackknife method tends to produce slightly higher variance estimates than the Taylor linearized method (Fuller, 2009).

Note that if the full sample estimate suffers from a monotone likelihood, then it is very likely that most replicate samples will also suffer from monotone likelihoods. This will results in many “unusable” replicate estimates.

Survey data analysis procedures in SAS/STAT support both Taylor linearized and replication variance estimation methods (Mukhopadhyay, An, Tobias and Watts, 2008).

• Previous
• Article main page
• Next