A likelihood-based approach for cure regression models

TL;DR

This paper introduces a likelihood-based method for estimating and selecting variables in cure regression models, effectively handling censored data and cured subjects without assuming specific survival time distributions.

Contribution

It presents a novel IPCW-based likelihood estimation framework with integrated variable selection for cure models, applicable under minimal assumptions.

Findings

Effective variable selection with L1 penalties.

Accurate estimation demonstrated via simulations.

Successful application to real data.

Abstract

We propose a new likelihood-based approach for estimation, inference and variable selection for parametric cure regression models in time-to-event analysis under random right-censoring. In this context, it often happens that some subjects are "cured", i.e., they will never experience the event of interest. Then, the sample of censored observations is an unlabeled mixture of cured and "susceptible" subjects. Using inverse probability censoring weighting (IPCW), we propose a likelihood-based estimation procedure for the cure regression model without making assumptions about the distribution of survival times for the susceptible subjects. The IPCW approach does require a preliminary estimate of the censoring distribution, for which general parametric, semi- or non-parametric approaches can be used. The incorporation of a penalty term in our estimation procedure is straightforward; in particular, we propose L1-type penalties for variable selection. Our theoretical results are derived under mild assumptions. Simulation experiments and real data analysis illustrate the effectiveness of the new approach.

Tables

Table 1. Table 1: Average bias and standard error (in brackets) of estimates ( π m = 0.4 , ρ = 0.1 formulae-sequence subscript 𝜋 𝑚 0.4 𝜌 0.1 \pi_{m}=0.4,\rho=0.1 )

		$n=100$			$n=300$			$n=1000$
Method	$\nu$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
Our	0	-0.05	0.10	0.12	-0.01	0.04	0.03	0.00	0.01	0.01
proposal		(0.36)	(0.49)	(0.47)	(0.18)	(0.24)	(0.22)	(0.10)	(0.12)	(0.12)
	2	-0.05	0.13	0.11	-0.02	0.04	0.02	0.00	0.01	0.01
		(0.37)	(0.53)	(0.44)	(0.19)	(0.25)	(0.21)	(0.10)	(0.13)	(0.11)
smcure	0	-0.08	0.09	0.09	-0.02	0.03	0.02	-0.01	0.01	0.01
		(0.33)	(0.40)	(0.38)	(0.17)	(0.20)	(0.19)	(0.09)	(0.11)	(0.10)
	2	-0.12	0.10	-0.09	-0.08	0.04	-0.13	-0.06	0.01	-0.13
		(0.35)	(0.40)	(0.33)	(0.19)	(0.21)	(0.17)	(0.10)	(0.11)	(0.09)

Table 2. Table 2: Standard error of estimates using different B 𝐵 B types ( π m = 0.4 , ρ = 0.1 formulae-sequence subscript 𝜋 𝑚 0.4 𝜌 0.1 \pi_{m}=0.4,\rho=0.1 )

		$n=100$			$n=300$			$n=1000$
$\nu$	$B$ Type	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
0	$B^{\ast }\left(S_C\right)$	(0.40)	(0.53)	(0.50)	(0.22)	(0.28)	(0.26)	(0.12)	(0.14)	(0.13)
	$B^{\ast }\left({\hat{S}}_C\right)$	(0.36)	(0.49)	(0.47)	(0.18)	(0.24)	(0.22)	(0.10)	(0.12)	(0.12)
	$B$	(0.25)	(0.32)	(0.30)	(0.14)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)
2	$B^{\ast }\left(S_C\right)$	(0.45)	(0.58)	(0.50)	(0.23)	(0.29)	(0.25)	(0.12)	(0.14)	(0.12)
	$B^{\ast }\left({\hat{S}}_C\right)$	(0.37)	(0.53)	(0.44)	(0.19)	(0.25)	(0.21)	(0.10)	(0.13)	(0.11)
	$B$	(0.26)	(0.30)	(0.30)	(0.14)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)

Table 3. Table 3: Correct zeros, incorrect zeros, and model degrees of freedom ( π m = 0.4 , ρ = 0.1 formulae-sequence subscript 𝜋 𝑚 0.4 𝜌 0.1 \pi_{m}=0.4,\rho=0.1 )

		$n=100$			$n=300$			$n=1000$
Type	$\nu$	C	IC	DF	C	IC	DF	C	IC	DF
oracle		4.00	0.00	3.00	4.00	0.00	3.00	4.00	0.00	3.00
alasso	0	3.34	0.17	3.49	3.52	0.00	3.48	3.69	0.00	3.31
	2	3.33	0.15	3.52	3.55	0.00	3.44	3.72	0.00	3.28

Table 4. Table 4: Colon cancer estimates

		Unpenalized			alasso		smcure
Covariate		Est.	95%CI	pval		Est.	Est.	95%CI	pval
Intercept		0.66	( 0.05, 1.36)	0.03		0.32	0.57	(-0.05, 1.15)	0.06
Treatment	Lev	0.42	(-0.11, 1.22)	0.14		0.00	0.19	(-0.21, 0.61)	0.33
	Lev+5FU	0.94	( 0.39, 1.73)	0.00		0.60	0.71	( 0.30, 1.15)	0.00
Surgery	$>20$days	-0.65	(-1.63,-0.11)	0.02		-0.41	-0.49	(-0.88,-0.09)	0.01
Age	Years	-0.01	(-0.03, 0.00)	0.10		0.00	-0.01	(-0.02, 0.00)	0.24
Sex	Male	-0.24	(-0.75, 0.16)	0.28		0.00	-0.11	(-0.46, 0.25)	0.60
Obstruction	Yes	-0.56	(-2.01, 0.05)	0.08		-0.19	-0.18	(-0.61, 0.20)	0.35
Adherence	Yes	-0.42	(-1.00, 0.07)	0.10		0.00	-0.69	(-1.50,-0.17)	0.01
Depth	Serosa	-0.81	(-1.51,-0.26)	0.01		-0.54	-0.71	(-1.27,-0.19)	0.02
Nodes	$>4$	-1.18	(-1.63,-0.82)	0.00		-0.94	-1.16	(-1.59,-0.80)	0.00
Age is mean-centered.

Table 5. Table 5: Average bias and standard error (in brackets) of estimates

				$n=100$			$n=300$			$n=1000$
Method	$\nu$	${\pi }_m$	$\rho$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
Our	0	0.2	0.1	-0.20	0.14	0.12	-0.06	0.05	0.03	-0.01	0.01	0.01
proposal				(0.65)	(0.62)	(0.58)	(0.29)	(0.29)	(0.26)	(0.15)	(0.14)	(0.14)
			0.2	-0.17	0.12	0.12	-0.06	0.05	0.05	-0.02	0.01	0.02
				(0.80)	(0.74)	(0.68)	(0.40)	(0.40)	(0.34)	(0.20)	(0.20)	(0.17)
		0.4	0.1	-0.05	0.10	0.12	-0.01	0.04	0.03	0.00	0.01	0.01
				(0.36)	(0.49)	(0.47)	(0.18)	(0.24)	(0.22)	(0.10)	(0.12)	(0.12)
			0.2	-0.03	0.16	0.17	-0.02	0.05	0.05	0.00	0.02	0.02
				(0.51)	(0.74)	(0.66)	(0.29)	(0.37)	(0.32)	(0.13)	(0.18)	(0.16)
	2	0.2	0.1	-0.16	0.12	0.11	-0.05	0.04	0.03	-0.01	0.01	0.01
				(0.62)	(0.60)	(0.51)	(0.29)	(0.29)	(0.24)	(0.15)	(0.15)	(0.13)
			0.2	-0.13	0.13	0.10	-0.09	0.07	0.06	-0.02	0.01	0.01
				(0.80)	(0.76)	(0.60)	(0.45)	(0.43)	(0.32)	(0.21)	(0.20)	(0.15)
		0.4	0.1	-0.05	0.13	0.11	-0.02	0.04	0.02	0.00	0.01	0.01
				(0.37)	(0.53)	(0.44)	(0.19)	(0.25)	(0.21)	(0.10)	(0.13)	(0.11)
			0.2	-0.03	0.15	0.15	-0.02	0.06	0.05	-0.01	0.02	0.01
				(0.55)	(0.77)	(0.60)	(0.30)	(0.41)	(0.27)	(0.15)	(0.20)	(0.13)
smcure	0	0.2	0.1	-0.20	0.12	0.08	-0.05	0.04	0.02	-0.02	0.01	0.01
				(0.54)	(0.48)	(0.46)	(0.25)	(0.23)	(0.22)	(0.14)	(0.12)	(0.12)
			0.2	-0.40	0.20	0.14	-0.09	0.04	0.05	-0.03	0.02	0.02
				(0.93)	(0.72)	(0.68)	(0.34)	(0.30)	(0.27)	(0.17)	(0.15)	(0.14)
		0.4	0.1	-0.08	0.09	0.09	-0.02	0.03	0.02	-0.01	0.01	0.01
				(0.33)	(0.40)	(0.38)	(0.17)	(0.20)	(0.19)	(0.09)	(0.11)	(0.10)
			0.2	-0.25	0.20	0.12	-0.06	0.05	0.03	-0.01	0.02	0.01
				(0.64)	(0.68)	(0.61)	(0.26)	(0.27)	(0.25)	(0.12)	(0.14)	(0.13)
	2	0.2	0.1	-0.16	0.11	-0.03	-0.03	0.02	-0.07	0.00	0.00	-0.09
				(0.54)	(0.48)	(0.42)	(0.25)	(0.23)	(0.20)	(0.13)	(0.12)	(0.11)
			0.2	-0.18	0.12	-0.10	-0.02	0.02	-0.14	0.04	-0.02	-0.16
				(0.73)	(0.59)	(0.48)	(0.33)	(0.27)	(0.23)	(0.17)	(0.14)	(0.11)
		0.4	0.1	-0.12	0.10	-0.09	-0.08	0.04	-0.13	-0.06	0.01	-0.13
				(0.35)	(0.40)	(0.33)	(0.19)	(0.21)	(0.17)	(0.10)	(0.11)	(0.09)
			0.2	-0.22	0.15	-0.17	-0.12	0.05	-0.20	-0.09	0.01	-0.20
				(0.55)	(0.57)	(0.42)	(0.28)	(0.27)	(0.20)	(0.14)	(0.13)	(0.10)

Table 6. Table 6: Absolute average bias and standard error (in brackets) of estimates: Six covariates

				$n=100$			$n=300$			$n=1000$
Method	$\nu$	${\pi }_m$	$\rho$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\text{-}2\right)}$	${\theta }_{\left(3\text{-}6\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\text{-}2\right)}$	${\theta }_{\left(3\text{-}6\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\text{-}2\right)}$	${\theta }_{\left(3\text{-}6\right)}$
Our	0	0.2	0.1	0.46	0.27	0.01	0.13	0.08	0.00	0.04	0.02	0.00
proposal				(0.92)	(0.75)	(0.55)	(0.33)	(0.30)	(0.24)	(0.16)	(0.15)	(0.12)
			0.2	0.37	0.24	0.01	0.21	0.12	0.01	0.05	0.03	0.00
				(1.00)	(0.85)	(0.67)	(0.56)	(0.46)	(0.33)	(0.22)	(0.20)	(0.16)
		0.4	0.1	0.11	0.23	0.01	0.03	0.07	0.00	0.01	0.02	0.00
				(0.47)	(0.60)	(0.43)	(0.20)	(0.26)	(0.20)	(0.10)	(0.12)	(0.10)
			0.2	0.05	0.27	0.01	0.05	0.13	0.01	0.01	0.04	0.00
				(0.62)	(0.84)	(0.59)	(0.32)	(0.43)	(0.30)	(0.15)	(0.20)	(0.15)
	2	0.2	0.1	0.45	0.25	0.01	0.15	0.09	0.01	0.03	0.02	0.00
				(0.89)	(0.74)	(0.55)	(0.36)	(0.32)	(0.24)	(0.16)	(0.15)	(0.12)
			0.2	0.24	0.14	0.02	0.23	0.14	0.01	0.08	0.05	0.00
				(0.97)	(0.81)	(0.67)	(0.65)	(0.51)	(0.37)	(0.30)	(0.25)	(0.18)
		0.4	0.1	0.12	0.25	0.02	0.03	0.07	0.01	0.01	0.02	0.00
				(0.44)	(0.61)	(0.44)	(0.21)	(0.27)	(0.20)	(0.10)	(0.13)	(0.10)
			0.2	0.05	0.24	0.02	0.04	0.11	0.01	0.02	0.04	0.00
				(0.61)	(0.85)	(0.62)	(0.38)	(0.47)	(0.33)	(0.20)	(0.24)	(0.17)
