An improved fully nonparametric estimator of the marginal survival function based on case-control clustered data

TL;DR

This paper introduces a new nonparametric estimator for the marginal survival function in case-control family studies, accounting for within-family correlations and right-censored data, with proven asymptotic properties and demonstrated effectiveness through simulations and real data application.

Contribution

It proposes a novel local linear nonparametric estimator for the marginal survival function in family-based case-control studies, addressing within-family correlation and censoring.

Findings

Estimator performs well in simulations

Asymptotic properties are established

Method applied successfully to prostate cancer data

Abstract

A case-control family study is a study where individuals with a disease of interest (case probands) and individuals without the disease (control probands) are randomly sampled from a well-defined population. Possibly right-censored age at onset and disease status are observed for both probands and their relatives. Correlation among the outcomes within a family is induced by factors such as inherited genetic susceptibility, shared environment, and common behavior patterns. For this setting, we present a nonparametric estimator of the marginal survival function, based on local linear estimation of conditional survival functions. Asymptotic theory for the estimator is provided, and simulation results are presented showing that the method performs well. The method is illustrated on data from a prostate cancer study. Keywords: case-control; family study; multivariate survival; nonparametric estimator; local linear

Tables

Table 1. Table 1: Prostate cancer case-control family data

t	The proposed Estimator	Gorfine et al.	Naive KM	SEER
50	0.9997 (0.0007)	0.9918 (0.0311)	0.9991 (0.0003)	0.9958
52	0.9997 (0.0010)	0.9801 (0.0340)	0.9986 (0.0005)	0.9930
54	0.9993 (0.0012)	0.9784 (0.0413)	0.9981 (0.0006)	0.9902
56	0.9990 (0.0023)	0.9784 (0.0451)	0.9963 (0.0008)	0.9843
58	0.9910 (0.0037)	0.9784 (0.0451)	0.9945 (0.0010)	0.9783
60	0.9934 (0.0054)	0.9784 (0.0461)	0.9881 (0.0015)	0.9703
62	0.9908 (0.0063)	0.9678 (0.0501)	0.9848 (0.0017)	0.9603
64	0.9895 (0.0085)	0.9423 (0.0577)	0.9813 (0.0019)	0.9504

Equations

$$S(t)=\text{mbox}pr(T_P>t)=\text{mbox}pr(T_R>t)$$

$$\text{mbox}pr(T_P>t,T_R>u)={\int }_t^{\infty }\text{mbox}pr(T_R>u|T_P=x)f(x)dx$$

$$\Rightarrow S(t)S_0(u|t)={\int }_t^{\infty }{S}_1(u|x)f(x)dx$$

$$\Rightarrow \frac{\partial }{\partial t}[S(t)S_0(u|t)]=-S_1(u|t)f(t)$$

$$\Rightarrow -f(t)S_0(u|t)+S(t)S_0^{\ast }(u|t)=-S_1(u|t)f(t)$$

$$\Rightarrow -\lambda (t)S_0(u|t)+S_0^{\ast }(u|t)=-S_1(u|t)\lambda (t)$$

$$\Rightarrow \lambda (t)(S_0(u|t)-S_1(u|t))=S_0^{\ast }(u|t)=-S_0(u|t){\Lambda }_0^{\ast }(u|t)$$

$$\Rightarrow \lambda (t)(S_0(u|t)-S_1(u|t){)}^2=-S_0(u|t)(S_0(u|t)-S_1(u|t)){\Lambda }_0^{\ast }(u|t)$$

$$\psi (u,t)={\left[{\int }_0^{\tau }(S_0(v|t)-S_1(v|t){)}^{2}dv\right]}^{-1}(S_0(u|t)-S_1(u|t))S_0(u|t)$$

$$\lambda (t)=-{\int }_0^{\tau }\psi (u,t){\Lambda }_0^{\ast }(u|t)du$$

$$\Lambda (t)=-{\int }_0^t{\int }_0^{\tau }\psi (u,s){\Lambda }_0^{\ast }(u|s)duds$$

