Firth’s penalized likelihood for proportional hazards regressions for complex surveys
Section 4. Summary

Table of contents

Firth’s penalized likelihood is useful for obtaining maximum likelihood estimates from a monotone likelihood from proportional hazards regression models. We proposed a weight scaling method and demonstrated that Firth’s penalized likelihood using the scaled weights have some desirable properties for complex surveys. A simulation study shows that estimated biases in point estimates and standard errors using the scaled weights are lower than estimated biases using the unscaled weights. Although Firth’s penalized likelihood produces “good” estimates in most simulated data sets, there are a few data sets in which the Firth penalty failed to produce “good” convergences. The Firth penalized likelihood that uses scaled weights successfully corrected for a monotone likelihood when we estimated hazard rates for heart attacks using a data set from the NHEFS. Although the numeric results are quite encouraging, further research is needed to derive asymptotic distributions of the estimators obtained by using Firth’s penalized likelihood.

We recommend the unpenalized likelihood when convergence is not an issue, but we recommend Firth’s penalized likelihood using the scaled weights when a monotone likelihood is encountered in fitting proportional hazards regression models for complex surveys.

Acknowledgements

I am grateful to Ying So, Randy Tobias, and Ed Huddleston at SAS Institute Inc. for their valuable assistance in the preparation of this article. I would also like to thank the two anonymous referees and the associate editor for their constructive suggestions.

Appendix 1

Consistency of the Firth penalized likelihood estimator

The estimators in Section 2 are defined as the solution to a system of equations that are constructed by using the score functions from proportional hazards regression models. In this appendix, we show that under certain regularity conditions these estimators are design consistent. Properties of estimators that are solutions to a set of estimating equations are well studied in the survey literature. For example, see Binder (1983), Godambe and Thompson (1986), and Fuller (2009, Section 1.3.4).

However, the estimating equations for proportional hazards regression models are more complex than the estimating equations for generalized linear models because the score functions involve weighted sums over the sampled units. Binder (1992) and Lin (2000) showed that the estimators obtained by solving the estimating equations for proportional hazards regression models are consistent. In this appendix, we follow arguments similar to those of Lin (2000) and Andersen and Gill (1982).

Several technical assumptions are necessary to show that the point estimates are consistent. We need assumptions about the estimating equations, the finite population, and the sample design – to whit:

The functions defining the estimating equations should be smooth and convex.
The finite population should be such that the moments for population quantities that are used in defining the estimating equations exists.
The sample design should be such that the Narain-Horvitz-Thompson (NHT) estimators (Rao, 2005) for the population totals are well behaved.

All these assumptions are common in the sample survey literature; for example, see Fuller (2009). The score functions for proportional hazards regression models involve ratios of means of exponential functions that are infinitely differentiable.

Let $U_{N}$ and $F_{N}$ denote, respectively, the index set and values for the $N^{th}$ finite population in a sequence of populations indexed by $N,$ and let $A_{N}$ be a sample of size $n$ from $U_{N} .$ To study large sample properties for sample-based estimators, we assume sequences of population and samples such that $N \to \infty$ and $(N - n) \to \infty,$ keeping the sampling fraction, $\frac{n}{N},$ fixed.

Assume $F_{N} = {(t_{i}, Δ_{i}, Z_{i} (\cdot))}_{i =1}^{N}$ is an independent random sample of size $N$ from the joint distribution of $(T, Δ, Z (\cdot)),$ where $t$ is the failure time or the censoring time, whichever is less; $Δ = 1$ if the failure time is less than the censoring time and 0 otherwise; and $Z (\cdot)$ is a vector of possibly time-varying explanatory variables.

Let $β$ be a set of regression parameters for the superpopulation that is defined by the joint distribution of $(T, Δ, Z (\cdot)) .$ Let $β_{N}$ be a set of finite population parameters obtained by solving the estimating equations when all $N$ units in the population are observed, and let ${\hat{β}}_{N}$ be an estimator of $β_{N}$ that is obtained by solving the weighted estimating equations by using only the sampled units. Our objective is to show that ${\hat{β}}_{N}$ approaches $β_{N}$ and that they both approach $β$ as the sample size and population size increase.