${\theta }_{\left(1\text{-}2\right)}$ pools the results for ${\theta }_{\left(1\right)}$ and ${\theta }_{\left(2\right)}$ which are numerically very close, and, similarly, ${\theta }_{\left(3\text{-}6\right)}$ pools the results for ${\theta }_{\left(3\right)}$, ${\theta }_{\left(4\right)}$, ${\theta }_{\left(5\right)}$, and ${\theta }_{\left(6\right)}$. Absolute bias is shown here (which differs from the other tables) so that the pooling of results can be carried out.

Table 7. Table 7: Average bias and standard error (in brackets) of estimates: Other covariate types

				$n=100$			$n=300$			$n=1000$
Covariates	$\nu$	${\pi }_m$	$\rho$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
Binary	0	0.2	0.1	-0.20	0.10	0.10	-0.07	0.04	0.03	-0.02	0.01	0.01
				(0.61)	(0.52)	(0.50)	(0.30)	(0.26)	(0.24)	(0.15)	(0.14)	(0.13)
			0.2	-0.13	0.08	0.10	-0.08	0.06	0.04	-0.02	0.01	0.01
				(0.68)	(0.59)	(0.58)	(0.38)	(0.34)	(0.29)	(0.19)	(0.17)	(0.15)
		0.4	0.1	-0.04	0.08	0.07	-0.02	0.02	0.02	0.00	0.01	0.01
				(0.34)	(0.40)	(0.38)	(0.18)	(0.21)	(0.19)	(0.10)	(0.11)	(0.11)
			0.2	-0.03	0.09	0.08	-0.01	0.04	0.04	0.00	0.02	0.01
				(0.45)	(0.55)	(0.52)	(0.25)	(0.31)	(0.29)	(0.13)	(0.15)	(0.15)
Continuous	0	0.2	0.1	-0.19	0.10	0.12	-0.06	0.04	0.04	-0.02	0.01	0.01
& Binary				(0.64)	(0.56)	(0.55)	(0.30)	(0.29)	(0.26)	(0.14)	(0.14)	(0.13)
			0.2	-0.19	0.13	0.12	-0.07	0.05	0.04	-0.02	0.01	0.01
				(0.80)	(0.75)	(0.64)	(0.38)	(0.37)	(0.31)	(0.18)	(0.18)	(0.15)
		0.4	0.1	-0.04	0.11	0.07	-0.01	0.03	0.02	0.00	0.01	0.00
				(0.34)	(0.46)	(0.37)	(0.18)	(0.24)	(0.19)	(0.10)	(0.12)	(0.10)
			0.2	-0.04	0.13	0.12	-0.02	0.06	0.05	0.00	0.02	0.01
				(0.51)	(0.70)	(0.58)	(0.26)	(0.37)	(0.29)	(0.13)	(0.18)	(0.14)

Table 8. Table 8: Empirical coverage of 95% bootstrapped confidence intervals

			$n=100$			$n=300$			$n=1000$
$\nu$	${\pi }_m$	$\rho$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
0	0.2	0.1	91.2	93.2	93.9	93.5	93.2	94.3	94.0	94.8	94.0
		0.2	94.7	95.0	94.8	94.8	94.5	93.2	93.8	93.0	93.5
	0.4	0.1	94.6	93.6	93.4	94.8	94.0	94.9	94.7	94.4	94.8
		0.2	94.7	94.2	93.8	93.9	94.4	93.8	94.7	94.0	94.2
2	0.2	0.1	92.9	94.3	93.4	93.9	93.5	94.8	94.2	94.3	94.9
		0.2	95.9	95.3	96.1	94.5	94.9	94.1	93.0	95.0	95.3
	0.4	0.1	93.3	93.0	92.4	93.9	94.6	94.6	95.0	93.8	95.0
		0.2	94.7	95.3	94.1	94.0	94.4	95.3	94.2	94.5	94.9

Table 9. Table 9: Correct zeros, incorrect zeros, and model degrees of freedom

				$n=100$			$n=300$			$n=1000$
Type	$\nu$	${\pi }_m$	$\rho$	C	IC	DF	C	IC	DF	C	IC	DF
oracle				4.00	0.00	3.00	4.00	0.00	3.00	4.00	0.00	3.00
lasso	0	0.2	0.1	2.86	0.37	3.77	2.39	0.00	4.61	2.31	0.00	4.69
			0.2	3.23	0.79	2.98	2.66	0.06	4.29	2.54	0.00	4.46
		0.4	0.1	2.52	0.11	4.37	2.33	0.00	4.67	2.30	0.00	4.70
			0.2	2.90	0.45	3.65	2.51	0.02	4.47	2.44	0.00	4.56
	2	0.2	0.1	2.76	0.34	3.90	2.40	0.00	4.60	2.28	0.00	4.72
			0.2	3.19	0.85	2.96	2.64	0.09	4.27	2.46	0.00	4.53
		0.4	0.1	2.58	0.10	4.32	2.37	0.00	4.63	2.24	0.00	4.76
			0.2	2.92	0.48	3.60	2.50	0.05	4.45	2.43	0.00	4.57
alasso	0	0.2	0.1	3.44	0.37	3.19	3.47	0.01	3.52	3.67	0.00	3.33
			0.2	3.50	0.74	2.76	3.50	0.09	3.41	3.69	0.00	3.31
		0.4	0.1	3.34	0.17	3.49	3.52	0.00	3.48	3.69	0.00	3.31
			0.2	3.40	0.47	3.13	3.45	0.05	3.50	3.67	0.00	3.32
	2	0.2	0.1	3.35	0.36	3.30	3.48	0.01	3.51	3.65	0.00	3.35
			0.2	3.47	0.77	2.76	3.46	0.13	3.42	3.65	0.00	3.35
		0.4	0.1	3.33	0.15	3.52	3.55	0.00	3.44	3.72	0.00	3.28
			0.2	3.38	0.47	3.15	3.42	0.07	3.50	3.64	0.01	3.35

Table 10. Table 10: Standard error of estimates using different B 𝐵 B types

				$n=100$			$n=300$			$n=1000$
$\nu$	${\pi }_m$	$\rho$	$B$ Type	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$	${\theta }_{\left(0\right)}$	${\theta }_{\left(1\right)}$	${\theta }_{\left(2\right)}$
0	0.2	0.1	$B^{\ast }\left(S_C\right)$	(0.89)	(0.72)	(0.72)	(0.46)	(0.39)	(0.36)	(0.21)	(0.19)	(0.17)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.65)	(0.62)	(0.58)	(0.29)	(0.29)	(0.26)	(0.15)	(0.14)	(0.14)
			$B$	(0.39)	(0.36)	(0.37)	(0.21)	(0.19)	(0.19)	(0.11)	(0.11)	(0.10)
		0.2	$B^{\ast }\left(S_C\right)$	(1.06)	(0.88)	(0.92)	(0.74)	(0.60)	(0.58)	(0.38)	(0.31)	(0.28)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.80)	(0.74)	(0.68)	(0.40)	(0.40)	(0.34)	(0.20)	(0.20)	(0.17)
			$B$	(0.38)	(0.37)	(0.36)	(0.21)	(0.20)	(0.19)	(0.11)	(0.10)	(0.10)
	0.4	0.1	$B^{\ast }\left(S_C\right)$	(0.40)	(0.53)	(0.50)	(0.22)	(0.28)	(0.26)	(0.12)	(0.14)	(0.13)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.36)	(0.49)	(0.47)	(0.18)	(0.24)	(0.22)	(0.10)	(0.12)	(0.12)
			$B$	(0.25)	(0.32)	(0.30)	(0.14)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)
		0.2	$B^{\ast }\left(S_C\right)$	(0.68)	(0.84)	(0.85)	(0.39)	(0.50)	(0.44)	(0.19)	(0.24)	(0.21)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.51)	(0.74)	(0.66)	(0.29)	(0.37)	(0.32)	(0.13)	(0.18)	(0.16)
			$B$	(0.26)	(0.31)	(0.30)	(0.14)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)
2	0.2	0.1	$B^{\ast }\left(S_C\right)$	(0.93)	(0.73)	(0.72)	(0.49)	(0.40)	(0.37)	(0.22)	(0.19)	(0.17)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.62)	(0.60)	(0.51)	(0.29)	(0.29)	(0.24)	(0.15)	(0.15)	(0.13)
			$B$	(0.37)	(0.36)	(0.35)	(0.21)	(0.19)	(0.20)	(0.11)	(0.11)	(0.11)
		0.2	$B^{\ast }\left(S_C\right)$	(1.12)	(0.86)	(0.87)	(0.73)	(0.59)	(0.53)	(0.41)	(0.33)	(0.28)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.80)	(0.76)	(0.60)	(0.45)	(0.43)	(0.32)	(0.21)	(0.20)	(0.15)
			$B$	(0.38)	(0.35)	(0.36)	(0.21)	(0.20)	(0.20)	(0.11)	(0.10)	(0.11)
	0.4	0.1	$B^{\ast }\left(S_C\right)$	(0.45)	(0.58)	(0.50)	(0.23)	(0.29)	(0.25)	(0.12)	(0.14)	(0.12)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.37)	(0.53)	(0.44)	(0.19)	(0.25)	(0.21)	(0.10)	(0.13)	(0.11)
			$B$	(0.26)	(0.30)	(0.30)	(0.14)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)
		0.2	$B^{\ast }\left(S_C\right)$	(0.73)	(0.91)	(0.83)	(0.38)	(0.47)	(0.36)	(0.20)	(0.25)	(0.18)
			$B^{\ast }\left({\hat{S}}_C\right)$	(0.55)	(0.77)	(0.60)	(0.30)	(0.41)	(0.27)	(0.15)	(0.20)	(0.13)
			$B$	(0.26)	(0.29)	(0.30)	(0.15)	(0.17)	(0.17)	(0.08)	(0.09)	(0.09)

Equations

$$T_0$$

$$B$$

$$\mathbb{E}\left(\frac{\Delta }{S_C(Y-\mid X)}|X\right)=\mathbb{E}(1-B\mid X)=1-\pi (X)$$

$$\mathbb{E}\left(\frac{\Delta r(Y,X)}{S_C(Y-\mid X)}|X\right)=\mathbb{E}\{r(T_0,X)\mid X\}\{1-\pi (X)\}$$

$$\ell (\theta )=\sum _{i=1}^n\left[B_i\log {\pi }_i(\theta )+(1-B_i)\log (1-{\pi }_i(\theta ))\right],$$

$$B_i^{\ast }(S_C)=B_i^{\ast }(Y_i,{\Delta }_i,X_i,S_C)=1-\frac{{\Delta }_i}{S_C(Y_i-\mid X_i)},$$

$${\ell }^{\ast }(\theta )=\sum _{i=1}^n\left[B_i^{\ast }(S_C)\log {\pi }_i(\theta )+\{1-B_i^{\ast }(S_C)\}\log (1-{\pi }_i(\theta ))\right]{\omega }_i.$$

$$U_i^{\ast }(\theta )=\frac{\left[B_i^{\ast }(S_C)-{\pi }_i(\theta )\right]{\omega }_i}{{\pi }_i(\theta )(1-{\pi }_i(\theta ))}\frac{\partial {\pi }_i(\theta )}{\partial \theta }\in {\mathbb{R}}^{p+1},$$

$$\hat{\theta }=({\hat{\theta }}_{(0)},{\hat{\theta }}_{(1)},\dots ,{\hat{\theta }}_{(p)}{)}^{\top }=\arg \underset{\theta }{\max }{\hat{\ell }}^{\ast }(\theta )$$

$${\hat{\ell }}^{\ast }(\theta )=\sum _{i=1}^n\left[B_i^{\ast }({\hat{S}}_C)\log {\pi }_i(\theta )+\{1-B_i^{\ast }({\hat{S}}_C)\}\log (1-{\pi }_i(\theta ))\right]{\omega }_i.$$

$${\hat{\theta }}_{\lambda }=({\hat{\theta }}_{\lambda ,(0)},{\hat{\theta }}_{\lambda ,(1)},\dots ,{\hat{\theta }}_{\lambda ,(p)}{)}^{\top }=\arg \underset{\theta \in \Theta }{\max }{\hat{\ell }}_{\lambda }^{\ast }(\theta )$$

$$\underset{\theta \in \Theta }{\sup }\left|\frac{1}{n}\sum _{i=1}^n\left[B_i^{\ast }({\hat{S}}_C)-B_i^{\ast }(S_C)\right]\omega (X_i)\log \left(\frac{{\pi }_i(\theta )}{1-{\pi }_i(\theta )}\right)\right|=o_{\mathbb{P}}(1).$$