$$\displaystyle\lambda_{q}(u|t)=\frac{\partial}{\partial u}\Lambda_{q}(u|t),\quad N_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)=\sum_{j=1}^{J}N_{Rij}(v),\quad Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)=\sum_{j=1}^{J}Y_{Rij}(v)$$

$$\displaystyle dM_{Rij}(v)=dN_{Rij}(v)-Y_{Rij}(v)\lambda_{Q_{i}}(v|X_{Pi})dv,\quad M_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)=\sum_{j=1}^{J}M_{Rij}(v)$$

$$\displaystyle\mathcal{Y}_{q}(s,v)=\frac{1}{n_{q}h}\sum_{i=1}^{n}\chi(Q_{i},q)Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K\left(\frac{X_{Pi}-s}{h}\right)$$

$$\displaystyle\bar{X}_{Pq}(s,v)=\mathcal{Y}_{q}(s,v)^{-1}\left[\frac{1}{n_{q}h}\sum_{i=1}^{n}\chi(Q_{i},q)Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K\left(\frac{X_{Pi}-s}{h}\right)X_{Pi}\right]$$

$$\displaystyle C_{q}(s,v)=\frac{1}{n_{q}h}\sum_{i=1}^{n}\chi(Q_{i},q)Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K\left(\frac{X_{Pi}-s}{h}\right)(X_{Pi}-\bar{X}_{Pq}(s,v))^{2}$$

$$\displaystyle\frac{dN_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)}{Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)}$$

$${\hat{\Lambda }}_q^{\ast }(dv|s)$$

$${\hat{\Lambda }}_q(dv|s)$$

$${\hat{\Lambda }}_q^{\ast }(u|s)$$

$${\hat{\Lambda }}_q(u|s)$$

$$\displaystyle-\int_{0}^{u}\frac{\bar{X}_{Pq}(s,v)-s}{C_{q}(s,v)}\left[\frac{1}{n_{q}h}\sum_{i}\chi(Q_{i},q)K\left(\frac{X_{Pi}-s}{h}\right)(X_{Pi}-\bar{X}_{Pq}(s,v))dN_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)\right]$$

$$\mathcal{A}(t)={\int }_0^t{\int }_0^{\tau }\psi (u,s)[{\hat{\Lambda }}_0^{\ast }(u|s)-{\Lambda }_0^{\ast }(u|s)]duds$$

$${\mathcal{A}}^{\prime }(t)={\int }_0^t{\int }_0^{\tau }[\hat{\psi }(u,s)-\psi (u,s))]{\Lambda }_0^{\ast }(u|s)duds$$

$${\mathcal{A}}^{\prime \prime }(t)={\int }_0^t{\int }_0^{\tau }[\hat{\psi }(u,s)-\psi (u,s)][{\hat{\Lambda }}_0^{\ast }(u|s)-{\Lambda }_0^{\ast }(u|s)]duds$$

$$\Gamma(s,v)=h^{-2}C_{0}(s,v)=\frac{1}{mh}\sum_{i=1}^{m}Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K\left(\frac{X_{Pi}-s}{h}\right)\left(\frac{X_{Pi}-\bar{X}_{P0}(s,v)}{h}\right)^{2}$$

$$\mathcal{A}(t)=-\frac{1}{m{h}^3}\sum _{i=1}^m{\int }_0^t{\int }_0^{\tau }{\int }_0^u$$

$$\displaystyle(dN_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)-Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)(X_{Pi}-\bar{X}_{P0}(s,v))\lambda_{0}^{*}(v|s)dv)\,du\,ds$$

$${\int }_0^u\mathcal{P}(v)d{M}_{Rij}(v)$$

$$E\left[\int_{0}^{u}P(v,X_{Pi})dM_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)\right]=0$$

$$\displaystyle\mathcal{A}_{1}(t)=\frac{1}{mh^{3}}\sum_{i=1}^{m}\int_{0}^{t}\int_{0}^{\tau}\int_{0}^{u}\psi(u,s)\Gamma(s,v)^{-1}Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K\left(\frac{X_{Pi}-s}{h}\right)(X_{Pi}-\bar{X}_{P0}(s,v))$$

