The PHREG Procedure

Proportional Rates/Means Models for Recurrent Events

Let $N(t)$ be the number of events experienced by a subject over the time interval $(0,t]$. Let $dN(t)$ be the increment of the counting process N over $[t,t+dt)$. The rate function is given by

\[  d\mu _{\bZ }(t)= E[dN(t) | Z(t)] = \mr {e}^{\bbeta \bZ (t)}d\mu _0(t)  \]

where $\mu _0(.)$ is an unknown continuous function. If the $\bZ $ are time independent, the rate model is reduced to the mean model

\[  \mu _ Z(t)= \mr {e}^{\bbeta Z} \mu _0(t)  \]

The partial likelihood for n independent triplets $(N_{i},Y_{i},\mb {Z}_{i}), i=1, \ldots , n$, of counting, at-risk, and covariate processes is the same as that of the multiplicative hazards model. However, a robust sandwich estimate is used for the covariance matrix of the parameter estimator instead of the model-based estimate.

Let $T_{ki}$ be the kth event time of the ith subject. Let $C_ i$ be the censoring time of the ith subject. The at-risk indicator and the failure indicator are, respectively,

\[  Y_ i(t) = I(C_ i\ge t) ~ ~ \mr {and}~ ~  \Delta _{ki}= I(T_{ki} \le C_ i)  \]


\[  S^{(0)}(\bbeta ,t) = \sum _{i=1}^ n Y_ i(t)\mr {e}^{\bbeta \bZ _ i(t)} ~ ~ \mr {and}~ ~  \bar{\bZ }(\bbeta ,t) = \frac{\sum _{i=1}^ n Y_ i(t)\mr {e}^{\bbeta \bZ _ i(t)} \bZ _ i(t)}{S^{(0)}(\bbeta ,t)}  \]

Let $\hat{\bbeta }$ be the maximum likelihood estimate of $\bbeta $, and let $\mc {I}(\hat{\bbeta })$ be the observed information matrix. The robust sandwich covariance matrix estimate is given by

\[  \mc {I}^{-1}(\hat{\bbeta })\sum _{i=1}^ n \biggl [ W_ i(\hat{\bbeta }) W’_ i(\hat{\bbeta }) \biggl ] \mc {I}^{-1}(\hat{\bbeta })  \]


\begin{eqnarray*}  W_ i(\bbeta ) & =&  \sum _ k \Delta _{ki}\biggl \{  Z_ i(T_{ki}) - \bar{\bZ }(\bbeta ,T_{ki}) \biggr \}  - \\ & &  \sum _{i=1}^ n \sum _ l \frac{\Delta _{lj}Y_ i(T_{lj}) \mr {e}^{\bbeta \bZ _ i(T_{lj})}}{S^{0}(\bbeta ,T_{lj}) } \biggl \{  Z_ i(T_{lj}) - \bar{\bZ }(\bbeta ,T_{lj}) \biggr \}  \end{eqnarray*}

For a given realization of the covariates $\bxi $, the Nelson estimator is used to predict the mean function

\[  \hat{\mu }_{\bxi }(t) = \mr {e}^{\hat{\bbeta }\bxi } \sum _{i=1}^ n \sum _ k \frac{I(T_{ki} \le t) \Delta _{ki}}{S^{(0)}(\hat{\bbeta },T_{ki})}  \]

with standard error estimate given by

\[  \hat{\sigma }^2(\hat{\mu }_{\bxi }(t)) = \sum _{i=1}^ n \biggl (\frac{1}{n}\hat{\Psi }_ i(t,\bxi ) \biggr )^2  \]


\begin{eqnarray*}  \frac{1}{n}\hat{\Psi }_ i(t,\bxi ) & = &  \mr {e}^{\hat{\bbeta }\bxi } \biggl \{  \sum _ k \frac{I(T_{ki}\le t)\Delta _{ik}}{S^{(0)}(\hat{\bbeta },T_{ki})} - \sum _{j=1}^ n\sum _ k \frac{Y_ i(T_{kj}) \mr {e}^{\hat{\bbeta }\bZ _ i(T_{kj})}I(T_{kj} \le t)\Delta _{kj}}{[S^{(0)}(\hat{\bbeta },T_{kj})]^2} - \\ & &  \biggl [ \sum _{i=1}^ n \sum _ k \frac{I(T_{ki}\le t)\Delta _{ik} [\bar{\bZ }(\hat{\bbeta },T_{ki}) - \bxi ]}{S^{(0)}(\hat{\bbeta },T_{ki})} \biggr ] \\ & &  \times \mc {I}^{-1}(\hat{\bbeta })\int _0^{\tau } [\bZ _ i(u) - \bar{\bZ }(\hat{\bbeta },u)]d\hat{M}_ i(u)\biggl \}  \end{eqnarray*}

Since the cumulative mean function is always nonnegative, the log transform is used to compute confidence intervals. The $100(1-\alpha )$% pointwise confidence limits for $ \mu _{\bxi }(t)$ are

\[  \hat{\mu }_{\bxi }(t) \mr {e}^{\pm z_{\alpha /2} \frac{\hat{\sigma }(\hat{\mu }_{\bxi }(t)) }{\hat{\mu }{\bxi }(t)}}  \]

where $z_{\alpha /2}$ is the upper $100\alpha /2$ percentage point of the standard normal distribution.