$$\underset{x\in \mathcal{X},\omega (x)>0}{\sup }\underset{y\le \tau (x)}{\sup }\left|{\hat{S}}_C(y-\mid x)-S_C(y-\mid x)\right|=o_{\mathbb{P}}(1),$$

$$\mathbb{E}\left[{\Vert \frac{\partial \pi ({\theta }_0)}{\partial \theta }\Vert }^2\right]<\infty$$

$$A({\theta }_0)=\mathbb{E}\left[\frac{\omega (X)}{\pi ({\theta }_0)[1-\pi ({\theta }_0)]}\frac{\partial \pi ({\theta }_0)}{\partial \theta }{\frac{\partial \pi ({\theta }_0)}{\partial \theta }}^{\top }\right]$$

$${\mathcal{F}}_{1,kl}$$

$${\mathcal{F}}_{2,kl}$$

$$\frac{1}{n}\sum _{1\le i\le n}\left[B_i^{\ast }({\hat{S}}_C)-B_i^{\ast }(S_C)\right]\phi (X_i)=\frac{1}{n}\sum _{1\le i\le n}{\mu }_C^{\phi }(Y_i,{\Delta }_i,X_i)+o_{\mathbb{P}}(n^{-1/2});$$

$$\phi (X)=\frac{\omega (X)}{\pi ({\theta }_0)[1-\pi ({\theta }_0)]}\frac{\partial \pi ({\theta }_0)}{\partial \theta }.$$

$$\begin{gathered} \hat{\theta }-{\theta }_0=A({\theta }_0{)}^{-1}\frac{1}{n}\sum _{i=1}^n\left\{\mu (Y_i,{\Delta }_i,X_i;{\theta }_0)+{\mu }_C^{\phi }(Y_i,{\Delta }_i,X_i;{\theta }_0)\right\} \\ +o_{\mathbb{P}}(n^{-1/2}). \end{gathered}$$

$$\mu (Y,\Delta ,X;{\theta }_0)=\frac{\left[B^{\ast }(S_C)-\pi ({\theta }_0)\right]\omega (X)}{\pi ({\theta }_0)[1-\pi ({\theta }_0)]}\frac{\partial \pi ({\theta }_0)}{\partial \theta }$$

$$\sqrt{n}\left(\hat{\theta }-{\theta }_0\right)⇝N_{p+1}\left(0,A({\theta }_0{)}^{-1}V({\theta }_0)A({\theta }_0{)}^{-1}\right)$$

$${\hat{\theta }}_{\mathcal{A},\lambda }-{\theta }_{\mathcal{A},0}$$

$$=A_{\mathcal{A}}({\theta }_0{)}^{-1}\frac{1}{n}\sum _{i=1}^n\left\{{\mu }_{\mathcal{A}}(Y_i,{\Delta }_i,X_i;{\theta }_0)+{\mu }_{\mathcal{A},C}^{\phi }(Y_i,{\Delta }_i,X_i;{\theta }_0)\right\}+o_{\mathbb{P}}(n^{-1/2}).$$

$${\mu }_{\mathcal{A}}(Y,\Delta ,X;{\theta }_0)=\frac{\left[B^{\ast }(S_C)-\pi ({\theta }_0)\right]\omega (X)}{\pi ({\theta }_0)[1-\pi ({\theta }_0)]}\frac{\partial \pi ({\theta }_0)}{\partial {\theta }_{\mathcal{A}}}$$

$${\phi }_{\mathcal{A}}(X)=\frac{\omega (X)}{\pi ({\theta }_0)[1-\pi ({\theta }_0)]}\frac{\partial \pi ({\theta }_0)}{\partial {\theta }_{\mathcal{A}}}.$$

$$\sqrt{n}\left({\hat{\theta }}_{\mathcal{A},\lambda }-{\theta }_{\mathcal{A},0}\right)⇝N_{p_0+1}\left(0,A_{\mathcal{A}}({\theta }_0{)}^{-1}{V}_{\mathcal{A}}({\theta }_0)A_{\mathcal{A}}({\theta }_0{)}^{-1}\right)$$

$$S_{T_0}(t\mid X)={\left\{\frac{\exp (-t^{\kappa })-\exp (-{\tau }^{\kappa })}{1-\exp (-{\tau }^{\kappa })}\right\}}^{\psi }$$

$$\mathbb{E}\left(\frac{\Delta r(Y,X)}{S_C(Y-\mid X)}|X\right)=\mathbb{E}\{r(T_0,X)\mid X\}\{1-\pi (X)\}.$$

$$\mathbb{E}(\Delta \mid T_0,X)$$

$$=\mathbb{E}\{\text{mathbbm}1(T_0\le C)\mathbb{E}(1-B\mid C,T_0,X)\mid T_0,X\}$$

$$=S_C(T_0-\mid X)\{1-\pi (X)\},$$

$$\mathbb{E}\left(\frac{\Delta r(Y,X)}{S_C(Y-\mid X)}|X\right)$$

Full text

Keywords. Binary regression; iid representation; Inverse probability censoring weighting; Penalized likelihood.

MSC 2010 subject classification: 62N01, 62N02, 62J07

1 Introduction

Standard survival models assume that all individuals experience the event of interest eventually (see Kalbfleisch and Prentice, (2002)). However, this assumption is not always tenable since, for example, some diseases might only activate when specific biological and/or lifestyle traits are present, or immunity may result from successful curative treatment (or a combination of both treatment and pre-treatment attributes); individuals who will never experience the event are referred to as cured or non-susceptible. In practice, the binary “cure status”, $B\in\{0,1\}$ (where $B=1\Leftrightarrow$ “the individual is cured”), typically cannot be measured directly, and, so, survival studies are required.

Let $T\in(0,\infty]$ denote the survival time, and note in particular that, in contrast to standard survival models, the support includes $T=\infty$ corresponding to cured individuals. Thus, $B=\mathbbm{1}(T=\infty)\sim\text{Bernoulli}(\pi)$ where $\pi$ is the cure probability. Furthermore, let $T_{0}\in(0,\infty)$, be the “latency” time (i.e., the survival time for an uncured individual) with survivor function $S_{T_{0}}(t)=\Pr(T_{0}>t)$. By construction, $T=T_{0}$ when $B=0$, and, since $\Pr(T>t\mid B)=(1-B)S_{T_{0}}(t)+B$, we have that $S_{T}(t)=\mathbb{E}\{\Pr(T>t\mid B)\}=$ $(1-\pi)S_{T_{0}}(t)+\pi$ such that $\lim_{t\rightarrow\infty}S_{T}(t)=\pi$. Most often, interest centres on modelling the effect of covariates, $X$, on $\pi$. Thus, covariates enter via $S_{T}(t\mid X)=\{1-\pi(X)\}S_{T_{0}}(t\mid X)+\pi(X)$ where $\pi(X)=\pi(X;\theta)\in(0,1]$ is the cure regression function, and $\theta$ is a vector of parameters; typically, although not necessarily, covariates enter through a linear predictor, i.e., $\pi(X;\theta)=\pi(X^{\top}\theta)$ where $X=(X_{(0)}\equiv 1,X_{(1)},\ldots,X_{(p)})^{\top}$ and $\theta=(\theta_{(0)},\theta_{(1)},\ldots,\theta_{(p)})^{\top}$ respectively. (Here and in the following, for any matrix $A$, $A^{\top}$ denotes its transpose.)

Parametric cure models have been considered by various authors (Berkson and Gage, , 1952; Farewell, , 1977; Peng et al., , 1998), but the $\theta$ estimates can be sensitive to the choice of latency model (Yu et al., , 2004). The semi-parametric cure model consists of a parametric cure regression model and a semi-parametric proportional hazards (PH) or accelerated failure time (AFT) latency model (Peng and Dear, , 2000; Sy and Taylor, , 2000; Li and Taylor, , 2002; Zhang and Peng, , 2007); estimation requires the EM algorithm (Dempster et al., , 1977) in combination with (modified) partial likelihood (Cox, , 1972, 1975) (PH case) or rank regression (Ritov, , 1990; Tsiatis, , 1990) (AFT case). However, structural model assumptions are made in the way that covariates enter the latency component (i.e., PH or AFT structure), and estimation of the cure component can be sensitive to such choices. (See Amico and Van Keilegom, (2018) for an excellent review of cure literature.)

Our approach differs from the aforementioned. First note that, if $B$ was directly observable, one would simply apply standard binary regression (e.g., logistic or probit) without modelling $S_{T_{0}}$. We define a “proxy” or “synthetic” variable which we denote by $B^{*}$, then replace $B$ by $B^{*}$ and proceed with classical binary regression. We construct $B^{*}$ using inverse probability censoring weighting (IPCW) arguments (Robins and Finkelstein, , 2000; Van der Laan et al., , 2003). This approach obviates the need for a latency model – although IPCW does require estimation of the censoring distribution. At first sight, we trade one missing data framework (EM) for another (IPCW), and latency estimation for censoring estimation. However, our novel procedure has key advantages: (i) unlike existing EM approaches, which are not so easily generalized beyond PH and AFT latency models, our framework straightforwardly permits a wide range of censoring models, (ii) the censoring model can be validated using standard model checking techniques since, unlike the latency model, it is directly identifiable from the observed data, and (iii) once $B^{*}$ has been computed (in a single initial step rather than in an iterative EM fashion), one may avail of standard, fast (penalized) GLM estimation procedures.

The remainder of this article is organized as follows. In Section 2, we introduce a new result which forms the basis of the proposed estimation procedure of Section 3, including a penalized version for the purpose of variable selection. Asymptotic theory is provided in Section 4, along with empirical evidence via simulation in Section 5. A real data example is given in Section 6, and we close with some remarks in Section 7. The proofs are postponed to the Appendix (Section 8).

2 Preliminaries

Before we proceed, we define the censoring time $C\in(0,\infty)$ with survivor function $S_{C}(t)$, whose support excludes infinity since, practically, observation windows are finite. Furthermore, let $Y=T\wedge C$ be the observed time ($\wedge$ is the minimum operator) and $\Delta=\mathbbm{1}(T\leq C)=(1-B)\mathbbm{1}(T_{0}\leq C)$ be the event indicator. We will assume the following:

[TABLE]

where (2.1) is the standard independence assumption made throughout survival literature, and (2.2) is introduced in the cure context to ensure that the cure regression model is identifiable. In particular, assumptions (2.1) and (2.2) guarantee that $T\perp C\mid X$ (see the Appendix, Lemma 8.2). With these assumptions in place, it can then be shown that

[TABLE]

where $S_{C}(t-\mid X)=\Pr(C\geq t\mid X).$ This result is the core of our estimation scheme which is described in Section 3. In fact, (2.3) is a special case of the more general result

[TABLE]

but with $r(\cdot,\cdot)=1$; indeed, it is (2.4) which is proved in Lemma 8.1. It is worth highlighting that (2.4) is an application of the Inverse-Probability Censoring Weighting (IPCW) approach (Robins and Finkelstein, , 2000) extended to the case where a cured proportion exists.

3 Estimation and inference

If independent and identically distributed (iid) replicates of $(B_{i},X_{i})$ were observed, we would use the Bernoulli log-likelihood,

[TABLE]

where $\pi_{i}(\theta)=\pi(X_{i};\theta)$. Of course, (3.5) is not operational in the current context since $B_{i}$ is unobserved, but it serves as our motivation for estimation and inference using iid replicates $(Y_{i},\Delta_{i},X_{i})$.

3.1 Likelihood-based estimation

We now define

[TABLE]

such that $\mathbb{E}\{B_{i}^{*}(S_{C})\mid X_{i}\}=\mathbb{E}(B_{i}\mid X_{i})$ follows from (2.3). Assuming initially that $S_{C}$ is known, we propose replacing $B_{i}$ with $B_{i}^{*}(S_{C})$ in (3.5) to obtain

[TABLE]

Unlike $\ell(\theta)$ however, $\ell^{*}(\theta)$ is formed using the observable quantities $Y_{i}$ and $\Delta_{i}$ rather than the unobservable $B_{i}$. Here, $\omega_{i}=\omega(X_{i})$ are positive weights which are introduced as a technical device when deriving general asymptotic results (but we anticipate that $\omega_{i}\equiv 1$ in practically all applications). In the sense that $\ell^{*}(\theta)$ is based on $\ell(\theta)$, we refer to our procedure as “likelihood-based”. While the phrase “pseudo-likelihood” has been used in similar contexts to ours (Xie and Liu, , 2005), we prefer to avoid this terminology due to its long-standing usage in settings where incorrect models can still yield valid inferences (Gourieroux et al., , 1984); our setting is rather different to this in that the model remains unchanged, but the unobservable $B_{i}$ is replaced with the observable surrogate $B_{i}^{*}$. The score function, $U^{*}(\theta)=\partial\ell^{*}(\theta)/\partial\theta=\sum_{i=1}^{n}U^{*}_{i}(\theta)$ where