$$({\lambda }_0(v|X_{Pi})-(X_{Pi}-{\bar{X}}_{P0}(s,v)){\lambda }_0^{\ast }(v|s))dvduds$$

$$\displaystyle\mathcal{A}_{2}(t)=\frac{1}{mh^{3}}\sum_{i=1}^{m}\int_{0}^{t}\int_{0}^{\tau}\int_{0}^{u}\psi(u,s)\Gamma(s,v)^{-1}K\left(\frac{X_{Pi}-s}{h}\right)(X_{Pi}-\bar{X}_{P0}(s,v))dM_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)\,du\,ds$$

$${\mathcal{A}}_2(t)=\frac{1}{m}\sum _{i=1}^m{\zeta }_i(t)+o_p((mh{)}^{-1/2})$$

$$\zeta_{i}(t)=\frac{1}{h}\int_{0}^{t}\int_{0}^{\tau}\int_{0}^{u}\psi(u,s)\gamma(s,v)^{-1}\left(\frac{X_{Pi}-s}{h}\right)\left[h^{-1}K\left(\frac{X_{Pi}-s}{h}\right)\right]dM_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)\,du\,ds$$

$$\text{mbox}pr(T_R>u|T_P\le C_P)\le S(u)\le \text{mbox}pr(T_R>u|T_P>C_P)$$

$$\underset{t\in \mathcal{T}(u)}{\inf }P(T_P>t,T_R>u)-P(T_P>t)P(T_R>u)>0$$

$$\bar{S}(t;h)$$

$$\text{mbox}{mse}_{est}(t,h)$$

Full text

Keywords: case-control; family study; multivariate survival; nonparametric estimator; local linear

1 Introduction

Many epidemiological and medical studies focus on disease events that are rare, so that a random sample from the population provides very few observations where failure has occurred during the monitoring time. In such situations, it is a common practice to consider a case-control strategy. Separate samples of $n_{1}$ individuals in whom the event has already occurred (case probands) and $n_{0}$ individuals in whom the event has not yet occurred (control probands) are obtained. Age at onset or age at censoring and disease status of each proband and of one or more of his/her relatives are recorded. For example, in a prostate cancer study, case probands are men diagnosed with prostate cancer, control probands are men without prostate cancer, and each proband is interviewed to obtain detailed disease history information of his relatives. The goals of such case-control family studies are to evaluate the effect of genetic and environmental factors on disease risk, and to estimate the distribution of the age at onset. The analysis of such data is complicated by the case-control selection scheme used to ascertain the families and by the within-family dependence. It is of interest to estimate the marginal survival function under dependent failure times of family members by a fully nonparametric estimator which avoids specific assumptions about the form of the distribution or the dependence structure among failure times with a family.

The estimator that first comes to mind is the Kaplan-Meier survival curve estimator based on the survival data on the relatives. This estimator, however, is biased because it does not take the case-control sampling and the within-family dependence into account. Gorfine et al. (2017) demonstrated the serious bias that can arise with this naive Kaplan-Meier estimator.

Accordingly, Gorfine et al. proposed a nonparametric estimator for this problem, based on a kernel smoothing approach. They presented a simulation study showing that their estimator performs well in terms of bias. The estimator of Gorfine et al. is based on the median of random variables with a complicated dependence structure. Consequently, the asymptotic properties of the estimator were not derived. Also, their use of the median causes some efficiency loss. In addition, their bandwidth selection procedure was not specifically targeted to the estimand of interest.

In the present paper we develop a new nonparametric estimator which performs well in terms of bias and much better in terms of variance than the estimator of Gorfine et al. In some scenarios the new estimator also outperforms Gorfine et al.  in terms of bias. We provide asymptotic theory for the estimator. We also present a bandwidth selection procedure that is specifically targeted to the estimand of interest.

2 Preliminaries