Consider the estimating equations that correspond to Firth’s penalized likelihood described in Section 2. For simplicity, we write these equations when there are no tied events. To further simplify notation, we write each component of the estimating equations separately. The finite population parameters, $β_{N},$ are a solution to the penalized partial likelihood score function, $U_{N} (β) = {(U_{N, 1} (β), \dots, U_{N, P} (β))}^{'},$ where

$\begin{array}{l} U_{N, p} (β) & = N^{- 1} \sum_{i \in U_{N}} Δ_{i} [Z_{i} (t_{i}) - \frac{S_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} + 0.5 tr (\begin{array}{l} {[N^{- 1} \sum_{i \in U_{N}} Δ_{i} {\frac{S^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})} - (\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'}}]}^{- 1} \\ {\begin{cases} (\frac{Q_{p}^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})}) \\ - (\frac{Q_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'} \\ - (\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{Q_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'} \end{cases}} \end{array})] \end{array}$

where

$S^{(a)} (β, t) = N^{- 1} \sum_{i \in U_{N}} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) {[Z_{i} (t)]}^{\otimes a}$

$Q_{p}^{(a)} (β, t) = N^{- 1} \sum_{i \in U_{N}} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) Z_{i, p} (t) {[Z_{i} (t)]}^{\otimes a}$

and where $a = 0, 1, 2;$ $Z_{i} (t) = {(Z_{i, 1} (t), \dots, Z_{i, p} (t))}^{'};$ $S^{(1)} (β, t) = (S_{1}^{(1)} (β, t), \dots, S_{P}^{(1)} (β, t))^{'};$ $tr (\cdot)$ denotes the trace of a matrix; $I (\cdot)$ denotes the indicator function; $p = 1, 2, \dots, P;$ and $P$ is the number of regression parameters. Note that $S^{(a)} (β, t)$ and $Q_{p}^{(a)} (β, t)$ depend on $N,$ although the notation does not reflect this for reasons of simplicity.

In defining the score function for the penalized likelihood, we assume that the information matrix for the finite population, $I_{N} (β, t),$ is always positive definite.

However, in any realistic situation, not all units in the finite population are available. Let a sample $A_{N}$ be selected by using a probability design that assigns a nonzero selection probability, $π_{i},$ to every unit in the population. Let $w_{i} = π_{i}^{- 1}$ be the design weight. A sample-based estimator, ${\hat{β}}_{N},$ is obtained by solving the estimated penalized partial likelihood score equations. Assuming that $N$ is known, a sample-based estimator for $U_{N, p} (β)$ is

${\hat{U}}_{N, p} (β) = N^{- 1} \sum_{i \in A_{N}} w_{i} Δ_{i} [Z_{i} (t_{i}) - \frac{{\hat{S}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} + 0.5 tr (\begin{array}{l} {[N^{- 1} \sum_{i \in A_{N}} w_{i} Δ_{i} {\frac{{\hat{S}}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - (\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'}}]}^{- 1} \\ {\begin{cases} (\frac{{\hat{Q}}_{p}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) \\ - (\frac{{\hat{Q}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'} \\ - (\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{Q}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'} \end{cases}} \end{array})]$

where

${\hat{S}}^{(a)} (β, t) = N^{- 1} \sum_{i \in A_{N}} w_{i} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) {[Z_{i} (t)]}^{\otimes a}$

${\hat{Q}}_{p}^{(a)} (β, t) = N^{- 1} \sum_{i \in A_{N}} w_{i} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) Z_{i, p} (t) {[Z_{i} (t)]}^{\otimes a}$

are the NHT estimators for $S^{(a)} (β, t)$ and $Q_{p}^{(a)} (β, t),$ respectively.

Because ${\hat{S}}^{(a)} (β, t)$ and ${\hat{Q}}_{p}^{(a)} (β, t)$ use weighted sums over sampled units, we need techniques defined in Lin (2000) to study large sample properties of these estimators. Define $G_{i} (t) = Δ_{i} I (t_{i} \leq t),$ $G (t) = N^{- 1} \sum_{i \in U_{N}} G_{i} (t),$ and $\hat{G} (t) = \sum_{i \in A_{N}} w_{i} G_{i} (t) .$ Then the finite population score functions can be written using stochastic integration,

$U_{N, p} (β) = N^{- 1} \sum_{i \in U_{N}} \int_{0}^{\infty} [Z_{i} (t_{i}) - \frac{S_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} + 0.5 tr (\begin{array}{l} {[N^{- 1} \sum_{i \in A_{N}} \int_{0}^{\infty} {\frac{S^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})} - (\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'}} d G (t_{i})]}^{- 1} \\ {\begin{cases} (\frac{Q_{p}^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(2)} (β, t_{i})}{S^{(0)} (β, t_{i})}) \\ - (\frac{Q_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'} \\ - (\frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})}) {(\frac{Q_{p}^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})} - \frac{Q_{p}^{(0)} (β, t_{i})}{S^{(0)} (β, t_{i})} \frac{S^{(1)} (β, t_{i})}{S^{(0)} (β, t_{i})})}^{'} \end{cases}} \end{array})] d G_{i} (t)$

and the sample-based score functions are

${\hat{U}}_{N, p} (β) = N^{- 1} \sum_{i \in U_{N}} \int_{0}^{\infty} I (i \in A_{n}) w_{i} [Z_{i} (t_{i}) - \frac{{\hat{S}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} + 0.5 tr (\begin{array}{l} {[N^{- 1} \sum_{i \in U_{N}} \int_{0}^{\infty} I (i \in A_{N}) w_{i} {\frac{{\hat{S}}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - (\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'}} d G_{i} (t)]}^{- 1} \\ {\begin{cases} (\frac{{\hat{Q}}_{p}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(2)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) \\ - (\frac{{\hat{Q}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'} \\ - (\frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})}) {(\frac{{\hat{Q}}_{p}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} - \frac{{\hat{Q}}_{p}^{(0)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})} \frac{{\hat{S}}^{(1)} (β, t_{i})}{{\hat{S}}^{(0)} (β, t_{i})})}^{'} \end{cases}} \end{array})] d G_{i} (t) .$

Note that the quantities $S^{(a)}$ and $Q_{p}^{(a)}$ are simply means over finite population quantities. Define the limits of these means as follows:

$\begin{array}{l} s^{(a)} (β, t) & := \lim_{N \to \infty} S^{(a)} (β, t) \\ = \lim_{N \to \infty} N^{- 1} \sum_{i \in U_{N}} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) {[Z_{i} (t)]}^{\otimes a} \\ q_{p}^{(a)} (β, t) & := Q_{p}^{(a)} (β, t) \\ = N^{- 1} \sum_{i \in U_{N}} I (t_{i} \geq t) \exp (β^{'} Z_{i} (t)) Z_{i, p} (t) {[Z_{i} (t)]}^{\otimes a} \\ g (t) & := \lim_{N \to \infty} N^{- 1} \sum_{i \in U_{N}} G_{i} (t) \\ α & := \lim_{N \to \infty} N^{- 1} \sum_{i \in U_{N}} \int_{0}^{\infty} Z_{i} (t) d G_{i} (t) . \end{array}$

Thus the finite population score function, $U_{N, p} (β),$ converges to the superpopulation score function $u_{N, p} (β),$ where

$u_{N, p} (β) = α - \int_{0}^{\infty} \frac{s_{p}^{(1)} (β, t)}{s^{(0)} (β, t)} d g (t) + 0.5 tr (\begin{array}{l} \int_{0}^{\infty} {[\int_{0}^{\infty} {\frac{s^{(2)} (β, t)}{s^{(0)} (β, t)} - (\frac{s^{(1)} (β, t)}{s^{(0)} (β, t)}) {(\frac{s^{(1)} (β, t)}{s^{(0)} (β, t)})}^{'}} d g (t)]}^{- 1} \\ {\begin{cases} {(\frac{q_{p}^{(2)} (β, t)}{s^{(0)} (β, t)} - \frac{q_{p}^{(0)} (β, t)}{s^{(0)} (β, t)} \frac{s^{(2)} (β, t)}{s^{(0)} (β, t)})}_{_{_{}}}^{^{^{}}} \\ - (\frac{q_{p}^{(1)} (β, t)}{s^{(0)} (β, t)} - \frac{q_{q}^{p} (β, t)}{s^{(0)} (β, t)} \frac{s^{(1)} (β, t)}{s^{(0)} (β, t)}) {(\frac{s^{(1)} (β, t)}{s^{(0)} (β, t)})}^{'} \\ - (\frac{s^{(1)} (β, t)}{s^{(0)} (β, t)}) {(\frac{q_{p}^{(1)} (β, t)}{s^{(0)} (β, t)} - \frac{q_{p}^{(0)} (β, t)}{s^{(0)} (β, t)} \frac{s^{(1)} (β, t)}{s^{(0)} (β, t)})}^{'} \end{cases}} d g (t) \end{array}) .$

Now assume that the population quantities, $Z_{i},$ that are used to define the score functions have finite moments and the sequence of sample designs is such that any smooth functions of NHT estimators are consistent. Because $U_{N} (β)$ is a smooth function of population totals, and each total is estimated by using a NHT estimator, ${\hat{U}}_{N} (β)$ is design-consistent for $U_{N} (β) .$ That is, $(U_{N} (β) - {\hat{U}}_{N} (β)) | F_{N} = o (1) .$ Therefore, by using arguments similar to Lin (2000) and Andersen and Gill (1982), it can be shown that $β_{N}$ and ${\hat{β}}_{N}$ converge to the same limit.

Because $n / N$ is fixed, $\sum_{A_{N}} w_{i}$ is the NHT estimator for $N,$ and $\hat{U} (β)$ is a consistent estimator (not necessarily unbiased) of 0, both $\hat{U} (β)$ and $(n / N) \hat{U} (β)$ converge to the same limit with the same order of convergence. It is straightforward to show that $(n / N) \hat{U} (β)$ and $(n / \sum_{A_{N}} w_{i}) \hat{U} (β),$ the estimating equations that use the scaled weights, have the same expectation.

Appendix 2

SAS program to obtain the Firth penalized likelihood estimates

The SAS statements at the end of this section fit a proportional hazards regression model using the scaled weights in Firth’s penalized likelihood. The PROC statement invokes the procedure, and the VARMETHOD = JK option requests the jackknife variance estimation method. You can also specify VARMETHOD = TAYLOR, VARMETHOD = BRR, or VARMETHOD = BOOT to request the Taylor series linearized, balanced repeated replication, or bootstrap replication variance estimation method, respectively. The DETAILS sub-option of the VARMETHOD = JK option prints estimates from each replicate sample along with the convergence status. The WEIGHT statement specifies the sampling weights, the STRATA statement specifies the strata, and the CLUSTER statement specifies the PSUs. The MODEL statement specifies the analysis model. The FIRTH option in the MODEL statement requests Firth’s penalized likelihood. The two HAZARDRATIO statements requests hazard ratios for blood cholesterol and smoking, respectively. The ODS OUTPUT statement stores replicate estimates and convergence status from each replicate in the SAS data set RepEstimatesFirth. This data set is useful for checking the convergence status of every replicate sample.

Start of text box

proc surveyphreg data = NHEFS varmethod=jk (details);

class

Gender HighBloodChol Race Smoker;

weight

ObservationWeight;

strata

Stratum;

cluster

PSU;

model

EventTime*HeartAttack(2) = Income HighBloodChol Smoker Race Gender Race*Gender / firth;

hazardratio

HighBloodChol;

hazardratio

Smoker;

ods output

repestimates=RepEtimatesFirth;

run;

End of text box

References

Andersen, P.K., and Gill, R.D. (1982). Cox’s regression model counting process: A large sample study. Annals of Statistics, 10, 1100-1120.

Bender, R., Augustin, T. and Blettner, M. (2005). Generating survival times to simulate Cox proportional hazards models. Statistics in Medicine, 24, 1713-1723.

Binder, D.A. (1983). On the variances of asymptotically normal estimators from complex surveys. International Statistical Review, 51, 279-292.

Binder, D.A. (1990). Fitting Cox’s proportional hazards models from survey data. Proceedings of the Survey Research Methods Section, American Statistical Association, 342-347.

Binder, D.A. (1992). Fitting Cox’s proportional hazards models from survey data. Biometrika, 79, 139-147.

Binder, D.A., and Patak, Z. (1994). Use of estimating functions for estimation from complex surveys. Journal of the American Statistical Association, 89(427),1035-1043.

Boudreau, C., and Lawless, J.F. (2006). Survival analysis based on the proportional hazards model and survey data. Canadian Journal of Statistics, 34, 203-216.

Breslow, N.E. (1974). Covariance analysis of censored survival data. Biometrics, 30, 89-99.

Chambers, R.L., and Skinner, C.J. (2003). Analysis of Survey Data. Chichester, UK: John Wiley & Sons, Inc.

Chambless, L.E., and Boyle, K.E. (1985). Maximum likelihood methods for complex sample data: Logistic regression and discrete proportional hazards models. Communications in Statistics - Theory and Methods, 14, 1377-1392.

Cox, D.R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society, Series B, (with discussion), 34, 187-220.

Cox, D.R. (1975). Partial likelihood. Biometrika, 62, 269-276.

Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80, 27-38.

Fuller, W.A. (1975). Regression analysis for sample survey. Sankhyā, Series C, 37, 117-132.

Fuller, W.A. (2009). Sampling Statistics. Hoboken, NJ: John Wiley & Sons, Inc.

Godambe, V.P., and Thompson, M.E. (1986). Parameters of superpopulation and survey population: Their relationships and estimation. International Statistical Review, 54, 127-138.

Heinze, G. (1999). The Application of Firth’s Procedure to Cox and Logistic Regression. Technical Report 10, updated January 2001, Department of Medical Computer Sciences, Section of Clinical Biometrics, University of Vienna.

Heinze, G., and Schemper, M. (2001). A solution to the problem of monotone likelihood in Cox regression. Biometrics, 51, 114-119.

Heinzel, G., Rüdiger, A. and Schilling, R. (2002). Spectrum and Spectral Density Estimation by the Discrete Fourier Transform (DFT), Including a Comprehensive List of Window Functions and Some New Flat-Top Windows. Tech. rep., Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Teilinstitut Hannover.

Kish, L., and Frankel, M.R. (1974). Inference from complex samples. Journal of the Royal Statistical Society, Series B, 36, 1-37.

Korn, E.L., and Graubard, B.I. (1999). Analysis of Health Surveys. New York: John Wiley & Sons, Inc.

Lee, E.W., Wei, L.J. and Amato, D.A. (1992). Cox-type regression analysis for large numbers of small groups of correlated failure time observations. Survival Analysis: State of the Art, (Eds., J.P. Klein and P.K. Goel), Dordrecht, Netherlands: Kluwer Academic, 237-247.

Lin, D.Y. (2000). On fitting Cox’s proportional hazards models to survey data. Biometrika, 87, 37-47.

Lin, D.Y., and Wei, L.J. (1989). The robust inference for the proportional hazards model. Journal of the American Statistical Association, 84, 1074-1078.

Mukhopadhyay, P.K. (2010). Not hazardous to your health: Proportional hazards modeling for survey data with the SURVEYPHREG procedure. Proceedings of the SAS Global Forum 2010 Conference. Cary, NC: SAS Institute Inc. http://support.sas.com/resources/papers/proceedings10/254-2010.pdf.

Mukhopadhyay, P.K., An, A.B., Tobias, R.D. and Watts, D.L. (2008). Try, try again: Replication-based variance estimation methods for survey data analysis in SAS 9.2. Proceedings of the SAS Global Forum 2008 Conference. Cary, NC: SAS Institute Inc. http://www2.sas.com/proceedings/forum2008/367-2008.pdf.

Pfeffermann, D. (1993). The role of sampling weights when modeling survey data. International Statistical Review, 61, 317-337.

Rao, J.N.K. (2005). Interplay between sample survey theory and practice: An appraisal. Survey Methodology, 31, 2, 117-138. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2005002/article/9040-eng.pdf.

SAS Institute Inc. (2018). SAS/STAT 15.1 User’s Guide. Cary, NC: SAS Institute Inc. http://go.documentation.sas.com/?docsetId=statug&docsetTarget=titlepage.htm&docsetVersion=15.1&locale=en.

Shao, J., and Tu, D. (1995). The Jackknife and Bootstrap. New York: Springer-Verlag.

Sitter, R.R. (1992). Comparing three bootstrap methods for survey data. Canadian Journal of Statistics, 20, 135-154.

Wolter, K.M. (2007). Introduction to Variance Estimation. 2^nd Ed. New York: Springer.

ISSN : 1492-0921

Editorial policy

Survey Methodology publishes articles dealing with various aspects of statistical development relevant to a statistical agency, such as design issues in the context of practical constraints, use of different data sources and collection techniques, total survey error, survey evaluation, research in survey methodology, time series analysis, seasonal adjustment, demographic studies, data integration, estimation and data analysis methods, and general survey systems development. The emphasis is placed on the development and evaluation of specific methodologies as applied to data collection or the data themselves. All papers will be refereed. However, the authors retain full responsibility for the contents of their papers and opinions expressed are not necessarily those of the Editorial Board or of Statistics Canada.

Submission of Manuscripts

Survey Methodology is published twice a year in electronic format. Authors are invited to submit their articles in English or French in electronic form, preferably in Word to the Editor, (statcan.smj-rte.statcan@canada.ca, Statistics Canada, 150 Tunney’s Pasture Driveway, Ottawa, Ontario, Canada, K1A 0T6). For formatting instructions, please see the guidelines provided in the journal and on the web site (www.statcan.gc.ca/SurveyMethodology).

Note of appreciation

Canada owes the success of its statistical system to a long-standing partnership between Statistics Canada, the citizens of Canada, its businesses, governments and other institutions. Accurate and timely statistical information could not be produced without their continued co-operation and goodwill.

Standards of service to the public

Statistics Canada is committed to serving its clients in a prompt, reliable and courteous manner. To this end, the Agency has developed standards of service which its employees observe in serving its clients.

Copyright

Published by authority of the Minister responsible for Statistics Canada.

Use of this publication is governed by the Statistics Canada Open Licence Agreement.

Catalogue No. 12-001-X

Frequency: Semi-annual

Ottawa

Date modified:: 2020-12-15

Language selection

Search and menus

Search

Firth’s penalized likelihood for proportional hazards regressions for complex surveys
Section 4. Summary

Acknowledgements

Appendix 1

Consistency of the Firth penalized likelihood estimator

Appendix 2

SAS program to obtain the Firth penalized likelihood estimates

References

Firth’s penalized likelihood for proportional hazards regressions for complex surveys Section 4. Summary

Acknowledgements

Appendix 1

Consistency of the Firth penalized likelihood estimator

Appendix 2

SAS program to obtain the Firth penalized likelihood estimates

References

Editorial policy

Submission of Manuscripts

Note of appreciation

Standards of service to the public

Copyright

Firth’s penalized likelihood for proportional hazards regressions for complex surveys
Section 4. Summary