[TABLE]

is unbiased due to property (2.3), and, in the case of a logistic cure regression model, $U_{i}^{*}(\theta)=[B_{i}^{*}(S_{C})-\pi_{i}(\theta)]X_{i}\omega_{i}$ with $\pi_{i}(\theta)=1/\{1+\exp(-\theta^{\top}X_{i})\}$. Furthermore, using standard inequalities, we show in the Appendix (Lemma 8.3) that, when the cure regression model is identifiable, $\mathbb{E}\{\ell^{*}(\theta)\}<\mathbb{E}\{\ell^{*}(\theta_{0})\}$ $\forall\,\theta\neq\theta_{0},$ where $\theta_{0}$ is the true parameter vector.

Note that we have explicitly written $B^{*}_{i}(S_{C})$ as a function of $S_{C}$ (while its dependence on $Y_{i}$, $\Delta_{i}$, and $X_{i}$ is implicit) since, typically, $S_{C}$ must be estimated in practice, i.e., we will use $B_{i}^{*}(\widehat{S}_{C})$; this is standard in IPCW applications (and recall that (2.3) is based on IPCW arguments). Furthermore, the asymptotic theory of Section 4 requires only that the estimator for $S_{C}$ has an iid representation. Thus, while $S_{T_{0}}$ is completely unspecified in our proposal, $S_{C}$ is estimated, and this can essentially be done in an arbitrarily flexible way. Therefore, for practical purposes, we propose the cure estimator defined as

[TABLE]

where

[TABLE]

The weights $\omega_{i}=\omega(X_{i})$ serve in theory to control the behavior of general estimates of $B_{i}^{*}(S_{C})$ in regions of low covariate density. Since (3.10) is the usual Bernoulli log-likelihood with $B_{i}$ replaced by $\widehat{B}_{i}^{*}=B_{i}^{*}(\widehat{S}_{C})$, standard GLM estimation procedures can be used, and the EM algorithm is avoided. Interestingly, as discussed further in Section 4.2, it appears in our case that replacement of $B_{i}$ with $\widehat{B}^{*}_{i}$ produces more efficient cure parameter estimates than with $B^{*}_{i}$. (Of course, $B_{i}$ itself would produce the most efficient estimates if it were available.) Note that the $\widehat{B}^{*}_{i}$ variables play a somewhat similar role to “synthetic observations” as used by Koul et al., (1981) in a different context (see also Delecroix et al., (2008)).

3.2 Variable selection

With the above in place, by construction, it is also straightforward to include a penalty for the purpose of variable selection. In particular, we consider the alasso (adaptive least absolute shrinkage and selection operator) penalty of Zou, (2006). Although, so far, we have not specified the functional form of $\pi(X;\theta)$, here, we will assume that $\pi(X;\theta)=\pi(X^{\top}\theta)$, and, thus, $X\in\mathbb{R}^{p+1}$ with $X_{(0)}\equiv 1$. Then, the alasso estimator is

[TABLE]

where $\widehat{\ell}^{*}_{\lambda}(\theta)=\widehat{\ell}^{*}(\theta)-\lambda\sum_{j=1}^{p}w_{j}|\theta_{(j)}|$ with tuning parameter $\lambda\geq 0$ and (potentially adaptive) weights $w_{j}\geq 0$ for $1\leq j\leq p$. Here, as is usual, the intercept, $\theta_{(0)}$, is not penalized in (3.11), and, furthermore, typically, the covariates are standardized. Setting $w_{j}=1$ $\forall j$ yields the lasso penalty, which penalizes all coefficients equally (Tibshirani, , 1996), while $w_{j}=1/|\widehat{\theta}_{(j)}^{(0)}|^{\gamma}$ for some $\gamma>0$ yields the alasso penalty (we will set $\gamma=1$ as is most common in practice). In the latter case, $\widehat{\theta}_{(j)}^{(0)}$ may be any consistent estimator of $\theta_{(j)}$, and, typically, $\widehat{\theta}_{(j)}^{(0)}=\widehat{\theta}_{(j)}$, where $\widehat{\theta}_{(j)}$ is the $j$th unpenalized estimator from (3.9). Details on implementation aspects of the alasso in our context (optimization procedure and tuning parameter selection) can be found in the Appendix (Section 8.3).

4 Asymptotic results

Our asymptotic results are proved under some minimal moment assumptions on the observed variables completed by some mild high-level assumptions on the cure regression model and on the model for the censoring variable. These conditions are quite natural in the context of right-censored data when covariates are present and are to be verified on a case by case basis according to the context of the application. In this section we use the notation $\pi(\theta)=\pi(X;\theta),$ $\partial\pi(\theta)/\partial\theta=\partial\pi(X;\theta)/\partial\theta$, and $B^{*}_{i}(S_{C})$ (hence, $B^{*}_{i}(\widehat{S}_{C})$) as defined in equation (3.6).

Assumption 4.1

(The data)* The observations $(Y_{i},\Delta_{i},X_{i}),$ $1\leq i\leq n,$ are independent replicates of $(Y,\Delta,X)\in\mathbb{R}\times\{0,1\}\times\mathcal{X}$, where $\mathcal{X}$ is some covariate space. Moreover, $\mathbb{E}[\Delta/S_{C}(Y-\mid X)]<\infty.$*

Let $\mathcal{M}=\{\pi(\theta):\theta\in\Theta\subset\mathbb{R}^{p+1}\}$ be a generic parametric cure regression model, that is a set of functions of the covariate vector $X$ indexed by $\theta$ in some parameter space $\Theta$, e.g., the logistic cure model is $\pi(\theta)=1/\{1+\exp(-\theta^{\top}X)\}.$

4.1 Consistency

For the definition of Glivenko-Cantelli function classes, we refer to van der Vaart, (2000).

Assumption 4.2

(The cure regression model)**

The weight $\omega(X)$ is bounded, almost surely nonnegative and has a positive expectation. 2. 2.

There exists $\theta_{0}\in\Theta$ such that $\Pr(T=\infty\mid X)=\pi(\theta_{0})$. Moreover, there exists $0<c<1/2$ such that, for any $\theta\in\Theta$, $\Pr(c\leq\pi(\theta)\leq 1-c)=1.$ 3. 3.

For any $\varepsilon>0,$ $\inf_{\|\theta-\theta_{0}\|>\varepsilon}\mathbb{E}[|\pi(\theta)-\pi(\theta_{0})|\omega(X)]>0$. 4. 4.

*The model $\mathcal{M}$ is a *$\mathbb{P}_{X}-$Glivenko-Cantelli class of functions of $X$ with constant envelope.

Assumption 4.3

(Uniform law of large numbers)* The estimator $\widehat{S}_{C}(\cdot\mid\cdot)$ satisfies the law of large numbers uniformly over the class of the logit transformations of the functions in $\mathcal{M}$:*

[TABLE]

For simplicity, we consider a bounded weight $\omega(X)$ and assume that $\pi_{i}(\theta)$ stays uniformly away from 0 and 1 $\forall i$. Note that Assumption 4.2.3 provides an indentifiability condition which guarantees that $\theta_{0}$ is a well-separated maximum of $\mathbb{E}\{\ell^{*}(\theta)\}$. It is satisfied, for example, by logistic or probit regression models when the covariates are not redundant, i.e., when $\mathbb{E}[XX^{\top}\omega(X)]$ is an invertible matrix.

Let us provide some mild sufficient conditions implying the uniform convergence of Assumption 4.3. These sufficient conditions involve a threshold , commonly used in cure literature, which is typically justified as representing a total follow-up of the study. Usually, this is assumed to be independent of the covariates, but we allow it to depend on the covariates in an arbitrary way.

Lemma 4.1

Assume that there exists $\tau(x)$ such that, for any $x$, $\Pr(T_{0}>\tau(x))=0$ and $\inf_{x\in\mathcal{X},\;\omega(x)>0}S_{C}(\tau(x)-\mid x)>0.$ Moreover, Assumptions 4.2.2 and 4.2.4 hold true. If

[TABLE]

then the uniform convergence in Assumption 4.3 holds true.

The common parametric, semiparametric and nonparametric estimators $\widehat{S}_{C}$ satisfy condition (4.12). Several examples are provided in the monographs by Anderson et al., (1993) and Kalbfleisch and Prentice, (2002), and, for convenience, some examples are recalled in the Appendix (Section 8.2). The consistency of our cure estimator is stated in the following result.

Theorem 4.1

Let Assumptions 4.1, 4.2 and 4.3 hold true. Then $\widehat{\theta}-\theta_{0}=o_{\mathbb{P}}(1)$.

4.2 Asymptotic normality

Assumption 4.4

(The cure regression model)**

For any $x\in\mathcal{X}$, the map $\theta\mapsto\pi(x;\theta)$ is twice continuously differentiable. 2. 2.

The true value $\theta_{0}$ is an interior point of $\Theta$,

[TABLE]

*and the *$(p+1)\times(p+1)-$matrix

[TABLE]

is positive definite. 3. 3.

*For any $0\leq k\leq l\leq p,$ the families of functions of $x$ indexed by $\theta$ *

[TABLE]

*are *$\mathbb{P}_{X}-$Glivenko-Cantelli classes of functions of $X$ with integrable envelopes.

Assumption 4.5

(I.I.D. representation)* Let $\varphi(X)$ be a vector-valued function such that $\mathbb{E}\{\|\varphi(X)\|^{2}\}<\infty$. Then there exists $\mu_{C}^{\varphi}(Y,\Delta,X)$ a zero-mean vector-valued function that depends on $\varphi(X)$, such that $\mathbb{E}\{\|\mu_{C}^{\varphi}(Y,\Delta,X)\|^{2}\}<\infty$ and*

[TABLE]

Assumption 4.4 introduces mild standard regularity conditions on the cure regression model. In particular, Assumption 4.4.3 yields the uniform law of large numbers and guarantees that the remainder terms in the standard Taylor expansion used to establish asymptotic normality are uniformly negligible. Such an assumption on the complexity of the classes of first and second order derivatives of the functions in the model are satisfied by the standard parametric models such as logit and probit models. As an alternative to Assumption 4.4.3, we could impose condition (4.12) and slightly stronger regularity conditions on the model $\mathcal{M}$. The details are provided in the proof of Theorem 4.2. Furthermore, note that an asymptotic representation as required in Assumption 4.5 is very common in survival analysis models with ‘nuisance’ parameters which may belong to a space of functions. In the Appendix (Section 8.2), we provide details for many standard survival estimators: Kaplan-Meier, conditional Kaplan-Meier (Beran), Cox model, transformation model, and proportional odds model.

In general, the expression of the function $\mu_{C}^{\varphi}(Y,\Delta,X)$ in Assumption 4.5 depends on the joint law of $(Y,\Delta,X)$. Furthermore, this function contributes to the asymptotic variance of our estimator $\widehat{\theta}$ proposed in (3.9) which, hence, will differ from that of the infeasible estimator defined with $B^{*}(S_{C})$ instead of $B^{*}(\widehat{S}_{C})$.

Theorem 4.2

Assume the conditions of Assumptions 4.1 and 4.2, and Lemma 4.1 are met. Moreover, let Assumptions 4.4 and 4.5 hold true, and

[TABLE]

Then

[TABLE]

where

[TABLE]

and $\mu_{C}^{\varphi}(Y,\delta,X)$ is the zero-mean vector-valued function from Assumption 4.5. In addition,

[TABLE]

with $V(\theta_{0})=\mbox{\rm Var}\left\{\mu(Y,\Delta,X;\theta_{0})+\mu_{C}^{\varphi}(Y,\Delta,X;\theta_{0})\right\}$. ($\rightsquigarrow$ denotes convergence in law.)