The setup is as in Gorfine et al. (2017). We denote the maximum observation time among the probands by $\tau_{0}$ and the maximum observation time among the relatives by $\tau$. We write $n=n_{0}+n_{1}$. We let $J_{i}$ denote the number of relatives in family $i$ and we let $J$ denote the maximum number of relatives for a given proband. We view $J_{i}$ as a random variable. We will write some of the formulas as if each proband has exactly $J$ relatives, with the extra relatives taken to be censored at time 0. For family $i$, $i=1,\ldots,n$, let $T_{Pi}$ denote the failure time of the proband and $T_{Rij},j=1,\ldots,J_{i}$ the failure times of the relatives. Let $C_{Pi}$ and $C_{Rij}$ denote the corresponding censoring times, $X_{Pi}=\min(T_{Pi},C_{Pi})$ and $X_{Rij}=\min(T_{Rij},C_{Pij})$ the corresponding observed times, and $\delta_{Pi}=I(T_{Pi}\leq C_{Pi})$ and $\delta_{Rij}=I(T_{Rij}\leq C_{Rij})$ the corresponding event indicators. We assume that the censoring is independent of the survival times. The sample includes $n_{1}$ case probands, with each case proband frequency matched with $a$ control probands, so that the total number of control probands is $n_{0}=an_{1}$. Thus, the data consist of $n_{1}$ independent and identically distributed matched sets comprising one case family and $a$ control families, and the observed data on family $i$ consists of $(X_{Pi},\delta_{Pi},X_{Ri1},\ldots,X_{RiJ_{i}},\delta_{Ri1},\ldots,\delta_{RiJ_{i}})$. We seek to estimate the marginal survival function

[TABLE]

which we are assuming is the same for probands and relatives, a common assumption in case-control family studies (Shih & Chatterjee, 2002; Chatterjee et al., 2006). We further assume that the bivariate survival function for the proband and a given relative is the same for all relatives. The marginal survival distribution is assumed to be continuous with density $f(t)$. In the case where $T_{P}$ and $T_{R}$ have different marginal distributions, our procedure yields a consistent estimate of the marginal distribution of $T_{P}$.

Denote $S_{0}(u|t)=\mbox{pr}(T_{R}>u|T_{P}>t)$ and $S_{1}(u|t)=\mbox{pr}(T_{R}>u|T_{P}=t)$. As in Gorfine et al., we have $\mbox{pr}(T_{R}>u|X_{P}=t,\delta_{P}=0)=S_{0}(u|t)$ and $\mbox{pr}(T_{R}>u|X_{P}=t,\delta_{P}=0)=S_{1}(u|t)$.

We now develop an expression for $S(t)$ in terms of $S_{0}(u|t)$ and $S_{1}(u|t)$. Let $\lambda(t)$ and $\Lambda(t)$ denote the hazard and cumulative hazard functions corresponding to $S(t)$, and define $S_{q}^{*}(u|t)=(\partial/\partial t)S_{q}(u|t)$, $q=0,1$. Also define $\Lambda_{q}(u|t)=-\log S_{q}(u|t)$ and $\Lambda_{q}^{*}(u|t)=(\partial/\partial t)\Lambda_{q}(u|t)$, $q=0,1$. We can write $S_{0}^{*}(u|t)=-S_{0}(u|t)\Lambda_{0}^{*}(u|t)$. We then have the following:

[TABLE]

Now define

[TABLE]

Note that the bracketed integral is nonzero provided that for every $t$ there exists a set of $u$ values of positive measure for which $S_{0}(u|t)\neq S_{1}(u|t)$. Integrating both sides of (1) and rearranging gives

[TABLE]

We use (4) to construct our estimator.

3 Estimation Procedure

In Gorfine et al., $S_{0}(u|t)$ and $S_{1}(u|t)$ were estimated using a generalized version of the kernel-smoothed Kaplan-Meier estimator proposed in an unpublished 1981 University of California at Berkeley technical report by R. Beran and examined in Dabrowska (1987), and the resulting estimators were used to construct an estimator of $S(t)$. Here, in light of (4), we work not only with $S_{0}(u|t)$ and $S_{1}(u|t)$ but also the derivative $\Lambda_{0}^{*}(u|t)$. Accordingly, we take a local linear estimation approach. Choose a symmetric kernel function $K$ and a bandwidth $h$. Let $N_{Rij}(t)=\delta_{Rij}I(T_{Rij}\leq t)$ and $Y_{Rij}(t)=I(T_{Rij}\geq t)$, and write $Q_{i}=\delta_{Pi}$. Let $\chi(q_{1},q_{2})=I(q_{1}=q_{2})$. Define (with $q=0,1$)