Following an anonymous reviewer’s suggestion, let us analyse the variance $V(\theta_{0})$. It is straightforward to show that $Var\left\{\mu(Y,\Delta,X;\theta_{0})\right\}$, the variance obtained by maximizing (3.7) which uses $B^{*}_{i}(S_{C})$, is larger than that obtained by maximizing by (3.5) which uses $B_{i}$. (Of course, neither of these are feasible since $B_{i}$ and $S_{C}$ are unknown.) However, what seems to be much less well known is that, in many cases, the contribution of $\mu_{C}^{\varphi}(Y,\Delta,X;\theta_{0})$, due to estimating $S_{C}$, has the effect of making $V(\theta_{0})$ smaller than $Var\left\{\mu(Y,\Delta,X;\theta_{0})\right\}$. This perhaps unexpected effect, explained by Hitomi et al., (2008), often occurs in semiparametric estimation problems as considered herein. More specifically, their Theorem 3 shows that that $Var\left\{\mu(Y,\Delta,X;\theta_{0})\right\}-V(\theta_{0})$ is positive semi-definite provided that $\mu_{C}^{\varphi}(Y,\Delta,X;\theta_{0})$ belongs to the tangent space $\mathcal{T}$ with respect to the ‘nuisance’ parameter $S_{C}$. The space $\mathcal{T}$ is defined as the mean-square closure of all linear combinations of scores for smooth parametric submodels passing through $S_{C}$. When there is no restriction on the form of $S_{C}$, this tangent space is typically large and very likely to include $\mu_{C}^{\varphi}(Y,\Delta,X;\theta_{0})$; from Remark 2 of Hitomi et al., (2008), this is also true for parametric $S_{C}$ models estimated using maximum likelihood. Finally, Hitomi et al., (2008) show (again in their Theorem 3) that $V(\theta_{0})$ is strictly less than $Var\left\{\mu(Y,\Delta,X;\theta_{0})\right\}$ if the projection of the score function (3.8) onto $\mathcal{T}$ is non-null. We claim that this is indeed the case for many possible $S_{C}$ models, and a formal investigation of this will be a focus of our future work. However, we provide numerical experiments which support this claim in Section 5.

To estimate the variance of $\widehat{\theta}$, we require estimates of $V(\theta_{0})$ and $A(\theta_{0})$, respectively. If estimates of the vectors $\mu(Y_{i},\Delta_{i},X_{i};\theta)$ and $\mu_{C}^{\varphi}(Y_{i},\Delta_{i},X_{i};\theta)$ are available, say, $\widehat{\mu}_{i}(\theta)$ and $\widehat{\mu}_{C,i}^{\varphi}(\theta)$, then $V(\theta_{0})$ can be estimated by the sample covariance, $n^{-1}\sum_{i=1}^{n}[\widehat{\mu}_{i}(\widehat{\theta})+\widehat{\mu}_{C,i}^{\varphi}(\widehat{\theta})]^{\bigotimes 2}$ where $a^{\bigotimes 2}=aa^{T}$. Meanwhile, $A(\theta_{0})$ can also be estimated by standard methods. However, the estimates $\widehat{\mu}_{i}(\theta)$ and $\widehat{\mu}_{C,i}^{\varphi}(\theta)$ are often quite intractable. Thus, one may alternatively use the nonparametric bootstrap; indeed, this approach works well empirically (see Appendix, Section 8.4) and is used in our real data analysis.

4.3 Oracle properties for the adaptive lasso

Next, we prove consistency in variable selection for the adaptive lasso proposed in Section 3.2. Moreover, we prove the asymptotic normality for the true subset of coefficients. Hence, we extend the Theorem 4 of (Zou, , 2006) to the cure regression context. Let $\theta_{0}=(\theta_{0,(0)},\theta_{0,(1)},\ldots,\theta_{0,(p)})^{\top}$ be the true cure regression parameter vector. Assume the true model has a sparse representation, where $\mathcal{A}=\left\{j:1\leq j\leq p,\theta_{0,(j)}\neq 0\right\}\cup\{0\}$. Without loss of generality, suppose $\mathcal{A}=\{0,1,\ldots,p_{0}\}$, $p_{0}<p.$ Below, the subscript $\mathcal{A}$ is used to define the subvectors or blocks in matrices with components corresponding to the indices in the set $\mathcal{A}$, i.e., $\theta_{\mathcal{A},0}$ is the subvector of the first $p_{0}+1$ components of $\theta_{0}$, $\partial\pi(\theta_{0})/\partial\theta_{\mathcal{A}}$ denotes the vector of partial derivatives with respect to the first $p_{0}+1$ components of $\theta$, and $A_{\mathcal{A}}(\theta_{0})$ is the upper-left block of dimension $(p_{0}+1)\times(p_{0}+1)$ of the matrix $A(\theta_{0})$ defined in Assumption 4.4.2.

Theorem 4.3

Assume the conditions of Theorem 4.2 are met and $\pi(X;\theta)$ is a given function of $X^{\top}\theta$. Let $\widehat{\theta}_{\lambda}$ be the estimator defined in (3.11) with $w_{j}=|\widehat{\theta}_{(j)}|^{-\gamma}$ for $\gamma>0$. Moreover, assume that $\lambda/\sqrt{n}\rightarrow 0$ and $\lambda n^{(\gamma-1)/2}\rightarrow\infty$. Let $\mathcal{A}_{n}=\left\{j:1\leq j\leq p,\widehat{\theta}_{\lambda,(j)}\neq 0\right\}\cup\{0\}.$ Then

$\lim_{n\rightarrow\infty}\mathbb{P}(\mathcal{A}_{n}=\mathcal{A})=1$.

** 2. 2.

[TABLE]

where

[TABLE]

and $\mu_{\mathcal{A},C}^{\varphi}(Y,\Delta,X;\theta_{0})$ is the zero-mean vector-valued function from Assumption 4.5 considered with

[TABLE]

In addition,

[TABLE]

with $V_{\mathcal{A}}(\theta_{0})=Var\left\{\mu_{\mathcal{A}}(Y,\Delta,X;\theta_{0})+\mu_{\mathcal{A},C}^{\varphi}(Y,\Delta,X;\theta_{0})\right\}$.

As was the case for Theorem 4.2, we can obtain Theorem 4.3 by imposing condition (4.12) and slightly stronger regularity conditions on the model $\mathcal{M}$ instead of Assumption 4.4.3.

5 Simulation studies

5.1 Setup

We first generate $B$, $T_{0}$ and $C$, from which we obtain $T=T_{0}$ when $B=0$ and $T=\infty$ otherwise, and, hence, the observed time, $Y=T\wedge C$, and censoring indicator, $\Delta=(1-B)\mathbbm{1}(T_{0}\leq C)$, respectively. The cure status is given by $B\sim\text{Bernoulli}(\pi)$ where $\pi(\theta)=1/\{1+\exp(-X^{\top}\theta)\}$, $X=(1,X_{(1)},X_{(2)})^{\top}$, and $X_{(1)}$ and $X_{(2)}$ are independent $\text{Normal}(0,1)$ variables. We set $\theta_{0}=(\theta_{0,(0)},1,1)^{\top}$ with $\theta_{0,(0)}\in\{-1.85,-0.55\}$ such that the marginal cure proportion $\pi_{m}=\mathbb{E}\{\pi(\theta)\}\in\{0.2,0.4\}$. Consider the survivor function

[TABLE]

which is that of a truncated Weibull whose support is $(0,\tau)$ with a rate parameter, $\psi$, and a shape parameter, $\kappa$. The latency time, $T_{0}$, was generated according to this distribution with $\psi=\exp(X^{\top}\beta_{T_{0}})$ and $\kappa=(1/\psi)^{\nu}$ where $\beta_{T_{0}}=(0,0,1)^{\top}$ and $\nu\in\{0,2\}$; the proportional hazards property holds when $\nu=0$. The value of $\tau$ was set at the 95th percentile of the marginal untruncated distribution, i.e., $\tau$ is the unique solution $\mathbb{E}\left\{\exp(-\psi\tau^{\kappa})\right\}=0.05$, and, clearly, $\tau$ depends on the value of $\nu$. Lastly, the censoring time, $C$, was generated from an exponential distribution with rate parameter $\psi_{C}=\exp(X^{\top}\beta_{C})$ where $\beta_{C}=(\beta_{C,(0)},0,1)^{\top}$. The value $\beta_{C,(0)}$ was chosen such that the overall censored proportion is given by $\pi_{\text{cen}}=\Pr(\Delta=0)=\pi_{m}+\rho$ where $\rho\in\{0.1,0.2\}$, and this depends on the values of $\theta_{(0)}$ and $\nu$; since $\pi_{m}\in\{0.2,0.4\}$, there are then four values for the censoring proportion, $\pi_{\text{cen}}\in\{0.3,0.4,0.5,0.6\}$.

It is worth highlighting that $X_{(1)}$ only affects cure probability (since $\beta_{T_{0},(1)}=\beta_{C,(1)}=0$), whereas $X_{(2)}$ affects all components of the data generating process (since $\theta_{(2)}=\beta_{T_{0},(2)}=\beta_{C,(2)}=1$). Sample sizes of $n\in\{100,300,1000\}$ were considered, and, with two values for each of $\theta_{(0)}$, $\nu$, and $\rho$, there are 24 scenarios altogether. Each simulation scenario was replicated 2000 times. Here we report only on the 6 scenarios where $\pi_{m}=0.4$ and $\rho=0.1$. The results for the remaining scenarios are broadly similar, as are other scenarios with binary covariates (see Appendix, Section 8.4).

5.2 Estimation procedure

We applied the estimation scheme described in Section 3 to the simulated data with $S_{C}$ estimated using a Cox model in which both covariates, $X_{(1)}$ and $X_{(2)}$, appear as predictors. Table 1 displays the average bias and standard error of estimates over simulation replicates. While the bias can be somewhat large when $n=100$, this vanishes as the sample size increases. Similarly, the standard errors also decrease with the sample size. Furthermore, the results do not change appreciably when $\nu$ is varied (i.e., the approach is not sensitive to the form of $S_{T_{0}}$).

By way of comparison, we also applied the EM approach of Peng and Dear, (2000) and Sy and Taylor, (2000) which has been implemented in the smcure (Chao et al., , 2012) package in R (R Core Team, , 2018). In contrast to our scheme, $S_{T_{0}}$, rather than $S_{C}$, must be estimated. Thus, $S_{T_{0}}$ was estimated using a Cox model in which both covariates, $X_{(1)}$ and $X_{(2)}$, appear as predictors. The results, also shown in Table 1, are similar to those of our proposal when $\nu=0$. However, when $\nu=2$ (i.e., $S_{T_{0}}$ does not have the proportional hazards property), we see bias in the smcure estimates which does not disappear with increasing sample size. In particular, the bias manifests through $\widehat{\theta}_{(0)}$ and $\widehat{\theta}_{(2)}$; interestingly, $\widehat{\theta}_{(1)}$ is unaffected (i.e., the coefficient of $X_{(1)}$, the covariate which only enters the cure component).

As discussed in Section 4.2, the setting we consider in this paper lies within the theory of Hitomi et al., (2008), suggesting that the variability of the estimates will be higher when the true $S_{C}$ is used for estimation as per (3.7) than when $\widehat{S}_{C}$ is used as per (3.10), and, in turn, it will be lower still if the true cure labels are used as per (3.5). Although only the $\widehat{S}_{C}$ case is feasible, we display the standard errors for all three cases in Table 2, and the results for the other simulation scenarios are given in the Appendix (Section 8.4). Separately, but not shown, we have estimated $S_{C}$ parametrically using maximum likelihood, and considered additional scenarios where the true $S_{C}$ does not depend on $X$ (estimating $S_{C}$ both parametrically and non-parametrically); in all cases we have considered, the results are as predicted by Hitomi et al., (2008).

5.3 Selection procedure

We simulated data as per Section 5.1, but with four additional independent $\text{Normal}(0,1)$ variables, $X_{(3)}$, $X_{(4)}$, $X_{(5)}$, and $X_{(6)}$. These variables do not affect the cure probability, i.e., their $\theta$ coefficients are zero, but we set $\beta_{T_{0},(3)}=\beta_{C,(3)}=1$ so that $X_{(3)}$ affects other aspects of the data generating process; the $\beta_{T_{0}}$ and $\beta_{C}$ coefficients for $X_{(4)}$, $X_{(5)}$, and $X_{(6)}$ are all zero. The bias and standard errors for this setup can be found in the Appendix (Section 8.4), and, expectedly, they are a little larger than the scenarios with two covariates. However, here, we focus on the results of variable selection using the alasso where $\lambda$ was chosen by minimizing the cross-validation error (see Appendix, Section 8.3 for details of this algorithm). Let $\lambda^{\text{CVE}}$ denote this $\lambda$ value, and define the following commonly-used metrics for assessing the selection performance: $\text{C}=\sum_{j=4}^{6}\mathbbm{1}(\widehat{\theta}_{\lambda^{\text{CVE}},(j)}=0)$, the number of coefficients correctly set to zero, $\text{IC}=\sum_{j=1}^{2}\mathbbm{1}(\widehat{\theta}_{\lambda^{\text{CVE}},(j)}=0)$, the number of coefficients incorrectly set to zero, and $\text{DF}=\sum_{j=0}^{6}\mathbbm{1}(\widehat{\theta}_{\lambda^{\text{CVE}},(j)}>0)$, the model degrees of freedom (i.e., the number of non-zero parameters); in our setup, for the oracle model, C $=4$, IC $=0$, and DF $=3$. These metrics, averaged over simulation replicates, are shown for the alasso in Table 3. In line with what we expect from Theorem 4.3, IC approaches zero as the sample size increases, while C approaches four; again the results are unaffected by $\nu$.

6 Data analysis: colon cancer

We consider a colon cancer dataset (contained in the survival package in R) which was collected as part of a well-known national intergroup randomized controlled trial (involving Eastern Cooperative Oncology Group, the North Central Cancer Treatment Group, the Southwest Oncology Group, and the Mayo Clinic). The aim of the study was to investigate the efficacy of the drugs levamisole and 5FU for the treatment of colon cancer following surgery; relapse-free survival was the outcome variable of interest, i.e., time from randomization until the earlier of cancer relapse or death. In total, 929 patients with stage C disease enrolled during the period March 1984 to October 1987, with a maximum follow-up time of nine years. These patients were randomized to the following treatments: observation (control / reference group), levamisole, and a combined treatment of levamisole and 5FU. In addition to the treatment variable, a variety of binary covariates were recorded (reference categories are shown first): days since surgery, $\{\leq 20,>20\}$; sex, $\{\text{female},\text{male}\}$; obstruction of colon by tumour, $\{\text{no},\text{yes}\}$; adherence to nearby organs, $\{\text{no},\text{yes}\}$; depth of invasion, $\{\text{submucosa or muscular layer},\text{serosa}\}$; positive lymph nodes, $\{\leq 4,>4\}$. Furthermore, the age of the patient was recorded, and we use a mean-centered version (mean age is 59.75 years). See Moertel et al., (1990) for further details. This dataset is a candidate for cure analysis based on its Kaplan-Meier (KM) curve which has a clear plateau at approximately 40% (Figure 1).

We estimate the cure parameters using our proposed procedure and, for comparison, apply smcure. We use a Cox model with all covariates for $\widehat{S}_{C}$ in our approach and for $\widehat{S}_{T_{0}}$ in smcure, and a logistic regression cure model in both cases; confidence intervals and p-values are produced using bootstrapping. We also carry out variable selection using the adaptive lasso where covariates are standardized for variable selection, after which the estimates are transformed back to correspond to the original scale. See Table 4.

First we consider the the unpenalized estimates. The effect of the levamisole treatment does not significantly increase the cure probability (compared with a patient receiving no treatment), while the combination of levamisole with 5FU does; indeed, the odds of being cured for this latter treatment are $2.56$ ($=\exp(0.94)$) with 95% confidence interval given by $(1.48,5.64)$. The effect of all other covariates is to reduce the cure probability, albeit sex is not statistically significant, and obstruction and adherence are only just significant at the 10% level. The results for smcure are broadly similar, apart from the fact that adherence is statistically significant. Now, turning to the penalized alasso estimates, several coefficients have been set to zero, and the retained variables are those with smaller p-values from the unpenalized model. The regularization paths for standardized cure coefficients (i.e., those corresponding to standardized covariates) provide useful information on the relative importance of each covariate; these are shown in Figure 2. We can see immediately that the Lev+5FU treatment is one of the most important features. The number of positive lymph nodes is also highly important, and the presence of more than four such lymph nodes reduces the chance of cure. Next, the timing of surgery and depth of the tumour have similar importance, followed by the presence of an obstruction.

7 Discussion

We have proposed an IPCW-likelihood-based estimation procedure for cure regression models; elsewhere IPCW has been advocated by Gerds et al., (2017) as a device for producing straightforward estimators in complex survival data. In contrast to current cure estimation procedures in the literature, our assumptions are placed on $S_{C}$ while $S_{T_{0}}$ is completely unspecified. Although we have considered a Cox model estimator for $S_{C}$ in the examples in this article, any arbitrarily flexible model can be used in practice as this simply “plugs in” to the likelihood function given in (3.10) without any added complexity to the estimation procedure. Moreover, our asymptotic results still hold once the estimator, $\widehat{S}_{C}$, permits an iid representation (and we have given many common examples in the Appendix, Section 8.2).

Except for the case of a fully nonparametric approach like in Xu and Peng, (2014) (which suffers from the curse of dimensionality), existing cure regression models impose assumptions on both the cure proportion and the law of the susceptible individuals, without satisfactory model diagnosis (besides ad-hoc efforts). In our approach, one can first use standard diagnosis procedures to validate the censoring model as this is identifiable from the observed data directly. For example, one could assess the proportional hazards assumption for $S_{C}$ using the test due to Grambsch and Therneau, (1994) which is implemented in the cox.zph function in the survival package in R. (Although not shown, this test supported the proportional hazards assumption in the application considered in Section 6.) Next one could consider model diagnostics for the cure regression. Furthermore, note that our theory is not limited to the logistic model choice used in our applications, and, more generally still, the functional form of the cure regression model $\pi(X^{\top}\theta)$ could itself be estimated (e.g., in a similar manner to Amico et al., (2018) who extended the EM approach in this way). Goodness-of-fit for the cure model and estimation of its functional form are beyond the scope of the current article.

Although the extension to penalized estimation is straightforward and computationally efficient in our setting, we note that penalized selection in the existing EM setting was considered by Liu et al., (2012). However, those authors remark on the computational intensiveness of the procedure (penalized estimation and tuning parameter selection are embedded in EM). Furthermore, their approach is limited to a Cox PH model for $S_{T_{0}}$, and lacks asymptotic theory.

Lastly, our focus has been on modelling the cure probability without specifying a model for the latency distribution. However, as pointed out by an anonymous reviewer, the latter might also be of interest in applications. Indeed, from the fact that $S_{T}(t\mid X)=\{1-\pi(X)\}S_{T_{0}}(t\mid X)+\pi(X)$, we have $S_{T_{0}}(t\mid X)=\{S_{T}(t\mid X)-\pi(X)\}/\{1-\pi(X)\}$ which immediately provides $\widehat{S}_{T_{0}}$ given $\widehat{\pi}$ and $\widehat{S}_{T}$. Of course, $S_{T}$ can be modelled in an arbitrarily flexible way as it is based directly on the observable data, and any of the standard survival models can simply “plug in” just as for the $S_{C}$ model in our framework.

8 Appendix

8.1 Proofs

Lemma 8.1

Let $r(Y,X)$ be an integrable real-valued function. Under conditions (2.1) and (2.2),

[TABLE]

Proof of Lemma 8.1.

First, we have that

[TABLE]

where $\mathbb{E}(1-B\mid C,T_{0},X)=\mathbb{E}(1-B\mid X)=1-\pi(X)$ follows from (2.2), and $\mathbb{E}[\mathbbm{1}(T_{0}\leq C)\mid T_{0},X]=S_{C}(T_{0}-\mid X)$ follows from (2.1). Thus,

[TABLE]

as required. (The first equality holds since only the $\Delta=1$ case contributes where $Y=T_{0}$.)

Lemma 8.2

Let $B$ be a Bernoulli random variable, $T_{0}$ a nonnegative random variable and let $T=T_{0}$ if $B=0$ and $T=\infty$ if $B=1$. Then

[TABLE]

Proof of Lemma 8.2.

By elementary properties of the conditional independence

[TABLE]

Next,

[TABLE]

The result follows from the fact that $T$ is completely determined by $B$ and $T_{0}$.

Lemma 8.3

Let

[TABLE]

If the cure regression model is identifiable,

[TABLE]

Proof of Lemma 8.3.

Let $\pi^{*}(X)=\mathbb{E}\{B_{i}^{*}(S_{C})\mid X\}$, which by construction lies between 0 and 1. Then, since for any $u>0$, $\log(u)\leq u-1$, we deduce

[TABLE]

and the inequality is strict unless $[\pi(X)-\pi^{*}(X)]\omega(X)=0$ almost surely. If the cure model is identifiable, this cannot happen and thus $\mathbb{E}\{\ell^{*}(\theta)\}<\mathbb{E}\{\ell^{*}(\theta_{0})\},$ $\forall\theta\neq\theta_{0}.$

Proof of Lemma 4.1.

Define the event

[TABLE]

On the set $\mathcal{E}_{n}$, given a measurable function $\phi(\cdot),$ we can write

[TABLE]

where $\widehat{R}_{i}=\widehat{S}_{C}(Y_{i}-\mid X_{i})-S_{C}(Y_{i}-\mid X_{i}).$ By the preservation of the Glivenko-Cantelli property for classes of functions (see Theorem 3 of van der Vaart and Wellner, (2000)), the set of logit transformations of $\pi(\theta)$ is a $\mathbb{P}_{X}-$Glivenko-Cantelli class of functions of $X$ with constant envelope, provided Assumption 4.2.2 and 4.2.4 hold true. The statement follows from the uniform law of large numbers for the empirical process indexed by the set of functions

$\left\{\omega(\cdot)S^{-2}_{C}(\cdot\;-\mid\cdot)|\log\{\pi(\cdot;\theta)/[1-\pi(\cdot;\theta)]\}:\theta\in\Theta\right\}$ with constant envelope, and condition (4.12) which implies $\sup_{1\leq i\leq n}|\widehat{R}_{i}|=o_{\mathbb{P}}(1)$ and $\Pr(\mathcal{E}_{n})\rightarrow 1$.

Proof of Theorem 4.1.

We apply Theorem 5.7 of van der Vaart, (2000). First, by construction and Assumption 4.2.3, $\theta_{0}$ is a well-separated maximum of the map $\theta\mapsto\mathbb{E}\{\ell^{*}(\theta)\},$ with $\ell^{*}(\theta)$ defined in equation (3.7). See the inequality of Lemma 8.3. Next, let us note that

[TABLE]

Since

[TABLE]

by our assumptions and the uniform law of large numbers for the empirical process,

[TABLE]

The consistency of $\widehat{\theta}$ follows from Theorem 5.7 of van der Vaart, (2000).

Proof of Theorem 4.2.

By the definition of $\widehat{\theta}$ we have the identity $\widehat{\ell}^{*}(\theta)$

[TABLE]

where $\overline{\theta}$ is some point on the segment between $\widehat{\theta}$ and $\theta_{0}$. Using the definition in equation (3.6) and the short notation $\pi_{i}(\theta)=\pi(X_{i};\theta)$ and $\partial\pi_{i}(\theta)/\partial\theta=\partial\pi(X_{i};\theta)/\partial\theta$, we can write

[TABLE]

and thus, by Assumption 4.5 we have

[TABLE]

Next, for any $\theta$ we have

[TABLE]

Here, $\partial^{2}\pi_{i}(\theta)/\partial\theta\partial\theta^{\top}$ denotes the matrix of the second order partial derivatives of $\pi(X_{i};\theta)$ with respect to $\theta$. By the preservation of the Glivenko-Cantelli property for classes of functions (see Theorem 3 of van der Vaart and Wellner, (2000)), and Assumptions 4.2.2 and 4.2.4, the sets

[TABLE]

and

[TABLE]

with $\mathcal{F}_{1,kl}$ and $\mathcal{F}_{2,kl},$ $0\leq k,l\leq p,$ defined in Assumption 4.4.3, are $\mathbb{P}_{X}-$Glivenko-Cantelli classes of functions of $X$ with integrable envelope. By Lemma 4.1, deduce that

[TABLE]

On the other hand, the uniform law of large numbers for the empirical process yields

[TABLE]

Lebesgue’s Dominated Convergence Theorem implies

[TABLE]

and

[TABLE]

Gathering facts, deduce

[TABLE]

and this completes the justification of the representation of $\widehat{\theta}-\theta_{0}$. The convergence in law of $\sqrt{n}(\widehat{\theta}-\theta_{0})$ is a direct consequence of this representation.

An alternative way to obtain (8.1) is to require that condition (4.1) holds true and impose more regularity on the regression functions in the cure regression model. More precisely, assume that, in addition to (4.1), there exists an integrable function $C(X)$ and a constant $a>0$ such that, for any $0\leq k\leq l\leq p$ and any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

and

[TABLE]

Then property (8.1) follows using arguments as in the proof of Lemma 4.1.