[TABLE]

We can write

[TABLE]

with E[dM_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)/Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)]=0. A first-order Taylor approximation gives the local linear representation

[TABLE]

We now carry out weighted linear least squares with response variable dN_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)/Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v), explanatory variable $X_{Pi}-s$, and weights \chi(Q_{i},q)Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)K((X_{Pi}-s)/h). This leads to the local linear estimators

[TABLE]

that is, for a given $u$,

[TABLE]

with the second equation leading to $\widehat{S}_{q}(u|s)=\exp(-\widehat{\Lambda}_{q}(u|s))$. We can now substitute $\widehat{S}_{0}(u|s)$, $\widehat{S}_{0}(u|s)$, and $\widehat{\Lambda}_{0}^{*}(u|s)$ into (3) and (4) to obtain estimators $\widehat{\lambda}(t)$ and $\widehat{\Lambda}(t)$ for $\lambda(t)$ and $\Lambda(t)$. We then take $\widehat{S}(t)=\exp(-\widehat{\Lambda}(t))$. When we want to emphasize the dependence on the bandwidth $h$, we will write $\widehat{S}(t;h)$.

4 Asymptotic Theory

We can write $\widehat{\Lambda}(t)-\Lambda(t)=\mathcal{A}(t)+\mathcal{A}^{\prime}(t)+\mathcal{A}^{\prime\prime}(t)$ with

[TABLE]

We will provide a detailed asymptotic analysis of $\mathcal{A}(t)$. Similar arguments can be used to show that $\mathcal{A}^{\prime}(t)$ converges in probability to zero at faster rate than $\mathcal{A}(t)$, and $\mathcal{A}^{\prime\prime}(t)$ is negligible in comparison with the other two terms.

In this section and in the appendix, we will write $m=n_{0}$ and assume that the indices have been arranged so that the control probands appear first in the list of probands, meaning that a sum over probands 1 to $m$ is a sum over the control probands.

Defining

[TABLE]

we have

[TABLE]

Now, the process $M_{Rij}$ is a martingale with respect to the filtration $\mathcal{F}_{ijv}=\sigma(X_{Pi},$ $\{N_{Rij}(d),$ $Y_{Rij}(d)$, $d\in[0,v]\})$. Accordingly, for any process $\mathcal{P}(v)$ that is predictable with respect to $\mathcal{F}_{ijv}$, the process

[TABLE]

is a mean-zero martingale. It follows, even though the process M_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ } does not have any martingale properties, that for any function $P(v,x)$ we have

[TABLE]

We now write $\mathcal{A}(t)=-(\mathcal{A}_{1}(t)+\mathcal{A}_{2}(t))$, where

[TABLE]

In the appendix it is shown, under the conditions listed in the appended supplementary document that $(mh)^{1/2}\mathcal{A}_{1}(t)\rightarrow 0$ in probability and that

[TABLE]

with

[TABLE]

where $\gamma$ is the limiting value of $\Gamma$. Define $\sigma_{\zeta}^{2}(t)=h\,\mbox{\rm Var}(\zeta_{i}(t))$. We then obtain the following theorem.

Theorem 1**.**

For each $t\in[0,\tau_{0}]$, $\sigma_{\zeta}^{2}(t)$ converges to a limit $\kappa(t)$ and $(mh)^{1/2}\,(\widehat{\Lambda}(t)-\Lambda(t))$ converges in distribution to $N(0,\kappa(t))$.

The proof of this theorem appears in the appendix. Details are given in the appended supplementary document.

5 Practical Implementation Details

In preliminary work, we found that the performance of the estimator of $\widehat{S}(t)$ can be improved dramatically by introducing a time transformation that makes the proband observation times approximately uniformly distributed. Along the lines of Doksum et al. (2017), we propose transforming according to an estimate of the distribution function of the proband observation times, which leads to a modified form of nearest neighbor regression. In the appended supplementary document, we show that the consistency and asymptotic normality is maintained under the time transformation if a smooth estimate of the distribution function is used. We believe that this result holds as well when the empirical distribution function is used. In our numerical work, we used the empirical distribution function.

In the context of family survival data, it is usually reasonable to assume that $\mbox{pr}(T_{P}>t,T_{R}>u)\geq\mbox{pr}(T_{P}>t)\mbox{pr}(T_{R}>u)=S(t)S(u)$ for all $t$ and $u$, i.e., $S_{0}(u|t)\geq S(u)$ for all $t$ and $u$. This condition implies that

[TABLE]

If we let $\widehat{S}_{KM,case}(t)$ and $\widehat{S}_{KM,control}(t)$ denote the Kaplan-Meier survival curve estimator based on the case relatives’ survival data and the control relatives’ survival data, respectively, the foregoing inequalities motivate modifying the estimator to the estimator $\tilde{S}(t)$ resulting from replacing $\widehat{S}(t)$ with $\widehat{S}_{KM,case}(t)$ if $\widehat{S}(t)\leq\widehat{S}_{KM,case}(t)$ and by $\widehat{S}_{KM,control}(t)$ if $\widehat{S}(t)\geq\widehat{S}_{KM,control}(t)$. We implemented this modification in our numerical work. The modification comes into play mainly when the event rate is extremely low or extremely high. Given the consistency of $\widehat{S}(t)$, if the inequalities in (8) are strict, then for large sample sizes the modification no longer comes into play. The inequalities in (8) are strict if the following mild condition holds: for each $u$ there exists a set $\mathcal{T}(u)$ such that $P(C_{P}\in\mathcal{T}(u))>0$ and

[TABLE]

For bandwidth selection, we propose a bootstrap procedure. Let $\widehat{S}^{C}(u)$ denote the Kaplan-Meier estimate of the survival function of the time to censoring among the relatives (which is the same for case relatives and control relatives). In each bootstrap replication $b=1,\ldots,B$, for each family $i$ we generate $J_{i}$ event times for proband $i$’s relatives according to the survival function $\widehat{S}_{Q_{i}}(u|T_{Pi})$ and $J_{i}$ censoring times for relatives according to the survival function $\widehat{S}^{C}(u)$. We then run our estimation procedure for a given $h$ on the resulting data, obtaining the estimate $\widehat{S}(u;h,b)$. Denote

[TABLE]

We evaluate $\mbox{{imse}}_{est}(h)$ over a grid of $h$ values in the range $(0,1]$ and choose the $h$ values with the minimum $\mbox{{imse}}_{est}(h)$.

To construct confidence intervals in finite samples with bandwidth selection, we use the percentile bootstrap method.

6 Simulation Study

We carried out a simulation study to evaluate the finite sample properties of the proposed estimator. Data were generated under frailty models in which the within-family dependence is expressed in terms of a shared frailty variate $W_{i}$, conditional on which the failure times of the family members are independent with hazard function $\lambda(t|W_{i})=W_{i}\lambda_{0}(t)$. We manipulated five design factors, as follows: (1) frailty distribution: gamma or positive stable, (2) cumulative end-of-study event rate: high (60%) or low (15%), (3) number of case probands: 500 or 1000 (with 1:1 matching of control probands to case probands), (4) number of relatives per family: 1 or 4, and (5) strength of within-family dependence: low (Kendall tau of 1/3 between the failure of times of two members of the same family) and moderate (Kendall tau of 1/2). We took $\lambda_{0}(t)=\nu(\mu t)^{p-1}$ with $p=4.6$, $\mu=0.01$, and $\nu$ chosen so as to obtain the desired cumulative end-of-study event rate. The end of study age was taken to be 110 years. The overall censoring rate, including both interim and end-of-study censoring, was about 60% in the high event rate case and 90% in the low event rate case. The number of case probands was taken to be 500, with 1:1 matching of control probands to case probands. The data generation was carried out in the same manner as in Gorfine et al. We carried out 1024 simulation replications for each of the 32 combinations of the design. For each replication, we carried out 30 inner replications for the bootstrap bandwidth selection procedure and 100 outer replications for the percentile bootstrap confidence interval procedure. The initial bandwidth was 0.5 and the bandwidth search was done in two stages. In the first stage, we searched over $[0.1,1]$ in steps of 0.1 and identified the $h$ value $h_{1}$ with the lowest $\mbox{{imse}}_{est}(h)$. In the second stage, we searched over $h_{1}-0.05,h_{1}$, and $\min(h_{1}+0.05,1)$ and chose the $h$ value with the lowest $\mbox{{imse}}_{est}(h)$ to be the final $h$ value. The kernel used was the triweight kernel $K(u)=(35/32)I(|u|\leq 1)(1-u^{2})^{3}$.