Proof of Theorem 4.3.

We follow along the lines of the proof of Theorem 4 in Zou, (2006). First consider the asymptotic normality part. Let $\theta=\theta_{0}+\mathbf{u}n^{-1/2}$ and define

[TABLE]

Let $\widehat{\mathbf{u}}=\arg\max_{\mathbf{u}}\Gamma_{n}(\mathbf{u}),$ such that $\widehat{\mathbf{u}}=\sqrt{n}(\widehat{\theta}_{\lambda}-\theta_{0}).$ By Taylor expansion

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The behaviour of $A_{1}^{(n)}$, $A_{2}^{(n)}$ and $A_{3}^{(n)}$ can be derived from the proof of Theorem 4.2. On the other hand, for $A_{4}^{(n)}$, by the same arguments as in the proof of Theorem 4 in Zou, (2006), we have

[TABLE]

in probability. To summarize, for every $\mathbf{u}$, $\Gamma_{n}(\mathbf{u})-\Gamma_{n}(\mathbf{0})$ converges in probability to the function

[TABLE]

where

[TABLE]

The iid representation, and hence the asymptotic normality for $\mathbf{u}_{\mathcal{A}}$, follows. Moreover, we deduce that $\widehat{\mathbf{u}}_{\mathcal{A}^{c}}\rightarrow 0,$ in probability.

Next, we investigate the consistency part. It remains to show that $\forall j^{\prime}\not\in\mathcal{A}$, we have $\Pr(j^{\prime}\in\mathcal{A}_{n})\rightarrow 0.$ Let us fix arbitrarily $j^{\prime}\in\mathcal{A}^{c}$ and consider the event $\{j^{\prime}\in\mathcal{A}_{n}\}$. Let $\mathbf{e}_{j^{\prime}}\in\mathbb{R}^{p+1}$ be the vector with zero components, except the $j^{\prime}$th component that is equal to 1. By the Karush-Kuhn-Tucker optimality conditions, we necessarily have

[TABLE]

Thus

[TABLE]

Next, by Taylor expansion, we can decompose

[TABLE]

with

[TABLE]

and

[TABLE]

with $\overline{\theta}_{\lambda}$ between $\widehat{\theta}_{\lambda}$ and $\theta_{0}$. By arguments that we have already used, we obtain

[TABLE]

Meanwhile, since $\sqrt{n}(\widehat{\theta}-\theta_{0})=O_{\mathbb{P}}(1)$ and $\lambda n^{(\gamma-1)/2}\rightarrow\infty$,

[TABLE]

Thus, $\Pr(j^{\prime}\in\mathcal{A}_{n})\rightarrow 0$ which completes the proof.

8.2 Uniform convergence and iid representations: examples

In this section we review several approaches for estimating the conditional law of the censoring time used for IPCW in such way that Assumption 4.3 and 4.5 hold true. Namely, we consider nonparametric estimators, such as the Kaplan-Meier and the conditional Kaplan-Meier estimators, and existing estimators in semiparametric models such as proportional hazards, proportional odds and transformation models. For all these models, we propose a guideline to derive iid representations under our assumptions from the existing asymptotic results.

Without loss of generality, we consider the case of a real-valued function $\varphi(X)$. Recall that we are interested in functions $\varphi$ such that $\omega(x)=0$ implies $\varphi(x)=0.$ In particular this implies

[TABLE]

For vector-valued functions it suffices to apply the results presented below for each component. For simplicity, we assume that the sample space $\mathcal{X}$ is a subset of a finite-dimensional space. Moreover, we assume $\Pr(T=C)=0$ and there exists $\tau\in\mathbb{R}$ such that $\Pr(T_{0}>\tau)=0$ and $\inf_{x\in\mathcal{X},\;\omega(x)>0}S_{C}(\tau\mid x)>0.$

8.2.1 Kaplan-Meier inverse probability weighting

When the law of the censoring variable $C$ does not depend on the covariates, the survivor function $S_{C}$ can be estimated by the Kaplan-Meier estimator

[TABLE]

Here $\widehat{H}(y-)$ is the left-sided limit of the estimate of marginal distribution function of the observed lifetimes $H(y)=\Pr(Y\leq y)$, $y\in\mathbb{R}$. Then the uniform convergence condition (4.12), and thus Assumption 4.3, is guaranteed by the uniform law of large numbers of Stute and Wang, (1993).

For Assumption 4.5, one can use the iid representation for Kaplan-Meier integrals, as stated in Theorem 1.1 of Stute, (1996). More precisely, for each $y\in\mathbb{R},$ define

[TABLE]

and

[TABLE]

with $H^{11}(x,y)=\Pr(X\leq x,Y\leq y,\Delta=1).$ Note that $H^{11}(x,y)=H^{11}(x,\tau)$, $\forall y>\tau$. Then Assumption 4.5 holds with

[TABLE]

8.2.2 Conditional Kaplan-Meier inverse probability weighting

The conditional Kaplan-Meier estimator, also called Beran estimator (Beran, , 1981) is defined as

[TABLE]

where

[TABLE]

Here, $b_{n}$ is a bandwidth sequence and $K(\cdot)$ is a multivariate kernel function. The uniform law of large numbers for this conditional Kaplan-Meier estimator was established in Corollary 2.1 of Dabrowska, (1989) under some regularity conditions on the density of $X$ and the functions $x\mapsto\Pr(Y\leq y,\Delta=j\mid X=x),$ $j=0,1$, $y\leq\tau$. To apply that corollary in our framework, one has to define a set on which the density of the covariate vector stays away from zero and to take the weight function $\omega(\cdot)$ in the definition of our estimators $\widehat{\theta}$ and $\widehat{\theta}_{\lambda}$ equal to zero outside this set.

The iid representation for Kaplan-Meier integrals was extended to conditional Kaplan-Meier integrals, see Lopez, (2011). However, this purely nonparametric approach suffers from the curse of dimensionality when the sample space $\mathcal{X}$ is multidimensional. Stronger regularity assumptions and high-order kernels are needed in such cases. For simplicity, following Lopez, (2011), assume that $C\perp X\mid Z$ where $Z=g(X)\in\mathbb{R}$ with $g(\cdot)$ a given function. For instance, $g(X)$ could be a component of $X$, or a given linear combination of components of $X$. The case where $g(\cdot)$ is known up to a finite-dimensional parameter which has to be estimated could be also considered, but would introduce an additional term in the iid representation in Assumption 4.5 which takes into account the estimation of $g(\cdot)$. For simplicity, herein we assume that $g(\cdot)$ is given. Moreover, for some small $\delta>0$, the weight function $\omega(\cdot)$ vanishes outside the set $\mathcal{X}_{\delta}=\{x\in\mathcal{X}:f_{g}(g(x))\geq\delta+\inf_{z}f_{g}(z)\}$. Here, $f_{z}$ denotes the density of $Z=g(X),$ which we assume exists and satisfies some differentiability conditions. Then, the Beran estimator with $Z_{i}=g(X_{i}),$ $1\leq i\leq n,$ is defined as

[TABLE]

where

[TABLE]

Here, $a_{n}$ is a bandwidth sequence converging to zero as $n$ tends to infinity, and $K(\cdot)$ is a univariate kernel function. If the bandwidth $a_{n}$ satisfies $\log(n)n^{-1}a_{n}^{-3}\rightarrow 0$ and $na_{n}^{4}\rightarrow 0,$

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

Here, $H_{0}(y\mid z)=\Pr(Y\leq y,\delta=0\mid g(X)=z)$ and $H^{11}(x,y\mid z)=\Pr(X\leq x,Y\leq y,\delta=1\mid g(X)=z).$

8.2.3 Semiparametric models

In this section we present a general method for guaranteeing Assumption 4.5. This method can, for example, be applied to the common semiparametric models used in survival analysis. (For all the models we mention, condition (4.12) is obviously satisfied, and thus Assumption 4.3.) Note that what we present here is a fuller version of the sketch which appears in Section 4.2 of the main paper, and, thus, there is some repetition which we maintain to improve readability herein. Assume that the conditions of Lemma 4.1 hold true and $\varphi(\cdot)$ is bounded. We can write

[TABLE]

where for any function $S$ depending on $y$ and $x$ that is càdlàg in $y$,

[TABLE]

Then Assumption 4.5 can be guaranteed in two steps. First, use the equicontinuity of Donsker classes and transform the sum with respect to $1\leq i\leq n$ to an expectation with respect to a generic triplet $(Y,\Delta,X)$ given the sample. Next, use the iid representation of $\widehat{S}_{C}(y-\mid x)$ which can be derived in common semiparametric models. For the first step, in the following, we introduce a general class of functions for which we prove the Donsker property. For several common semiparametric models and estimators, it can be shown that $q(\cdot,\cdot,\cdot;S_{C})$ and $q(\cdot,\cdot,\cdot;\widehat{S}_{C})$ belong to this class with probability tending to 1. For other models, it may be necessary to define alternative Donsker classes to be used in the first step of the method we propose herein.

Let us suppose that $X$ is composed of continuous and discrete components, that is $X=(X_{d}^{\top},X_{c}^{\top})^{\top}\in\mathcal{X}_{d}\times\mathcal{X}_{c}\subset\mathbb{R}^{p_{d}}\times\mathbb{R}^{p_{c}}.$ Then, each vector $x$ in the support of $X$ can be split in the subvectors $x_{d}\in\mathcal{X}_{d}$ and $x_{c}\in\mathcal{X}_{c}.$ For simplicity, assume that the support $\mathcal{X}_{d}$ is finite.

Let $BV_{M}[0,\tau]$ be the set of all real-valued càdlàg functions defined on $[0,\tau]$ with the total variation bounded by $M$. Let $\mathcal{C}_{M}^{\alpha}(\mathcal{X}_{c})$ be the set of all continuous functions $f:\mathcal{X}_{c}\rightarrow\mathbb{R}$ with $\|f\|_{\alpha}\leq M,$ where $\|\cdot\|_{\alpha}$ is the usual uniform norm defined on the class of functions with uniformly bounded partial derivatives up to order $\underline{\alpha}$ (the greatest integer smaller than $\alpha$) and whose highest partial derivatives are Lipschitz of order $\alpha-\underline{\alpha}$. See chapter 2.7.1 in van der Vaart and Wellner, (1996). We follow Lopez, (2011) and, for $\alpha,M>0$, define the following class of functions defined on $[0,\tau]\times\mathcal{X}$ :

[TABLE]

where for any vector $\mathbf{k}=(k_{1},\ldots,k_{p_{c}})$ of $p_{c}$ integers, $D^{\mathbf{k}}$ is the differential operator

[TABLE]

The next result shows that $\mathcal{G}$ is a general Donsker class. For this it suffices to show that $\mathcal{G}$ has finite bracketing integral, which here is tantamount to $\int_{0}^{1}\sqrt{\log N_{[\;]}(\varepsilon,\mathcal{G},L_{2})}\;d\varepsilon<\infty,$ where $N_{[\;]}(\varepsilon,\mathcal{G},L_{2})$ is the $\varepsilon-$bracketing number of the class $\mathcal{G}$ with respect to the $L_{2}-$norm.

Lemma 8.4

Assume that $\mathcal{X}_{c}$ is a bounded, open and convex subset of $\mathbb{R}^{p_{c}}$ and $\mathcal{X}_{d}$ is finite. If $M<\infty$ and $\alpha>p_{c}$, $\mathcal{G}$ is a Donsker class.

Proof of Lemma 8.4.

We provide a sketch proof. The idea is to show that

[TABLE]

for some constant $C$. Since $\int_{0}^{1}(1/\varepsilon)^{p_{c}/\alpha}\log^{1/2}(1/\varepsilon)d\varepsilon<\infty$ when $p_{c}/\alpha<1$, this will guarantee that $\mathcal{G}$ has finite bracketing integral. With a finite set $\mathcal{X}_{d}$, by the permanence of the Donsker property (see Theorem 2.10.6 in van der Vaart and Wellner, (1996)), it suffices to show property (8.4) separately for each value $x_{d}\in\mathcal{X}_{d}.$ In other words, it suffices to consider that $x=x_{c}$.