The results for 500 case probands are summarized in Fig. 1 and 2. The left two columns of each figure show the true survival curve, along with Gorfine et al.’s estimator and the new estimator. The finite-sample bias of the new estimator tends to be smaller, and in some settings, such as the gamma frailty model with very low event rates, its finite-sample bias is dramatically smaller. The right two columns of each figure summarize the point-wise $95\%$ coverage rates of the percentile-bootstrap confidence interval of the proposed estimator, along with the standard errors of the Gorfine et al. estimator and the proposed estimators. Clearly, the proposed estimator substantially outperforms the old estimator in terms of efficiency. In general, the coverage rates are reasonably close to $95\%$, except at very early ages with small number of observed events. Similar results were obtained with 1000 case probands.

7 Example

In this section we illustrate our method by re-analysing the data presented in Gorfine et al. (2017), population-based case-control family study of early onset prostate cancer (Stanford et al., 1999). Briefly, case participants were identified from the Seattle-Puget Sound Surveillance, Epidemiology, and End Results (SEER) cancer registry. Cases were those with age at diagnosis between 40 and 64 years. Controls were identified by use of random-digit dialing and they were frequency matched to case participants by age. The information collected on the relatives is the age at diagnosis for prostate cancer if the relative had prostate cancer or age at the last observation if the relative did not have prostate cancer. Here we use the information about age at onset or age at censoring and disease status that was observed for the probands and their fathers, brothers, and uncles. The following analysis is based on 730 prostate-cancer case probands, 693 control probands, and a total of 7316 relatives. Out of the 3793 case-probands’ relatives, 211 had prostate cancer, and out of the 3523 control-probands’ relatives, 102 had prostate cancer. The age range of the relatives with prostate cancer was 40–93. Given that frequency matching was used rather than exact age matching, and that the number of relatives per proband varied across the probands, we carried out the time transformation based on the empirical distribution of the proband observation times across all 7316 relatives in the data set. For bandwidth selection we used the same two-stage procedure as in the simulations.

Figure 3 and Table 1 present the estimates of prostate-cancer marginal survival function using the naive Kaplan-Meier estimator based on the relatives’ data, Gorfine et al.’s estimator with nearest-neighbor smoothing and the median operator, the SEER survival curve based on the SEER Cancer Statistics Review 1975–2012, and the proposed estimator. In this dataset, Gorfine et al.’s estimator is closer to the SEER survival curve, but with very large point-wise standard errors compared to the proposed estimator. The standard errors reported in Table 1 are much larger than those reported in Gorfine et al. due to an error in the bootstrap code applied back then.

Acknowledgements

The authors would like to thank Dr. Stanford for her generosity in sharing the prostate cancer dataset that was used for illustrating the method. Malka Gorfine gratefully acknowledges support from the U.S.-Israel Binational Science Foundation in carrying out this work.

Supplementary material

We append a document with details of the proof of the asymptotic properties of the estimator. R code used to carry out the simulations and R code for applying the method to a data set may be found at https://github.com/david-zucker/marginal-survival/.

Appendix: Asymptotic Theory

A. Preliminaries

We first present some definitions. Let $g(x)$ denote the density of $X_{Pi}$. Define y(s,v)=E[Y_{Ri\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ }(v)|X_{Pi}=s,\delta_{Pi}=0] and $\varphi(s,v)=g(s)y(s,v)$. In addition, define

[TABLE]

Note that, by symmetry of $K$, $\mu_{k}(1)=0$ for all odd $k$.

We now present two lemmas, whose proofs appear in the online Supplemental Materials. The notations $O$ and $o$, and similarly $O_{p}$ and $o_{p}$, should be understood as being uniform in the relevant arguments.

Lemma 1**.**

For $k$ even we have

[TABLE]

and for $k$ odd we have

[TABLE]

Lemma 2**.**

For $s\in\mathcal{I}$ we have

[TABLE]

and for $s\in\mathcal{U}$ we have

[TABLE]

B. Analysis of $\mathcal{A}_{1}(t)$

We can write $\mathcal{A}_{1}(t)$ as

[TABLE]

with

[TABLE]

By Taylor expansion, we can write

[TABLE]

in which $|\mathcal{R}(s,v,x)|\leq\mathcal{R}^{*}|s-x|^{3}$, with $\mathcal{R}^{*}=\sup_{s,v}|\lambda_{0}^{***}(v|s)|/6$, where $\lambda_{0}^{**}(v|s)$ and $\lambda_{0}^{***}(v|s)$ denote, respectively, the second and third partial derivatives of $\lambda_{0}(v|s)$ with respect to $s$ We then have

[TABLE]

The term in square brackets does not depend on $i$. Since

[TABLE]

we get

[TABLE]

We can then write $\mathcal{S}(s,v)=\mathcal{S}_{1}(s,v)+\mathcal{S}_{2}(s,v)$, where

[TABLE]

We have

[TABLE]

We now consider separately the case of $s\in\mathcal{I}$ and $s\in\mathcal{U}$. For $s\in\mathcal{I}$, the results of Lemmas 1 and 2 imply that $\mathcal{S}_{1}(s,v)=O_{p}(h^{2})$ and $\mathcal{S}_{2}(s,v)=O_{p}(h^{2})$, so that $\mathcal{S}(s,v)=O_{p}(h^{2})$ and

[TABLE]

For $s\in\mathcal{U}$, the results of Lemmas A1 and A2 imply that $\mathcal{S}_{1}(s,v)=O_{p}(h)$ and $\mathcal{S}_{2}(s,v)=O_{p}(h)$, so that $\mathcal{S}(s,v)=O_{p}(h)$ and

[TABLE]

(recalling that the length of $\mathcal{U}$ is 2$h$). We thus obtain $\mathcal{A}_{1}=O_{p}(h^{2})$, so that $(mh)^{1/2}\mathcal{A}_{1}(t)=o_{p}(1)$, since $\nu>1/4$.

C. Analysis of $\mathcal{A}_{2}(t)$

Let $\rho(s,u,v)=\psi(u|s)/\varphi(s,v)$ and define

[TABLE]

Define further

[TABLE]

We can then write $\mathcal{A}_{2}(t)=\mathcal{B}(t)+\mathcal{B}_{1}(t)+\mathcal{B}_{2}(t)$, where

[TABLE]

Our claim is that $(mh)^{1/2}\mathcal{B}(t)$ is asymptotically mean-zero normal, and that $(mh)^{1/2}\mathcal{B}_{1}(t)$ and $(mh)^{1/2}\mathcal{B}_{2}(t)$ are both $o_{p}(1)$.

By (7), $E[\zeta_{i}(t)]=0$. Thus, $\mathcal{B}$ is the sum of i.i.d. mean-zero terms. Accordingly, to show that $(mh)^{1/2}\,\mathcal{B}(t)$ is asymptotically mean-zero normal we need to show that $h\,\mbox{\rm Var}(\zeta_{i}(t))$ converges to a limit $\kappa(t)$ and that $\zeta_{i}(t)$ satisfies the Lindeberg condition

[TABLE]

where

[TABLE]

The appended supplementary document provides a proof of the above two assertions, along with a proof that $(mh)^{1/2}\mathcal{B}_{1}(t)$ is $o_{p}(1)$. The proof of the corresponding result for $\mathcal{B}_{2}(t)$ is similar.