First, for any $t\in[0,\tau],$ let $\mathcal{G}_{t}=\{x\mapsto h(t,x):h\in\mathcal{G}\}.$ Then, by Corollary 2.6.2 of van der Vaart and Wellner, (1996), we have that $\log N_{[\;]}(\varepsilon,\mathcal{G}_{t},L_{2})\leq C_{1}\varepsilon^{-p_{c}/\alpha},$ for each $t\in[0,\tau]$, where $C_{1}$ is a constant depending only on $M$, $\alpha,$ $\text{diam}(\mathcal{X}_{c})$ and $p_{c}$ (and is independent of $t$). Meanwhile, $\log N_{[\;]}(\varepsilon,BV_{M}[0,\tau],L_{\infty})\leq C_{2}(1/\varepsilon)\log(1/\varepsilon),$ where $C_{2}$ is a constant depending only on $M$ and $\tau$. This is a simple bound for the bracketing number, and can be obtained by constructing brackets which are piecewise constant on a regular grid. See also the beginning of the proof of Theorem 2.7.5 in van der Vaart and Wellner, (1996). In the following, without loss of generality, we consider $M=1$. To prove (8.3), it suffices to reconsider the proof of Theorem 2.7.1 (and Corollary 2.7.2) of van der Vaart and Wellner, (1996). More precisely, in the the proof of Theorem 2.7.1, replace the entries of the matrices $Af$ by functions of $t$ with total variation bounded by $M$. In van der Vaart and Wellner, (1996), the entries of $Af$ are defined for $f$ ranging over $\mathcal{C}_{1}^{\alpha}(\mathcal{X}_{c})$ and having the rows built using the values of the partial derivatives of $f$, up to order equal to the greatest integer strictly smaller than $\alpha$, discretized on a grid with the mesh controlled by $\varepsilon$. In our case, such entries may depend also on $t$ and be taken as a set of brackets covering $BV_{M}[0,\tau]$. Deduce a set of brackets that cover $\mathcal{G}$ and the cardinality of the set of brackets is of order

[TABLE]

Taking the logarithm of the above, we recover the bound $(1/\varepsilon)^{2p_{c}/\alpha}\log(1/\varepsilon)$ for the order of the bracketing entropy of $\mathcal{G}$.

As a first step for guaranteeing Assumption 4.5, using the properties of the model and estimator considered, one can check that

[TABLE]

in order to deduce, using the permanence properties of Donsker classes, that $q(\cdot,\cdot,\cdot;S_{C})$ and $q(\cdot,\cdot,\cdot;\widehat{S}_{C})$ belong to a Donsker class with probability tending to 1. By the asymptotic equicontinuity of Donsker classes (see van der Vaart and Wellner, (1996)), and our Lemma 8.1,

[TABLE]

where

[TABLE]

and $F(x,t)=\Pr(X\leq x,T_{0}\leq t)$. In the second step, it remains to show that

[TABLE]

using the properties of the model considered for $S_{C}$.

In the following we study conditions (8.4) and (8.5) in several common (semi)parametric models in survival analysis. See Guo and Zeng, (2014) for an illuminating survey on semiparametric models. In these models, the properties of $\widehat{S}_{C}$ are usually derived from the properties of the conditional cumulative hazard function of $C$. Note that by the Duhamel identity (Gill and Johansen, , 1990) we have

[TABLE]

where $\Lambda_{C}(\cdot\mid x)$ is the conditional cumulative hazard function of $C$ given $X=x$ and $\widehat{\Lambda}_{C}(\cdot\mid x)$ is the estimator of $\Lambda_{C}(\cdot\mid x)$ in the model for $S_{C}$. When $S_{C}(\cdot\mid x)$ is continuous,

[TABLE]

Parametric models. One can achieve flexible modeling for $S_{C}$ using parametric models, such as the Weibull model, where the parameters are replaced by functions of $X$ depending on some unknown vector of coefficients $\beta_{C}$. Then, in general, all the survivor functions in the model satisfy the regularity conditions defining the class $\mathcal{G}$, so that condition (8.4) is automatically met. Condition (8.5) follows by Taylor expansion and the asymptotic linear expansions of $\widehat{\beta}_{C}$, the estimator of $\beta_{C}$ considered. Such asymptotic linear expansions for $\widehat{\beta}_{C}$, based on the so-called influence functions, are available for all common estimators $\widehat{\beta}_{C}$.

Cox’s proportional hazard model. In this case

[TABLE]

where $\Lambda_{0,C}(\cdot)$ is the so-called baseline cumulative hazard function. Clearly, $\Lambda_{C}(t\mid x)$ belongs to $\mathcal{G}$. Any estimate of $\Lambda_{C}(t\mid x)$ belongs to $\mathcal{G}$ once $\Lambda_{0,C}(\cdot)$ is estimated by a function with total variation bounded by some suitable $\widetilde{M}$ depending on $M$, the compact parameter set for $\beta_{C}$, and the bounded $\mathcal{X}$. Thus, condition (8.4) is easily granted. Next, for any consistent estimator $\widehat{\beta}_{C}$ and uniformly consistent $\widehat{\Lambda}_{0,C}(\cdot)$ we can write

[TABLE]

In the case where $\widehat{\beta}_{C}$ is the maximum partial likelihood estimator, and $\widehat{\Lambda}_{C}$ is the associated Breslow estimator, condition (8.5) follows by the asymptotic results of Andersen and Gill, (1982).

Transformation model. We consider the class of transformation models investigated by Zeng and Lin, (2006), where

[TABLE]

where $G(\cdot)$ is a given smooth, strictly increasing transformation function with $G(0)=0,$ $G^{\prime}(0)>0$ and $G(\infty)=\infty$. Box-Cox transformations, and log transformations are possible examples. Zeng and Lin, (2006) extend the partial likelihood idea and the Breslow estimator from Cox’s model and introduce the estimators $\widehat{\beta}_{C}$ and $\widehat{\Lambda}_{0,C}$ for which they derive a Gaussian limit. Since

[TABLE]

we can easily adapt the arguments used for the case of Cox’s model and, based on definition of $\widehat{\beta}_{C}$ and $\widehat{\Lambda}_{0,C}(\cdot)$ and the asymptotic results of Zeng and Lin, (2006), guarantee conditions (8.4) and (8.5).

Proportional odds model. In such a model we have

[TABLE]

with $H(\cdot)$ some càdlàg function with $H(0)=0$. See Murphy et al., (1997). Clearly, $S_{C}\in\mathcal{G}$ and $\widehat{S}_{C}\in\mathcal{G},$ with $\widehat{S}_{C}(t\mid x)=\exp(-x^{\top}\widehat{\beta}_{C})/\{\widehat{H}(t)+\exp(-x^{\top}\widehat{\beta}_{C})\}$, once $\widehat{H}(t)\in BV_{\widetilde{M}}[0,\tau]$ for some suitable $\widetilde{M}$. For the required iid representation, we can use the maximum likelihood estimators for $\beta_{C}$ and $H(\cdot)$, linearize the expression of $S_{C}$ with respect to the parameters and use the asymptotic representations of $\sqrt{n}(\widehat{H}-H)$ and $\sqrt{n}(\widehat{\beta}_{C}-\beta_{C})$ established in the proof of the Theorem 2.2 in Murphy et al., (1997). That result is established under an additional condition which, in our setup, means that $\Pr(T_{0}=\tau)=\Pr(T_{0}\geq\tau)>0$. This additional technical constraint is quite usual in the cure regression literature. See, for instance, Lu, (2008) and Fang et al., (2005).

8.3 Implementation of the alasso

8.3.1 Optimization

A variety of algorithms have been implemented for solving non-differentiable lasso problems, e.g., quadratic programming (Tibshirani, , 1996), least angle regression (LARS) (Efron et al., , 2004), and co-ordinate descent (Friedman et al., , 2007). However, we prefer the use of a differentiable penalty since standard gradient-based optimization procedures can then be utilized. Therefore, we propose the use of

[TABLE]

where $a_{\epsilon}(x)=(x^{2}+\epsilon^{2})^{1/2}-\epsilon$ is an extension of the absolute value function such that $\lim_{\epsilon\rightarrow 0}a_{\epsilon}(x)=|x|$, and which is differentiable for $\epsilon>0$. Clearly, smaller $\epsilon$ values bring the penalty closer to the alasso, but also bring (8.6) closer to being non-differentiable. In our work, we have found that $\epsilon=10^{-4}$ works well.

8.3.2 Tuning parameter selection

For the purpose of selecting the tuning parameter, $\lambda$, we consider cross-validation; in particular, we aim to minimize the $k$-fold cross-validation error. Since we have $\mathbb{E}\{B_{i}^{*}(S_{C})-\pi_{i}(\theta)\mid X_{i}\}=0$, we may define the error term, $B_{i}^{*}(\widehat{S}_{C})-\pi_{i}(\widehat{\theta})$. Then, for a partition $F_{1},\ldots,F_{K}$ of the set $\{1,\ldots,n\}$, the mean-squared error for the $j$th fold, $F_{j}$ is given by $\sum_{i\in F_{j}}\{B_{i}^{*}(\widehat{S}_{C})-\pi_{i}(\widehat{\theta}_{\lambda}^{-j})\}^{2}$ where $\widehat{\theta}_{\lambda}^{-j}$ is the penalized estimate with the $j$th fold removed. Thus, the $k$-fold cross-validation error is

[TABLE]

where we use $k=10$ as is standard in practice. Minimizing (8.7) with respect to $\lambda$ can be achieved by profiling over a range of $\lambda$ values or by using a one-dimensional optimizer, e.g., golden search; we will define $\lambda^{\text{CVE}}$ to be the minimizer of (8.7). One might also consider a BIC-type criterion of the form $-2\widehat{\ell}^{*}(\widehat{\theta}_{\lambda})+\text{dim}(\widehat{\theta}_{\lambda})\log(n)$. We have also tested this in our simulation study where it produced very similar results to the cross-validation approach; thus, in Section 5, we only present cross-validation.

8.4 Additional simulation results

Our simulations setup, described in Section 5.1 of the main paper, comprises 24 scenarios, whereas, only the results for the 6 scenarios where $\pi_{m}=0.4$ and $\rho=0.1$ were reported on in detail. Therefore, here, all results are displayed in Table 5. The estimates improve with the sample size (as seen in the main paper), disimprove with increased censoring proportion (controlled via $\rho$), and improve with an increased cure proportion, $\pi_{m}$. The estimates from our approach and smcure are similar, with smcure being slightly more efficient; of course, in the case where $\nu=2$ (non-PH $S_{T_{0}}$) the smcure estimates are biased. In Section 5.3, we described a simulation setup with six covariates. However, our focus in that setting was on the performance of our variable selection procedure. For completeness, Table 6 displays the results of estimation in that setting which are quite similar to the case with two covariates (albeit with slightly increased bias and variability as expected).

As suggested by an anonymous reviewer, it is also of interest to consider simulation setups with binary covariates. Thus, we considered two further setups: one where $X_{(1)}$ and $X_{(2)}$ are both binary, and another where $X_{(1)}$ is continuous while $X_{(2)}$ is binary. The results, shown in Table 7, are numerically very close to those of Table 5. We also considered these additional setups with six covariates, but, to avoid repetition, we do not present the results as they are very close to those of Table 6.

Theorem 4.2 establishes the asymptotic normality of our proposed estimator. However, as mentioned in Section 4.2 of the main paper, the asymptotic covariance can be difficult to estimate in general, and, therefore, we suggest using bootstrapping. Table 8 shows the empirical coverage for 95% confidence intervals constructed using bootstrapping with 399 replicates for the simulated data; we find that the empirical coverage is close to the nominal level. The same is true for the cases with six covariates and binary covariates, but we do not present them here for the sake of brevity.

Table 9 displays the results of the adaptive lasso variable selection procedure for all simulation scenarios, as well as results for the lasso procedure. The IC values tend towards zero as the sample size increases for both the lasso and alasso. However, the lasso tends to select a more complex model than the alasso as indicated by the smaller C values and larger DF values; it is well known that the lasso exhibits this behaviour which is why alasso is preferred (Fan and Lv, , 2010). Overall, the alasso works well with C approaching the oracle value of four as the sample size increases. The results are unaffected by the value of $\nu$ as we might expect (since $S_{T_{0}}$ is unspecified in our approach), whereas, when $n=100$, increased censoring proportion, $\rho$, or decreased cure proportion, $\pi_{m}$, both lead to fewer variables being selected.

Table 10 displays the standard errors of the estimates (in the case with two continuous covariates) when estimation proceeds based on (i) knowing the true $S_{C}$ and using $B_{i}^{*}(S_{C})$, (ii) estimating $S_{C}$ and using $B_{i}^{*}(\widehat{S}_{C})$, and (iii) directly observing the cure labels, $B_{i}$. The theory of Hitomi et al., (2008) predicts that the variability will be sequentially lower over these three cases, and this is clearly the case from Table 10. Most importantly, there is efficiency gain rather than loss by estimating $S_{C}$.