Generalized G-estimation and Model Selection

TL;DR

This paper extends G-estimation methods for dynamic treatment regimes to handle non-additive effects and non-continuous outcomes, introducing a model selection criterion and demonstrating its effectiveness through simulations and depression treatment data.

Contribution

It introduces a generalized G-estimation approach using iteratively-reweighted least squares and develops an information criterion for blip model selection.

Findings

Effective G-estimation extension for complex models

Proposed model selection criterion improves blip model choice

Successful application to depression treatment data

Abstract

Dynamic treatment regimes (DTRs) aim to formalize personalized medicine by tailoring treatment decisions to individual patient characteristics. G-estimation for DTR identification targets the parameters of a structural nested mean model known as the blip function from which the optimal DTR is derived. Despite considerable work deriving such estimation methods, there has been little focus on extending G-estimation to the case of non-additive effects, non-continuous outcomes or on model selection. We demonstrate how G-estimation can be more widely applied through the use of iteratively-reweighted least squares procedures, and illustrate this for log-linear models. We then derive a quasi-likelihood function for G-estimation within the DTR framework, and show how it can be used to form an information criterion for blip model selection. These developments are demonstrated through application to a variety of simulation studies as well as data from the Sequenced Treatment Alternatives to Relieve Depression study.

Tables

Table 1. Table 1: Mean blip parameter estimates (standard errors) from 1000 simulation runs for log-linear outcome model via iteratively-reweighted least squares. True blip parameters ( ψ j ​ 0 , ψ j ​ 1 ) subscript 𝜓 𝑗 0 subscript 𝜓 𝑗 1 (\psi_{j0},\psi_{j1}) of ( 0.5 , − 0.5 ) 0.5 0.5 (0.5,-0.5) .

$n$	$P\left(Y=0\right)$	${\hat{\psi }}_{10}$ (SE)	${\hat{\psi }}_{11}$ (SE)	${\hat{\psi }}_{20}$ (SE)	${\hat{\psi }}_{21}$ (SE)
50	5%	0.485 (0.337)	-0.349 (0.260)	0.487 (0.320)	-0.437 (0.372)
	10%	0.511 (0.370)	-0.345 (0.295)	0.500 (0.348)	-0.449 (0.407)
	20%	0.523 (0.418)	-0.339 (0.343)	0.496 (0.388)	-0.445 (0.461)
100	5%	0.503 (0.209)	-0.427 (0.155)	0.497 (0.217)	-0.475 (0.235)
	10%	0.501 (0.226)	-0.428 (0.165)	0.496 (0.232)	-0.470 (0.247)
	20%	0.520 (0.248)	-0.419 (0.191)	0.504 (0.266)	-0.483 (0.287)
200	5%	0.501 (0.144)	-0.471 (0.101)	0.504 (0.153)	-0.486 (0.166)
	10%	0.501 (0.154)	-0.471 (0.109)	0.502 (0.163)	-0.486 (0.171)
	20%	0.510 (0.174)	-0.470 (0.124)	0.501 (0.186)	-0.484 (0.198)
500	5%	0.498 (0.089)	-0.485 (0.062)	0.502 (0.095)	-0.496 (0.102)
	10%	0.497 (0.095)	-0.485 (0.066)	0.505 (0.100)	-0.495 (0.110)
	20%	0.502 (0.107)	-0.485 (0.077)	0.500 (0.114)	-0.495 (0.120)

Table 2. Table 2: Model selection for a variety of sample sizes ( n 𝑛 n ). Numbers indicate proportion of 1000 simulation runs where the correct model was selected by the corresponding method. (F) and (B) denote forwards and backwards selection, respectively. Bold indicates the most successful approach for each setup.

$n$	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
50	$x_{11}$	0.285	0.276	0.303	0.297
	$x_{11},x_{12}$	0.232	0.238	0.137	0.160
	$x_{11},x_{12},x_{13}$	0.218	0.259	0.060	0.108
	$x_{21}$	0.225	0.216	0.263	0.248
	$x_{21},x_{22}$	0.203	0.213	0.119	0.145
	$x_{21},x_{22},x_{23}$	0.181	0.214	0.068	0.09
100	$x_{11}$	0.437	0.423	0.479	0.471
	$x_{11},x_{12}$	0.373	0.380	0.284	0.318
	$x_{11},x_{12},x_{13}$	0.372	0.419	0.181	0.254
	$x_{21}$	0.372	0.366	0.422	0.418
	$x_{21},x_{22}$	0.353	0.365	0.263	0.285
	$x_{21},x_{22},x_{23}$	0.357	0.381	0.150	0.205
200	$x_{11}$	0.533	0.530	0.606	0.601
	$x_{11},x_{12}$	0.569	0.573	0.484	0.519
	$x_{11},x_{12},x_{13}$	0.587	0.628	0.379	0.451
	$x_{21}$	0.479	0.477	0.581	0.577
	$x_{21},x_{22}$	0.505	0.507	0.439	0.453
	$x_{21},x_{22},x_{23}$	0.538	0.557	0.334	0.371

Table 3. Table 3: Stage 1 model selection ( n = 100 𝑛 100 n=100 ) when the stage 2 model is correctly specified (‘Correct’), selected by the corresponding approach (‘Recommended’), or an intercept-only model (‘Intercept’). Numbers indicate proportion of 1000 simulation runs where the correct model was selected by the corresponding method. (F) and (B) denote forwards and backwards selection, respectively. Bold indicates the most successful approach for each setup.

Selection	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
Correct	$x_{11}$	0.437	0.423	0.479	0.471
	$x_{11},x_{12}$	0.373	0.380	0.284	0.318
	$x_{11},x_{12},x_{13}$	0.372	0.419	0.181	0.254
Recommended	$x_{11}$	0.420	0.409	0.485	0.471
	$x_{11},x_{12}$	0.369	0.375	0.283	0.316
	$x_{11},x_{12},x_{13}$	0.372	0.417	0.175	0.250
Intercept	$x_{11}$	0.431	0.417	0.479	0.469
	$x_{11},x_{12}$	0.373	0.380	0.283	0.311
	$x_{11},x_{12},x_{13}$	0.373	0.411	0.167	0.245

Table 4. Table 4: Model selection ( n = 200 𝑛 200 n=200 ) for discrete outcome case. Entries are selection rates for the model indicated by the column heading. Bold indicates correct model selection rate.

True	Intercept	$x_{j1}$	$x_{j1},x_{j2}$	$x_{j1},x_{j2},x_{j3}$
$x_{11}$	0.017	0.663	0.169	0.151
$x_{11},x_{12}$	0.001	0.004	0.799	0.196
$x_{11},x_{12},x_{13}$	0.001	0.002	0.078	0.919
$x_{21}$	0.005	0.617	0.189	0.189
$x_{21},x_{22}$	0.001	0.006	0.836	0.157
$x_{21},x_{22},x_{23}$	0.000	0.002	0.076	0.921

Table 5. Table 5: Model selection ( n = 100 𝑛 100 n=100 ) with true blip parameters of 1, 0.5, or 0.1, under the setup described in section 3.2. (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

${\psi }_{jp}$	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
1	$x_{11}$	0.437	0.423	0.479	0.471
	$x_{11},x_{12}$	0.373	0.380	0.284	0.318
	$x_{11},x_{12},x_{13}$	0.372	0.419	0.181	0.254
	$x_{21}$	0.372	0.366	0.422	0.418
	$x_{21},x_{22}$	0.353	0.365	0.263	0.285
	$x_{21},x_{22},x_{23}$	0.357	0.381	0.150	0.205
0.5	$x_{11}$	0.267	0.259	0.215	0.208
	$x_{11},x_{12}$	0.160	0.169	0.071	0.086
	$x_{11},x_{12},x_{13}$	0.128	0.165	0.033	0.063
	$x_{21}$	0.232	0.229	0.219	0.214
	$x_{21},x_{22}$	0.162	0.177	0.077	0.089
	$x_{21},x_{22},x_{23}$	0.132	0.156	0.028	0.049
0.1	$x_{11}$	0.122	0.120	0.073	0.072
	$x_{11},x_{12}$	0.043	0.048	0.005	0.012
	$x_{11},x_{12},x_{13}$	0.021	0.037	0.000	0.002
	$x_{21}$	0.138	0.138	0.076	0.073
	$x_{21},x_{22}$	0.053	0.063	0.016	0.02
	$x_{21},x_{22},x_{23}$	0.030	0.045	0.003	0.003

Table 6. Table 6: Model selection ( n = 100 𝑛 100 n=100 ) with non-treatment covariates drawn from multivariate normal distributions with correlations of 0 (‘none’), 0.25 (‘medium’) and 0.5 (‘strong’), under the setup described in section 3.2. (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

$n$	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
None	$x_{11}$	0.437	0.423	0.479	0.471
	$x_{11},x_{12}$	0.373	0.380	0.284	0.318
	$x_{11},x_{12},x_{13}$	0.372	0.419	0.181	0.254
	$x_{21}$	0.372	0.366	0.422	0.418
	$x_{21},x_{22}$	0.353	0.365	0.263	0.285
	$x_{21},x_{22},x_{23}$	0.357	0.381	0.150	0.205
Medium	$x_{11}$	0.357	0.353	0.412	0.408
	$x_{11},x_{12}$	0.331	0.338	0.247	0.262
	$x_{11},x_{12},x_{13}$	0.325	0.346	0.150	0.173
	$x_{21}$	0.351	0.348	0.388	0.382
	$x_{21},x_{22}$	0.293	0.301	0.212	0.223
	$x_{21},x_{22},x_{23}$	0.289	0.307	0.129	0.141
Strong	$x_{11}$	0.351	0.347	0.375	0.374
	$x_{11},x_{12}$	0.297	0.299	0.194	0.196
	$x_{11},x_{12},x_{13}$	0.237	0.244	0.089	0.093
	$x_{21}$	0.316	0.313	0.337	0.335
	$x_{21},x_{22}$	0.234	0.236	0.151	0.152
	$x_{21},x_{22},x_{23}$	0.196	0.206	0.065	0.066

Table 7. Table 7: Non-aggregated selection ( n = 100 𝑛 100 n=100 ), under the setup described in section 3.2. (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

Stage 1	Stage 2	QICG (F)	QICG (B)	Wald (F)	Wald (B)
		Stage 1 Selection
$x_{11}$	$x_{21}$	0.437	0.422	0.492	0.481
	$x_{21},x_{22}$	0.443	0.429	0.476	0.467
	$x_{21},x_{22},x_{23}$	0.432	0.417	0.469	0.465
$x_{11},x_{12}$	$x_{21}$	0.383	0.388	0.287	0.328
	$x_{21},x_{22}$	0.367	0.373	0.288	0.321
	$x_{21},x_{22},x_{23}$	0.370	0.378	0.276	0.306
$x_{11},x_{12},x_{13}$	$x_{21}$	0.392	0.434	0.189	0.274
	$x_{21},x_{22}$	0.366	0.418	0.183	0.247
	$x_{21},x_{22},x_{23}$	0.359	0.406	0.171	0.241
		Stage 2 Selection
$x_{11}$	$x_{21}$	0.373	0.370	0.426	0.424
$x_{11},x_{12}$		0.375	0.366	0.418	0.411
$x_{11},x_{12},x_{13}$		0.367	0.362	0.421	0.419
$x_{11}$	$x_{21},x_{22}$	0.364	0.375	0.275	0.299
$x_{11},x_{12}$		0.351	0.364	0.255	0.278
$x_{11},x_{12},x_{13}$		0.345	0.356	0.258	0.278
$x_{11}$	$x_{21},x_{22},x_{23}$	0.358	0.386	0.162	0.215
$x_{11},x_{12}$		0.361	0.386	0.142	0.202
$x_{11},x_{12},x_{13}$		0.351	0.371	0.147	0.198

Table 8. Table 8: Model selection for a variety of sample sizes ( n 𝑛 n ), under the setup described in section 3.2 but with normal (rather than log-normal) errors ϵ ∼ N ​ ( 0 , 1 ) similar-to italic-ϵ 𝑁 0 1 \epsilon\sim N(0,1) . (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

$n$	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
50	$x_{11}$	0.377	0.365	0.482	0.475
	$x_{11},x_{12}$	0.388	0.399	0.294	0.341
	$x_{11},x_{12},x_{13}$	0.412	0.470	0.163	0.253
	$x_{21}$	0.297	0.286	0.375	0.364
	$x_{21},x_{22}$	0.291	0.300	0.205	0.238
	$x_{21},x_{22},x_{23}$	0.291	0.329	0.126	0.167
100	$x_{11}$	0.556	0.543	0.767	0.763
	$x_{11},x_{12}$	0.643	0.646	0.644	0.705
	$x_{11},x_{12},x_{13}$	0.796	0.842	0.513	0.652
	$x_{21}$	0.489	0.481	0.637	0.639
	$x_{21},x_{22}$	0.574	0.582	0.539	0.571
	$x_{21},x_{22},x_{23}$	0.676	0.706	0.410	0.485
200	$x_{11}$	0.624	0.618	0.856	0.850
	$x_{11},x_{12}$	0.774	0.774	0.887	0.908
	$x_{11},x_{12},x_{13}$	0.994	0.995	0.937	0.974
	$x_{21}$	0.614	0.608	0.820	0.820
	$x_{21},x_{22}$	0.764	0.766	0.857	0.871
	$x_{21},x_{22},x_{23}$	0.958	0.964	0.865	0.903

Table 9. Table 9: Stage 1 model selection ( n = 100 𝑛 100 n=100 ) when the stage 2 model is correctly specified (‘Correct’), selected by the corresponding approach (‘Recommended’), or an intercept-only model (‘Intercept’), under the setup described in section 3.2 but with normal (rather than log-normal) errors ϵ ∼ N ​ ( 0 , 1 ) similar-to italic-ϵ 𝑁 0 1 \epsilon\sim N(0,1) . (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

Selection	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
Correct	$x_{11}$	0.556	0.543	0.767	0.763
	$x_{11},x_{12}$	0.643	0.646	0.644	0.705
	$x_{11},x_{12},x_{13}$	0.796	0.842	0.513	0.652
Recommended	$x_{11}$	0.555	0.541	0.750	0.748
	$x_{11},x_{12}$	0.636	0.641	0.627	0.692
	$x_{11},x_{12},x_{13}$	0.784	0.831	0.492	0.629
Intercept	$x_{11}$	0.552	0.539	0.721	0.720
	$x_{11},x_{12}$	0.628	0.633	0.576	0.640
	$x_{11},x_{12},x_{13}$	0.730	0.776	0.419	0.557

Table 10. Table 10: Model selection ( n = 100 𝑛 100 n=100 ) with true blip parameters of 1, 0.5, or 0.1, under the setup described in section 3.2 but with normal (rather than log-normal) errors ϵ ∼ N ​ ( 0 , 1 ) similar-to italic-ϵ 𝑁 0 1 \epsilon\sim N(0,1) . (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

${\psi }_{jp}$	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
1	$x_{11}$	0.556	0.543	0.767	0.763
	$x_{11},x_{12}$	0.643	0.646	0.644	0.705
	$x_{11},x_{12},x_{13}$	0.796	0.842	0.513	0.652
	$x_{21}$	0.489	0.481	0.637	0.639
	$x_{21},x_{22}$	0.574	0.582	0.539	0.571
	$x_{21},x_{22},x_{23}$	0.676	0.706	0.410	0.485
0.5	$x_{11}$	0.387	0.379	0.378	0.377
	$x_{11},x_{12}$	0.308	0.318	0.164	0.201
	$x_{11},x_{12},x_{13}$	0.265	0.322	0.071	0.127
	$x_{21}$	0.297	0.293	0.310	0.305
	$x_{21},x_{22}$	0.273	0.280	0.142	0.158
	$x_{21},x_{22},x_{23}$	0.217	0.247	0.060	0.088
0.1	$x_{11}$	0.130	0.127	0.075	0.076
	$x_{11},x_{12}$	0.039	0.044	0.006	0.012
	$x_{11},x_{12},x_{13}$	0.014	0.028	0.002	0.006
	$x_{21}$	0.132	0.131	0.104	0.105
	$x_{21},x_{22}$	0.058	0.064	0.022	0.024
	$x_{21},x_{22},x_{23}$	0.027	0.032	0.000	0.003

Table 11. Table 11: Model selection ( n = 100 𝑛 100 n=100 ) with non-treatment covariates drawn from multivariate normal distributions with correlations of 0 (‘none’), 0.25 (‘medium’) and 0.5 (‘strong’), under the setup described in section 3.2 but with normal (rather than log-normal) errors ϵ ∼ N ​ ( 0 , 1 ) similar-to italic-ϵ 𝑁 0 1 \epsilon\sim N(0,1) . (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

Correlation	Model	QICG (F)	QICG (B)	Wald (F)	Wald (B)
None	$x_{11}$	0.556	0.543	0.767	0.763
	$x_{11},x_{12}$	0.643	0.646	0.644	0.705
	$x_{11},x_{12},x_{13}$	0.796	0.842	0.513	0.652
	$x_{21}$	0.489	0.481	0.637	0.639
	$x_{21},x_{22}$	0.574	0.582	0.539	0.571
	$x_{21},x_{22},x_{23}$	0.676	0.706	0.410	0.485
Medium	$x_{11}$	0.510	0.502	0.666	0.661
	$x_{11},x_{12}$	0.588	0.594	0.550	0.570
	$x_{11},x_{12},x_{13}$	0.714	0.732	0.451	0.493
	$x_{21}$	0.468	0.460	0.585	0.584
	$x_{21},x_{22}$	0.509	0.517	0.443	0.459
	$x_{21},x_{22},x_{23}$	0.578	0.597	0.324	0.363
Strong	$x_{11}$	0.470	0.466	0.586	0.587
	$x_{11},x_{12}$	0.525	0.526	0.445	0.449
	$x_{11},x_{12},x_{13}$	0.563	0.574	0.282	0.293
	$x_{21}$	0.407	0.402	0.494	0.488
	$x_{21},x_{22}$	0.410	0.413	0.316	0.322
	$x_{21},x_{22},x_{23}$	0.402	0.410	0.190	0.178

Table 12. Table 12: Non-aggregated selection ( n = 100 𝑛 100 n=100 ), under the setup described in section 3.2 but with normal (rather than log-normal) errors ϵ ∼ N ​ ( 0 , 1 ) similar-to italic-ϵ 𝑁 0 1 \epsilon\sim N(0,1) . (F) and (B) denote forwards and backwards selection, respectively. Bold indicates most successful approach for each setup.

Stage 1	Stage 2	QICG (F)	QICG (B)	Wald (F)	Wald (B)
		Stage 1 Selection
$x_{11}$	$x_{21}$	0.551	0.535	0.775	0.773
	$x_{21},x_{22}$	0.561	0.546	0.766	0.764
	$x_{21},x_{22},x_{23}$	0.557	0.547	0.762	0.756
$x_{11},x_{12}$	$x_{21}$	0.645	0.647	0.650	0.712
	$x_{21},x_{22}$	0.637	0.641	0.634	0.699
	$x_{21},x_{22},x_{23}$	0.646	0.649	0.651	0.706
$x_{11},x_{12},x_{13}$	$x_{21}$	0.799	0.839	0.525	0.663
	$x_{21},x_{22}$	0.794	0.842	0.507	0.644
	$x_{21},x_{22},x_{23}$	0.794	0.846	0.511	0.653
		Stage 2 Selection
$x_{11}$	$x_{21}$	0.498	0.490	0.642	0.643
$x_{11},x_{12}$		0.491	0.481	0.646	0.643
$x_{11},x_{12},x_{13}$		0.477	0.473	0.623	0.631
$x_{11}$	$x_{21},x_{22}$	0.591	0.600	0.558	0.589
$x_{11},x_{12}$		0.582	0.589	0.541	0.576
$x_{11},x_{12},x_{13}$		0.548	0.556	0.519	0.548
$x_{11}$	$x_{21},x_{22},x_{23}$	0.702	0.729	0.443	0.512
$x_{11},x_{12}$		0.667	0.700	0.408	0.484
$x_{11},x_{12},x_{13}$		0.658	0.688	0.380	0.460

Equations

$${\gamma }_j({\mathbf{h}}_j,a_j)=E[Y({\overline{\mathbf{a}}}_{j-1},a_j,{\underline{\mathbf{a}}}_{j+1}^{opt})-Y({\overline{\mathbf{a}}}_{j-1},a_j^0,{\underline{\mathbf{a}}}_{j+1}^{opt})|{\mathbf{h}}_j]$$

$$E[Y({\underline{\mathbf{a}}}_J)]=E[Y^{opt}]-\sum _{j=1}^J\left[{\gamma }_j({\mathbf{h}}_j,a_j^{opt})-{\gamma }_j({\mathbf{h}}_j,a_j)\right]$$

$$G_j({\underline{\mathbf{\psi }}}_j)$$

$$U_j({\underline{\mathbf{\psi }}}_j;{\mathbf{\beta }}_j;{\mathbf{\alpha }}_j)=\left\{S_j(A_j)-E[S_j(A_j)|{\mathbf{h}}_j;{\mathbf{\alpha }}_j]\right\}\{G_j({\underline{\mathbf{\psi }}}_j)-E[G_j({\underline{\mathbf{\psi }}}_j)|{\mathbf{h}}_j;{\mathbf{\beta }}_j]\}$$

$$E[Y(a_1,a_2)]=E[Y(a_1^{opt},a_2^{opt})]\frac{E[Y(a_1,a_2^{opt})|{\mathbf{h}}_1]}{E[Y(a_1^{opt},a_2^{opt})|{\mathbf{h}}_1]}\frac{E[Y(a_1,a_2)|{\mathbf{h}}_2]}{E[Y(a_1,a_2^{opt})|{\mathbf{h}}_2]}$$

$$E[Y({\underline{\mathbf{a}}}_J)|{\underline{\mathbf{h}}}_J]=E[Y^{opt}]\prod _{j=1}^J\left[{\gamma }_j({\mathbf{h}}_{\psi j},a_j)/{\gamma }_j({\mathbf{h}}_{\psi j},a_j^{opt})\right]$$

$$\log (E[Y({\underline{\mathbf{a}}}_J)|{\underline{\mathbf{h}}}_J])=\log (E[Y^{opt}])-\sum _{j=1}^J\left[\log \left({\gamma }_j({\mathbf{h}}_{\psi j},a_j^{opt})/{\gamma }_j({\mathbf{h}}_{\psi j},a_j)\right)\right],$$

$${\tilde{y}}_j=y\times \prod _{k=j+1}^J\left[\left({\gamma }_k({\mathbf{h}}_{\psi k},a_k^{opt};{\mathbf{\psi }}_k)/{\gamma }_k({\mathbf{h}}_{\psi k},a_k;{\mathbf{\psi }}_k)\right)\right],$$

$$\log \left\{E[G_j({\underline{\mathbf{\psi }}}_j)|{\mathbf{h}}_{\beta j};{\mathbf{\beta }}_j]\right\}={\mathbf{h}}_{\beta j}{\mathbf{\beta }}_j\log {\gamma }_j({\mathbf{h}}_{\psi j},a_j;{\mathbf{\psi }}_j)=a_j{\mathbf{h}}_{\psi j}{\mathbf{\psi }}_j,$$

$$U_j({\underline{\mathbf{\psi }}}_j;{\mathbf{\beta }}_j;{\mathbf{\alpha }}_j)$$

$$0=\sum _{i=1}^n{d}_i{\mathbf{h}}_{\psi i}({\tilde{y}}_i-{\mu }_i(\mathbf{\beta },\mathbf{\psi })),$$

$${\tilde{y}}_{J-1}=0.001\times \left({\gamma }_J({\mathbf{h}}_{\psi J},{\hat{a}}_J^{opt};{\hat{\mathbf{\psi }}}_J)/{\gamma }_J({\mathbf{h}}_{\psi J},a_J;{\hat{\mathbf{\psi }}}_J)\right)$$

$${\hat{\mathbf{\psi }}}^{(1)}$$

$${\hat{\mathbf{\beta }}}^{(1)}$$

$${\hat{\mathbf{\beta }}}_j({\mathbf{\psi }}_j,{\underline{\hat{\mathbf{\psi }}}}_{j+1})={\left[{\mathbf{h}}_{\beta j}^{\top }{\mathbf{h}}_{\beta j}\right]}^{-1}{\mathbf{h}}_{\beta j}^{\top }({\tilde{\mathbf{y}}}_j-{\mathbf{A}}_j{\mathbf{h}}_{\psi j}{\mathbf{\psi }}_j).$$

$$U(\mathbf{\psi })=(\mathbf{D}{\mathbf{h}}_{\psi }{)}^{\top }(\tilde{\mathbf{y}}-{\mathbf{h}}_{\beta }\hat{\mathbf{\beta }}(\mathbf{\psi })-\mathbf{A}{\mathbf{h}}_{\psi }\mathbf{\psi })$$

$$\hat{\mathbf{\psi }}={\left({\mathbf{h}}_{\psi }^{\top }\mathbf{W}\mathbf{A}{\mathbf{h}}_{\psi }\right)}^{-1}{\mathbf{h}}_{\psi }^{\top }\mathbf{W}\tilde{\mathbf{y}}.$$

$$\frac{\partial Q}{\partial \mathbf{\psi }}=\frac{\partial Q}{\partial \mathbf{\mu }}\frac{\partial \mathbf{\mu }}{\partial \mathbf{\psi }}={\mathbf{h}}_{\psi }^{\top }\mathbf{W}(\tilde{\mathbf{y}}-\mathbf{\mu }),$$

$$Q(\mathbf{\psi })$$

$$0=\sum _{i=1}^n{d}_i{\mathbf{h}}_{\psi i}(y_i^{\ast }-{\mu }_i^{\ast }(\mathbf{\psi }))$$

$$\frac{\partial Q}{\partial \mathbf{\psi }}=\frac{\partial Q}{\partial {\mathbf{\mu }}^{\ast }}\frac{\partial {\mathbf{\mu }}^{\ast }}{\partial \mathbf{\psi }}={\mathbf{h}}_{\psi }^{\top }\mathbf{D}({\tilde{\mathbf{y}}}^{\ast }-{\mathbf{\mu }}^{\ast }),$$

$$\sqrt{n}({\hat{\mathbf{\psi }}}_{(m)}-{\mathbf{\psi }}_{(m,\ast )})\stackrel{d}{⟶}\text{Normal}(0,\mathcal{V}({\mathbf{\psi }}_{(m,\ast )}))$$

$$\mathcal{I}({\mathbf{\psi }}^{\prime })={E\left[-\frac{{\partial }^2{Q}_1(\mathbf{\psi })}{\partial \mathbf{\psi }\partial {\mathbf{\psi }}^{\top }}\right]|}_{\mathbf{\psi }={\mathbf{\psi }}^{\prime }}\mathcal{J}({\mathbf{\psi }}^{\prime })={E\left[\left\{\frac{\partial Q_1(\mathbf{\psi })}{\partial \mathbf{\psi }}\right\}{\left\{\frac{\partial Q_1(\mathbf{\psi })}{\partial \mathbf{\psi }}\right\}}^{\top }\right]|}_{\mathbf{\psi }={\mathbf{\psi }}^{\prime }}.$$

$$\frac{\partial Q_1(\mathbf{\psi })}{\partial \mathbf{\psi }}=D_1{A}_1{\mathbf{h}}_{\psi 1}(Y_1^{\ast }-{\mu }_1^{\ast }(\mathbf{\psi }))$$

$$\Delta (m)$$

$${\text{QIC}}_G(m)=\hat{\Delta }(m)=-2Q({\hat{\mathbf{\psi }}}_{(m)})+2\text{tr}\{\mathbf{J}({\hat{\mathbf{\psi }}}_{(m)})\mathbf{I}({\hat{\mathbf{\psi }}}_{(m)}{)}^{-1}\}$$

$${\text{QIC}}_{\text{G}}=-2Q(\hat{\mathbf{\psi }})+2\text{tr}\{\mathbf{I}(\hat{\mathbf{\psi }})\hat{\mathbf{V}}(\hat{\mathbf{\psi }}))\}.$$

$$U(\mathbf{\psi })$$

$$\gamma ({\mathbf{h}}_{\psi },a;{\mathbf{\psi }}_1,{\mathbf{\psi }}_2)=a{\mathbf{h}}_{{\psi }_1}{\mathbf{\psi }}_1+a^2{\mathbf{h}}_{{\psi }_2}{\mathbf{\psi }}_2$$

$$U(\psi )$$

$$Q(\mathbf{\psi })$$

Full text

Keywords: Adaptive treatment strategies; Dynamic treatment regimes; Iteratively-reweighted least squares;   Quasi-likelihood Information Criterion; Structural nested models.

1 Introduction

Dynamic treatment regimes (DTRs) - sequences of decision rules that take patient information as input and output recommended treatments - are part of a rapidly expanding literature on personalized medicine (Chakraborty and Moodie 2013). By tailoring treatments to individual patient characteristics, DTRs are able to improve long-term outcomes for a population when compared with more traditional non-tailored approaches. Identification of the optimal regime (which maximizes expected outcome) is a major challenge due to, for example, delayed treatment effects and covariate-dependent treatment assignment.

Numerous methods have been proposed for optimal DTR estimation. A general class of DTR estimation approaches relies on structural nested mean models (SNMMs, Robins 1994; Vansteelandt and Goetghebeur 2003). In our formulation, the SNMM parameterizes the difference between the conditional expectation of the outcome following observed treatment with that of a counterfactual outcome under a (potentially unobserved) treatment regime. By estimating the parameters of this model we are then able to identify the optimal DTR, i.e. the sequence of treatment decisions that maximizes the expected outcome across all patients. This general approach of parameterizing and estimating components of the outcome mean model is used in a variety of specific DTR estimation methods, including Q-learning (Watkins 1989; Sutton and Andrew 1998), dynamic weighted least squares (Wallace and Moodie 2015), and G-estimation (Robins 2004), the last of which is the focus of this paper.

Almost all of the methodological developments for DTR estimation have focused on continuous outcomes and additive effects of treatment on the expected counterfactual outcome, with time-to-event outcomes included as a special case. Estimation for discrete outcomes, or for effects of treatment on non-additive scales has received little attention. A recent exception is the work of Moodie et al. (2014) who used generalized additive models to apply Q-learning in this setting. Thall et al. (2000), meanwhile, considered a likelihood-based approach in the case of a binary outcome. Such examples are rare, however, and typically grounded in methods that do not offer a great deal of flexibility or robustness in modeling. The presentation (and implementation) of G-estimation primarily for continuous outcomes therefore represents an important limitation of the approach.

Even for continuous outcomes, there has been little focus on model selection in the context of G-estimation. The methods listed above all implicitly assume that the SNMMs upon which they rely are correctly (or possibly over-) specified. Very little work has been published related to the problem of choosing between a set of candidate models or model checking. Exceptions include the diagnostic plots of Rich et al. (2010) and the method of Wallace et al. (2016) that exploits the so-called double-robustness property (discussed below) for model assessment. Neither of these, however, assess the component of the model quantifying the effect of treatment – i.e., the blip model – alone, and can at best assess the validity of both the blip model and another component model simultaneously.

In this paper, we present two generalizations of G-estimation. First, we derive and illustrate how iteratively-reweighted least squares (IRLS) may be used to implement G-estimation in a discrete-outcome scenario using log-linear models. We then present a new approach to model selection when using G-estimation for DTRs based on a Quasi-likelihood Information Criterion (QIC). Our QIC formulation is applicable to G-estimation procedures in general, but we demonstrate how the QIC can be applied to G-estimation in the DTR setting, encompassing multiple stages of treatment for both the cases of continuous and count outcomes.

2 DTRs and G-estimation

We establish notation by considering G-estimation in its conventional form, where effects of exposure are additive and a linear model is presumed for the counterfactual outcomes. We consider a cohort of subjects on whom data are gathered at fixed intervals (such as visits to a physician) or at fixed clinical decision points (diagnosis, remission, and so on), with a treatment decision made at each of these time points. Our objective is to identify the sequence of treatment decision rules (the DTR) which maximizes a subject’s long-term expected outcome (defined such that larger values are preferred). We assume there are a total of $J$ successive treatment decisions (or stages): $y$ denotes observed patient outcome; $a_{j}$ denotes the stage $j$ treatment decision ($j=1,...,J$), with $a_{j}^{0}$ denoting “no treatment” (such as a control or standard care); $\bm{h}_{j}$ denotes the covariate matrix containing patient information (history) prior to the $j^{th}$ treatment decision. The history can include previous treatments $a_{1},...,a_{j-1}$ along with non-treatment information $\bm{x}_{j}$. In addition, over- and underline notation is used to indicate the past and future, respectively. For example $\overline{\bm{a}}_{j}$ denotes the vector of treatment decisions up to and including the stage $j$ decision, while $\underline{\bm{a}}_{j+1}$ denotes the last $J-j$ decisions (from stage $j+1$ up to and including stage $J$). The optimal treatment at any given stage is denoted ${a}_{j}^{opt}$.

The (stage $j$) optimal blip-to-reference (or simply blip) function is defined as

[TABLE]

which is the expected difference in outcome when using a reference treatment ${a}_{j}^{0}$ instead of ${a}_{j}$ at stage $j$, in subjects with history $\bm{h}_{j}$ who receive optimal treatment across the remaining $J-j$ intervals ($\underline{\bm{a}}_{j+1}^{opt}$). The optimal treatment at stage $j$ maximizes the blip. Under additive local rank preservation (see 2.1.3 of Chakraborty and Moodie 2013), we can decompose the expectation of the observed potential outcome as

[TABLE]

where ${Y}^{opt}$ can be thought of as the optimal outcome that would be observed if the optimal treatment was followed at every stage. The observed outcome ${y}$ is then equal in expectation to the optimal outcome minus the difference in outcome between optimal and observed treatment at each stage.

In practice, we assume $\gamma_{j}(.)$ takes a known parametric form ${\gamma}_{j}(\bm{h}_{j},{a}_{j};\bm{\psi}_{j})$ with parameters $\bm{\psi}_{j}$. We then estimate $\bm{\psi}_{j}$, and identify the optimal treatment regime by choosing, for each subject, the treatment that maximizes the estimated blip. G-estimation is one method which may be used to estimate $\bm{\psi}_{j}$, and relies on two standard assumptions: the stable unit treatment value assumption and the assumption of no unmeasured confounding (or sequential randomization). The former means that a subject’s outcome is not influenced by other subjects’ treatment allocation (Rubin 1980) and that the counterfactual outcome under a particular treatment is equal to the observed outcome under that treatment; the latter states that the treatment received at stage $j$ is independent of any future (potential) covariate or outcome, conditional on history $\bm{h}_{j}$.

Writing $\bm{\underline{\psi}}_{j}=(\bm{\psi}_{j},\bm{\psi}_{j+1},...,\bm{\psi}_{J})$, we define for each $j$, ${G}_{j}(\bm{\underline{\psi}}_{j})=\widetilde{{y}}_{j}-{\gamma}_{j}(\bm{h}_{\psi j},{a}_{j};\bm{\psi}_{j})$ where $\widetilde{{y}}_{j}={y}+\sum_{k=j+1}^{J}\left[{\gamma}_{k}(\bm{h}_{\psi k},{a}_{k}^{opt};\bm{\psi}_{k})-{\gamma}_{k}(\bm{h}_{\psi k},{a}_{k};\bm{\psi}_{k})\right]$ can be viewed as a pseudo-outcome which we compute at each stage based on those $\widehat{\bm{\psi}}_{k}$ ($k>j$), and hence $\widehat{a}_{k}^{opt}$ already estimated. Therefore

[TABLE]

and we can regard ${G}_{j}(\bm{\underline{\psi}}_{j})$ as being equal to the expected outcome with the effects of stage $j$ treatment ‘removed’ and the difference between optimal and observed treatment thereafter ‘added’. Under the above assumptions we have that $E[{G}_{j}(\bm{\underline{\psi}}_{j})|\bm{h}_{j}]=E[Y(\overline{\bm{a}}_{j-1},{a}_{j}^{0},\underline{\bm{a}}_{j+1}^{opt})|\bm{h}_{j}]$ which represents the expected outcome for a subject who receives treatment history $\overline{\bm{a}}_{j-1}$ up to stage $j-1$, no treatment at stage $j$, and optimal treatment thereafter. We refer to ${G}_{j}$ as the stage $j$ treatment-free outcome.

To estimate the blip parameters $\bm{\psi}_{j}$, G-estimation considers the set of functions

[TABLE]

where typically ${S}_{j}({A}_{j})={a}_{j}\bm{h}_{j}$. A fully efficient form of ${S}_{j}({A}_{j})$ has been proposed (Robins 2004), but requires knowledge of the variance of ${G}_{j}(\bm{\underline{\psi}}_{j})$, which is rarely available in practice. The estimating functions require the specification of a number of models, namely the stage $j$ blip model: ${\gamma}_{j}(\bm{h}_{\psi j},{a}_{j};\bm{\psi}_{j})$; the stage $j$ treatment-free model: $E[{G}_{j}(\bm{\underline{\psi}}_{j})|\bm{h}_{\beta j};\bm{\beta}_{j}]$; and the stage $j$ treatment model: $E[{A}_{j}|\bm{h}_{\alpha j};\bm{\alpha}_{j}]$ (or, more generally, $E[{S}_{j}({A}_{j})|\bm{h}_{j};\bm{\alpha}_{j}]$), where $\bm{h}_{\psi j}$, $\bm{h}_{\beta j}$ and $\bm{h}_{\alpha j}$ are subsets of patient history that feature in the blip, treatment-free, and treatment models, respectively. An important property of G-estimation is its double-robustness: if the blip is correctly specified, then as long as at least one of the treatment and treatment-free models is also correctly specified the resulting blip parameter estimators will be consistent.

Because the functions ${G}_{j}(\bm{\underline{\psi}}_{j})$ depend on the observed outcome ${y}$ and each blip model from stage $j$ onwards, G-estimation proceeds recursively, starting at the final stage $J$ and working backwards to stage 1. At each stage the above models are specified, and then the following three steps are carried out:

Estimate the treatment model parameters $\widehat{\bm{\alpha}}_{j}$ by regressing the stage $j$ treatment ${a}_{j}$ on the treatment model covariates $\bm{h}_{\alpha j}$. 2. 2.

Estimate the treatment-free model parameters $\bm{\beta}_{j}$ by ‘regressing’ ${G}_{j}(\bm{\underline{\psi}}_{j})$ on $\bm{h}_{\beta j}$, where rather than conducting a standard least squares regression, we instead solve the corresponding least squares equation to give $\widehat{\bm{\beta}}_{j}(\bm{\psi}_{j},\bm{\underline{\widehat{\psi}}}_{j+1})$ in terms of the stage $j$ blip parameters $\bm{\psi}_{j}$ and the estimated blip parameters $(\bm{\underline{\widehat{\psi}}}_{j+1})$ from previous stages. 3. 3.

Using the estimates $\widehat{\bm{\alpha}}_{j}$ and $\widehat{\bm{\beta}}_{j}(\bm{\psi}_{j},\bm{\underline{\widehat{\psi}}}_{j+1})$ from steps 1 and 2, solve the equation $E_{n}[{U}_{j}(\bm{\psi}_{j};\widehat{\bm{\beta}}_{j},\widehat{\bm{\alpha}}_{j})]=0$ to estimate $\bm{\psi}_{j}$, where $E_{n}$ denotes the mean over all subjects.

We can then use the resulting blip parameter estimates $\widehat{\bm{\psi}}_{j}$ to estimate the optimal stage $j$ treatment ${a}_{j}^{opt}$ for each subject, and hence the function ${G}_{j-1}(\bm{\underline{\psi}}_{j-1})$, and repeat the above steps until estimates are obtained for every stage of the analysis.

3 G-estimation for generalized linear models

The framework in section 2 is standard for G-estimation applications. It assumes a continuous outcome, and that the treatment modifies the expected outcomes additively, that is, the blip acts additively on the original outcome scale. The construction can be modified to be applicable to discrete outcomes, but relaxing an assumption of additivity of the treatment effect needs more care. In this section, we demonstrate how G-estimation may be generalized to handle other effect types, and show how estimation can be achieved using standard computational approaches. Specifically, we will apply G-estimation for generalized linear models by using iteratively-reweighted least squares.

3.1 G-estimation for multiplicative effects

For an arbitrary counterfactual outcome $Y(a)$, the effect of exposure may be framed in terms of the average potential outcome $E[Y(a)]$, and contrasts comparing this average for different exposures. For example, for a binary exposure we might consider the ratio of expectations $E[Y(1)]/E[Y(0)]$ rather than the expected ratio $E[Y(1)/Y(0)]$; this focuses on population- rather than individual-level contrasts and avoids identifiability issues associated with attempting to specify a joint model for $\{Y(0),Y(1)\}$.

For G-estimation, consider for illustration the two interval case; our approach will focus on constructing models for $E[Y(a_{1},a_{2})]$ using the decomposition

[TABLE]

that is, using a multiplicative modification of the optimal outcome, and making a multiplicative rank preserving assumption. In this context, the blip function ${\gamma}_{j}(\bm{h}_{\psi j},{a}_{j})$ may be defined as the ratio of expected counterfactual outcomes $E[{Y}(\overline{\bm{a}}_{j-1},a_{j},\underline{\bm{a}}_{j+1}^{opt})|\bm{h}_{\psi j}]/E[{Y}(\overline{\bm{a}}_{j-1},{a}_{j}^{0},\underline{\bm{a}}_{j+1}^{opt})|\bm{h}_{\psi j}]$ and the expected counterfactual outcome may be computed as

[TABLE]

or equivalently

[TABLE]

giving rise to a stage-$j$ pseudo-outcome analogous to that in the continuous outcome, linear model setting as

[TABLE]

and hence ${G}_{j}(\bm{\underline{\psi}}_{j})=\widetilde{{y}}_{j}/{\gamma}_{j}(\bm{h}_{\psi j},{a}_{j};\bm{\psi}_{j})$. We then propose log-linear models for the treatment-free and blip models

[TABLE]

from which, via some rearrangement, the G-estimating functions become

[TABLE]

Again suppressing stage-specific notation, and introducing subscript-$i$ notation for subject $i$, G-estimation at each stage thus solves

[TABLE]

with $d_{i}=a_{i}-E[{A}_{i}|\bm{h}_{\alpha i};\bm{\alpha}]$ and $\mu_{i}(\bm{\beta},\bm{\psi})=\exp(\bm{h}_{\beta i}\bm{\beta}+{a}_{i}\bm{h}_{\psi i}\bm{\psi})$. In section 3.2, we present an IRLS algorithm to estimate blip parameters at each stage in the usual recursive manner.

Note that when the observed $y$ value is zero, the pseudo-outcomes in (1) will also be zero unless a further adjustment is made. A simple approach to this issue is to assume that when $y=0$, it is drawn from a Poisson distribution with mean 0.001 (or some other small value), and replace $y$ with its expectation. For example, for stage $J-1$, after parameters $\bm{\psi}_{J}$ are estimated, the usual adjustment $\widetilde{{y}}_{J-1}=y\times\left({\gamma}_{J}(\bm{h}_{\psi J},\widehat{a}_{J}^{opt};\widehat{\bm{\psi}}_{J})/{\gamma}_{J}(\bm{h}_{\psi J},{a}_{J};\widehat{\bm{\psi}}_{J})\right)$ becomes

[TABLE]

if $y$ is zero; $\widetilde{{y}}_{J-1}$ is guaranteed non-negative.

3.2 Iteratively-reweighted least squares

We now demonstrate how G-estimation may proceed for log-linear models using IRLS. Without loss of generality we consider a single-stage example allowing us to suppress stage-specific notation. The G-estimation equations, as written in (2) are of a standard form from which IRLS may be used to estimate $\bm{\psi}$. Suppose ${y}$ has a mean function $\mu$ modeled using link function $g(\cdot)$ and linear predictor vector ${\eta}=\bm{h}_{\beta}\bm{\beta}+{a}\bm{h}_{\psi}\bm{\psi}$ such that ${\mu}=g^{-1}({\eta})$. Denote the variance function $V(\mu)$. We can then estimate $\bm{\psi}$ via IRLS using the following algorithm:

1. Set initial parameters $\widehat{\bm{\beta}}^{(0)}$, $\widehat{\bm{\psi}}^{(0)}$ and compute for each subject the initial linear predictor $\widehat{{\eta}}^{(0)}=\bm{h}_{\beta}\widehat{\bm{\beta}}^{(0)}+a\bm{h}_{\psi}\widehat{\bm{\psi}}^{(0)}$ and mean value $\widehat{{\mu}}^{(0)}=g^{-1}(\widehat{{\eta}}^{(0)})$.

2. Set $\widehat{z}_{i}^{(1)}=\widehat{\eta}_{i}^{(0)}+(y_{i}-\widehat{\mu}_{i}^{(0)})\dot{g}(\widehat{\mu}_{i}^{(0)})$, $\widehat{w}_{i}^{(1)}=w_{i}/[\{\dot{g}(\widehat{\mu}_{i}^{(0)})\}^{2}V(\widehat{\mu}_{i}^{(0)})]$. Denote by $\bm{D}^{(1)}$ the diagonal matrix with $(i,i)$ element $\widehat{d}_{i}^{(1)}$, and by $\bm{A}$ the diagonal matrix with $(i,i)$ element $a_{i}$, the observed treatment for subject $i$.

3. Apply the G-estimation procedure to re-estimate $\bm{\beta}$ and $\bm{\psi}$:

[TABLE]

where $\bm{h}_{\beta D}=\bm{h}_{\beta}(\bm{h}_{\beta}^{\top}\bm{D}^{(1)}\bm{h}_{\beta})^{-1}\bm{h}_{\beta}^{\top}\bm{D}^{(1)}$, and $\mathbf{I}_{n}$ is the $n\times n$ identity matrix.

4. Define vectors $\widehat{\bm{\eta}}^{(1)}=\bm{h}_{\beta}\widehat{\bm{\beta}}^{(1)}+\bm{A}\bm{h}_{\psi}\widehat{\bm{\psi}}^{(1)}$ and $\widehat{\bm{\mu}}^{(1)}=g^{-1}\left(\widehat{\bm{\eta}}^{(1)}\right)$.

5. Return to 2 and iterate through 2-5 using $\widehat{\bm{\mu}}^{(1)}$ and $\widehat{\bm{\eta}}^{(1)}$ as the updated starting values, obtaining ($\widehat{\bm{\mu}}^{(2)}$,$\widehat{\bm{\eta}}^{(2)}$), then repeat to generate ($\widehat{\bm{\mu}}^{(3)}$,$\widehat{\bm{\eta}}^{(3)}$), and so on.

6. Repeat until $\widehat{\bm{\mu}}^{(t)}$ and $\widehat{\bm{\eta}}^{(t)}$ satisfy $\left|\widehat{\bm{\mu}}^{(t)}-\widehat{\bm{\mu}}^{(t-1)}\right|<\epsilon_{\mu}$ and/or $\left|\widehat{\bm{\eta}}^{(t)}-\widehat{\bm{\eta}}^{(t-1)}\right|<\epsilon_{\eta}$ for tolerances $\epsilon_{\mu}$ and $\epsilon_{\eta}.$

The sequence of estimates produced by this algorithm converges to the solution of the G-estimating equations.

4 G-estimation and quasi-likelihood

Inference for SNMMs using G-estimation, unlike inference for more conventional models such as generalized linear models (GLMs), is not likelihood-based and so established model selection approaches such as Akaike’s Information Criterion (AIC, Akaike 1973) cannot be directly used within the DTR framework. However, we shall reframe the preceding presentations of G-estimation to illustrate how quasi-likelihood theory may be applied.

4.1 Linear case

We assume the treatment-free model is linear in $\bm{h}_{\beta j}$, i.e. that $E[{G}_{j}(\bm{\underline{\psi}}_{j})|\bm{h}_{\beta j};\bm{\beta}_{j}]=\bm{h}_{\beta j}\bm{\beta}_{j}$. Then by ordinary least squares, we may estimate $\bm{\beta}_{j}$ as

[TABLE]

For convenience we again suppress the subscript-$j$ notation – all of what follows may be applied on a stage-by-stage basis - and write $\widehat{\bm{h}}_{\beta}=\bm{h}_{\beta j}\left[\bm{h}_{\beta j}^{\top}\bm{h}_{\beta j}\right]^{-1}\bm{h}_{\beta j}^{\top}$. Substituting the estimate of $\widehat{\bm{\beta}}$ in terms of $\bm{\psi}$, we may rewrite the estimating function vector as

[TABLE]

where $\bm{W}=\bm{D}^{\top}(\mathbf{I}_{n}-\widehat{\bm{h}}_{\beta})$. From here, the estimation of $\bm{\psi}$ follows by

[TABLE]

The form of this estimator is straightforward (and is almost identical to a standard weighted ordinary least squares estimator). This affords greater simplicity in implementation, as well as giving a clear indication that quasi-likelihood methods may be easily applied.

We follow Wedderburn (1974) in defining the quasi-likelihood of $\bm{\psi}$ by writing $\bm{\mu}=\bm{A}\bm{h}_{\psi}\bm{\psi}$ and solving

[TABLE]

yielding

[TABLE]

where $\bm{m}=\bm{h}_{\psi}^{\top}\bm{W}\widetilde{\bm{y}}$, $\bm{M}=\bm{h}_{\psi}^{\top}\bm{W}\bm{A}\bm{h}_{\psi}=\bm{h}_{\psi}^{\top}\bm{D}(\mathbf{I}_{n}-\widehat{\bm{h}}_{\beta})\bm{A}\bm{h}_{\psi}$, and we ignore the constant term. Because $\mathbf{I}_{n}-\widehat{\bm{h}}_{\beta}$ is positive definite, $\bm{M}$ is positive semi-definite in expectation, and thus provided $n$ is large, this quasi-likelihood is uniquely maximized at $\widehat{\bm{\psi}}$ in large samples. Furthermore, given the stage-by-stage, recursive nature of the G-estimation approach within the DTR setting, we may derive this quasi-likelihood at each stage of an analysis.

4.2 Log-linear case

In the log-linear case we first reformulate (2), dividing through by $\exp(\bm{h}_{\beta}\bm{\beta})$ to give

[TABLE]

where ${y}_{i}^{*}={y}_{i}\exp(-\bm{h}_{\beta i}\bm{\beta})$ and $\mu_{i}^{*}(\bm{\psi})=\exp(a_{i}\bm{h}_{\psi i}\bm{\psi})$. This moves the nuisance parameters $\bm{\beta}$ into a pseudo-outcome $y^{*}$, framing the estimating equations more explicitly in terms of the target blip parameters $\bm{\psi}$, as in the linear case. This allows us to return to the theory of Wedderburn and proceed as before by solving

[TABLE]

but this does not yield a quasi-likelihood in a simple way. However, by appealing to the IRLS procedure, and a recursive calculation, we compute a quasi-likelihood suitable for model comparison by considering the sequence of linear approximations to the log-linear estimating equations implied by (4). The IRLS procedure produces a solution to (6) by utilizing a quadratic approximation to the actual quasi-likelihood at the maximizing value; by standard theory the solution is an $o(1)$ approximation to the actual maximizing value of the quasi-likelihood. This strategy appeals to the common approach of defining a quasi-likelihood from estimating equations by considering the dual quadratic minimization problem (see, for example, Green 1984; McCullagh 1991).

4.3 The quasi-likelihood information criterion

We now address selection of the blip model. Based on the preceding derived quasi-likelihoods, we propose an information criterion whose general form builds on standard likelihood theory, where the Kullback-Leibler divergence between a proposed model and the true, data-generating model is minimized.

Under standard regularity conditions on the quasi-likelihood function $Q(.)$, inference proceeds in the usual way for misspecified models. Let $f(y)$ denote the true, data-generating distribution, and let $\gamma(y;\bm{\psi}_{(m)})$ denote a proposed blip model which, combined with treatment and treatment-free models, fully specifies the G-estimating quasi-likelihood, $Q(y;\bm{\psi}_{(m)})$, and the corresponding density $f_{m}(y)\equiv f(y;\bm{\psi}_{(m)})$. The proposed blip model is taken from a class of candidate models, $\mathcal{M}(m)=\{\gamma(y;\bm{\psi}_{(m)})|\bm{\psi}_{(m)}\in\Psi(m)\}$ with fitted models $\gamma(y;\widehat{\bm{\psi}}_{(m)})$. The divergence between $f(y)$ and $f_{m}(y)$ estimated using the observed data and $\widehat{\bm{\psi}}_{(m)}$ is given (up to an additive constant) by $\delta(\widehat{\bm{\psi}}_{(m)})=E[-2Q(Y;\bm{\psi})]|_{\bm{\psi}_{(m)}=\widehat{\bm{\psi}}_{(m)}}$, computed with $\widehat{\bm{\psi}}_{(m)}$ fixed, and the expected divergence is given by $\Delta(m)=E[\delta(\widehat{\bm{\psi}}_{(m)})]$; in the latter expression, the expectation is over the distribution of the estimator $\widehat{\bm{\psi}}_{(m)}$. All expectations are taken with respect to the true distribution $f(y)$ by considering independent copies of the data.

Let $\bm{\psi}_{(m,*)}=\arg\min_{\bm{\psi}_{(m)}\in\Psi(m)}Q(y;\bm{\psi}_{(m)})$. If the true blip function is parametric and contained in $\mathcal{M}(m)$, then $\bm{\psi}_{(m,*)}$ is the “true” parameter; if the set of candidate models does not contain the true blip, then $\bm{\psi}_{(m,*)}$ is the value such that $Q(y;\bm{\psi}_{(m,*)})$ provides the best approximation to $f(y)$ in the sense of minimizing the expected Kullback-Leibler divergence. Under standard regularity conditions on $Q(.)$, we have that $\widehat{\bm{\psi}}_{(m)}$ is consistent for $\bm{\psi}_{(m,*)}$, and

[TABLE]

where $\mathcal{V}$ is a positive definite matrix given by $\mathcal{V}(\bm{\psi})=\mathcal{I}(\bm{\psi})^{-1}\mathcal{J}(\bm{\psi})\mathcal{I}(\bm{\psi})^{-1}$ where, for $\bm{\psi}^{\prime}\in\mathcal{N}$, an open neighborhood of $\bm{\psi}_{(m,*)}$

[TABLE]

and

[TABLE]

is the G-estimating function inspired by (5) for the first data point.

Theorem: Suppose that $Q(\bm{\psi})$ is twice continuously differentiable with bounded expectation of its second derivative in an open neighborhood $\mathcal{N}$ of $\bm{\psi}_{(m,*)}$. Then, under the stable unit treatment value and no unmeasured confounding assumptions (detailed in Subsection 2.1), the expected divergence $\Delta(m)$ can be approximated as

[TABLE]

which is consistently estimated by

[TABLE]

where $\bm{I}(.)$ and $\bm{J}(.)$ are the observed (empirical) versions of $\mathcal{I}$ and $\mathcal{J}$. Thus, the model selection procedure that chooses a model by minimizing $\text{QIC}_{G}(m)$ across $\mathcal{M}(m)$ identifies the model that minimizes $\Delta(m)$ with probability 1 as $n\longrightarrow\infty$.

Proof: see Supplementary Material.

This result gives rise to our quasi-likelihood information criterion, which in terms of the estimator of the asymptotic variance $\mathcal{V}$, $\widehat{\bm{V}}(\widehat{\bm{\psi}})=n\bm{I}(\widehat{\bm{\psi}})^{-1}\bm{J}(\widehat{\bm{\psi}})\bm{I}(\widehat{\bm{\psi}})^{-1}$, may be written

[TABLE]

In their derivation of a related criterion, Taguri et al. (2014) use a direct sandwich estimator for $\mathcal{V}$. However, while this allows a slight simplification of expression (7), we note that this approach fails to accommodate all sources of uncertainty. The estimation of the parameters ${\bm{\alpha}}$ of the treatment model at each interval should be acknowledged, and as we move through stages recursively estimation of all previous parameters should be similarly accommodated; this is achieved through the application of Taylor expansions to the estimating function ${U}(\bm{\psi})$ (Robins 2004; Moodie 2009), although in our experience such corrections make little difference to the resulting variance estimates.

Our derivation of the QIC also differs from that of Taguri et al. (2014) in two substantial ways. First, the estimation of $\widehat{\bm{\beta}}$ is corrected for automatically by its substitution in the estimation of $\widehat{\bm{\psi}}$, that is, using implicit forms $\bm{\beta}(\bm{\psi})$ in the linear model or the IRLS recursion for the log-linear model. That is, we do not estimate the treatment-free model parameters $\bm{\beta}$ in a separate calculation using only the untreated individuals. Secondly, our derivation of the quasi-likelihood matches that of Taguri et al. (2014) in the linear case, however their approach cannot be extended to the log-linear case.

The form of (7) is typical in information criterion-style approaches (Takeuchi 1976), and writing $K=\text{tr}\{\bm{I}(\widehat{\bm{\psi}})\widehat{\bm{V}}(\widehat{\bm{\psi}}))\}$ we may present it as $\text{QIC}_{\text{G}}=-2Q(\widehat{\bm{\psi}})+2K$ to more clearly evoke this similarity. This criterion may be applied at each stage of the G-estimation process with the blip model returning the lowest criterion value being recommended, as in a more typical analysis. Note, however, that it is necessary to assume that at all but the first stage of treatment an at-worst overspecified blip model is contained within the set of candidate models as otherwise poor parameter estimation can have a cumulative effect. Similarly, we must assume that at least one of the treatment or treatment-free models is correctly specified, so that the resulting blip parameter estimators are consistent. These assumptions are necessary for any recursive procedure.

The above theory extends to the case of continuous treatments. The primary complication is that the blip function is extended to include a quadratic treatment term, so that the optimal treatment at any given stage may lie inside the range of possible values it may take (Rich et al. 2014). After this modification, we can proceed to define an equivalent quasi-likelihood (and quasi-likelihood information criterion) at each stage of treatment. Full details are included in the Supplementary Material.

5 Analysis

In this section, we first use simulations to demonstrate the IRLS approach to G-estimation for a count outcome and to demonstrate the performance of the QICG in a continuous outcome scenario. We then proceed to apply both IRLS and the quasi-likelihood information criterion to an empirical analysis, performing analyses which treat the (discrete) outcome as either continuous or as a count.

5.1 Simulation study: IRLS for a log-linear SNMM

First, we present an illustration of the IRLS algorithm for G-estimation in the case of a log-linear outcome model. Simulating a two-stage example, we generate data as follows:

• •

stage 1 patient information: $X_{1}\sim N(0,1)$;

• •

stage 1 treatment: $a_{1}\in\{0,1\}$, ${P}(A_{1}=1|\bm{h}_{1})=\text{expit}(x_{1})$;

• •

stage 2 patient information: $X_{2}\sim N(a_{1},1)$;

• •

stage 2 treatment: $a_{2}\in\{0,1\}$, ${P}(A_{2}=1|\bm{h}_{2})=\text{expit}(x_{2})$;

• •

stage $j$ blip: $\gamma_{j}(a_{j},x_{j})=a_{j}(\psi_{j0}+\psi_{j1}x_{j1})$ such that $a_{j}^{opt}={1}_{\{\psi_{j0}+\psi_{j1}x_{j1}>0\}}$;

• •

outcome: $P(Y=k)=\lambda^{k}e^{-k}/k!$,

with $\lambda=\exp\left[\beta_{0}+\log(|x_{1}|)-\sum_{j=1}^{2}[\gamma_{j}(a_{j}^{opt},x_{j})-\gamma_{j}(a_{j},x_{j})]\right]$.

For all simulations we set $(\psi_{j0},\psi_{j1})=(0.5,-0.5)$, $j=1,2$. As described in section 3.1, one concern in extending G-estimation to discrete outcomes is that of zero values in the response, and the effect they can have on the stage-specific pseudo-outcomes. Replacing the 0 with 0.001 when computing the pseudo-outcome, we investigate the performance of the algorithm in three sets of simulations with varying values for $\beta_{0}$, chosen to yield outcomes with approximately 5%, 10% and 20% zeros. Our analyses correctly specified the treatment model, but mis-specified the treatment-free model, supposing it was linear in $x_{1}$ in contrast to the true $\log(|x_{1}|)$ term. Initial simulation runs had unacceptably high rates of failure to converge. To address this we adjusted the IRLS algorithm detailed above slightly, introducing step-halving whereby the initial parameter estimates at each iteration were the mean of the previous two stages, and ignoring the termination condition dependent on $|\bm{\mu}^{(t)}-\bm{\mu}^{(t-1)}|$, instead terminating when only $|\bm{\eta}^{(t)}-\bm{\eta}^{(t-1)}|$ dropped below a given tolerance. This reduced convergence failure rates to 3% or lower.

We generated 1000 simulated datasets per setup, setting the tolerance $\epsilon_{\eta}$ to 0.001 and limiting the number of iterations at each stage to 1000. Exploratory analyses of smaller simulation runs with lower tolerances and larger iteration limits did not yield substantially different results for parameter estimates or failure rates. Results are summarized (Table 1), where stage 1 estimates for the covariate-by-treatment interaction are slightly (though not statistically significantly) biased in small samples. Bias does not appear to be related to the probability of zero-outcomes, although standard errors appear to increase with it.

5.2 Simulation study: QICG

Next, we demonstrate the use of QICG in the DTR framework with a variety of simulated two-stage examples from the continuous outcome setting (we present results for the discrete-outcome setting in the Supplementary Material). We generate data as follows:

• •

stage 1 patient information: $X_{1k}\sim N(0,1)$ for $k=1,2,3$;

• •

stage 1 treatment: $a_{1}\in\{0,1\}$, ${P}(A_{1}=1|\bm{h}_{1})=\text{expit}(x_{11}+x_{12}+x_{13})$;

• •

stage 2 patient information: $X_{2k}\sim N(a_{1},1)$ for $k=1,2,3$;

• •

stage 2 treatment: $a_{2}\in\{0,1\}$, ${P}(A_{2}=1|\bm{h}_{2})=\text{expit}(x_{21}+x_{22}+x_{23})$;

• •

stage $j$ blip: $\gamma_{j}(a_{j},\bm{h}_{j})=a_{j}(1+\psi_{j1}x_{j1}+\psi_{j2}x_{j2}+\psi_{j3}x_{j3})$

such that $a_{j}^{opt}=1_{\{1+\psi_{j1}x_{j1}+\psi_{j2}x_{j2}+\psi_{j3}x_{j3}>0\}}$;

• •

outcome: $Y=-\sum_{j=1}^{2}[\gamma_{j}(a_{j}^{opt},\bm{h}_{j})-\gamma_{j}(a_{j},\bm{h}_{j})]+\epsilon$, with $\epsilon\sim\text{log-normal}(0,1)-e^{0.5}$;

where $\text{expit}(x)=\left[1+\text{exp}(-x)\right]^{-1}$ is the expit or inverse-logit function. We have used skewed errors in our generation of the outcome (centralized to have mean zero) to better illustrate the potential benefits of the QICG approach. Results using normal errors are included in the Supplementary Material for reference. In our first analyses, we consider datasets of size $n=50,100$ and $200$, and set the blip parameters to $(\psi_{j1},\psi_{j2},\psi_{j3})=(1,0,0)$, $(1,1,0)$ or $(1,1,1)$ giving a range of models including one, two, or all three variables at each stage.

We conducted a G-estimation analysis of 1000 simulated datasets considering eight different blip models corresponding to each of the possible combinations of the predictors at each stage (that is, using none, one, two or all three). The treatment models were always correctly specified (and modeled using logistic regression), while the treatment-free models $E[{G}_{j}(\bm{\underline{\psi}}_{j})|\bm{h}_{\beta j};\bm{\beta}_{j}]$ were linear with covariates $(1,x_{11},x_{12},x_{13})$ at stage 1 and $(1,x_{11},x_{12},x_{13},a_{1}x_{11},a_{1}x_{12},a_{1}x_{13},x_{21},x_{22},x_{23})$ at stage 2 (i.e. using all available covariates).

Using QICG, we performed forward and backward selection within the set of candidate models as in a standard AIC-type stepwise analysis. For comparison we also conducted forward and backward selection based on Wald test p-values at a 0.05 significance level. For these initial simulations, stage 1 results are based on analysis carried out following fitting of the correct model at stage 2. Results (Table 2) indicate QICG outperforms the Wald-type approaches except for the smallest true models (although even then it largely remains competitive). Furthermore, while QICG shows a slight tendency to overfit, the Wald-type approaches show considerably more bias towards underfitting. In addition, we note a slightly greater consistency between the forward and backward QICG results than between the Wald test results, suggesting QICG may be more robust to choice of selection direction.

In these multi-stage simulations, all methods perform better in selection of the first stage blip model rather than the second due to the first stage being a simpler model fitting problem: at stage 2 both first and second stage covariates affect analysis. We note however that the recursive nature of the G-estimation approach means that the analysis of the first stage happens subsequent to analysis of the second stage, using estimates obtained from that stage which we might expect to affect stage 1 model assessment. It seems that any impact this aspect of the estimation process might have is small compared to the inherent simplicity that earlier stages involve fewer covariates.

Moreover, in the results for stage 1 model assessment we have ignored the problem of stage 2 model selection, and instead fixed the stage 2 model as the correct one. While in reality this is an implicit assumption we must make, it seems prudent to investigate what impact this may have on model selection. We conduct analyses identical to those above (limited to $n=100$) but with two other approaches to model selection. In addition to the ‘best-case’ scenario of selecting the correct stage 2 model, we also investigate the consequences of choosing the model recommended by the model selection procedures themselves, and a ‘worst-case’ scenario where intercept-only models are used (Table 3). These non-optimal approaches do result in a drop in performance, but this is not particularly dramatic (and not statistically significant) even in the worst-case scenario. We note, however, that in more complex setups it is likely mis-specification of the stage 2 models could have more dramatic consequences for stage 1 model selection.

We also investigated the impact of other aspects of our data generation setup. Weaker effect sizes (simulated by setting the blip parameters to 0.1 and 0.5 instead of 1 as above) predictably resulted in lower success rates across all the methods under consideration. However, the quasi-likelihood approach was much more resilient to these effects. Introducing a correlation structure among the non-treatment covariates also resulted in worse model selection, but again the quasi-likelihood approach appeared to be slightly more robust to these changes. Results of these additional scenarios are included in the Supplementary Material.

These results were aggregated across the various model setups at each stage to afford greater simplicity in the presentation of our results. For example, the results corresponding to the stage 1 blip model which only included $x_{11}$ were taken from simulations across the three different stage 2 blip models. Non-aggregated results (see Supplementary Material) show little evidence of stage 2 model complexity affecting stage 1 selection or vice-versa.

5.3 Investigating the trace term

We have framed our $\text{QIC}_{G}$ in terms of the quasi-likelihood and a ‘penalty term’ defined in terms of the trace $K=\text{tr}\{\bm{I}(\widehat{\bm{\psi}})\widehat{\bm{V}}(\widehat{\bm{\psi}})\}$. In the likelihood-based setting it is known that this trace term may, if the model under consideration is a good approximation of the truth, be approximated by the dimension of that model (Burnham and Anderson 2002; Shibata 1989). Here, the quasi-likelihood is grounded in estimating equation theory and so, if the error terms in the outcome generating model were normally distributed, we might expect a similar result. This is complicated by our estimation of two other models besides the blip, as well as the recursive multi-stage nature of the G-estimation framework.

We summarize estimates of the term $K$ from our first set of simulations above with $n=100$ in a figure in the supplement, with the modification that we generate the error term $\epsilon$ in the outcome from a standard normal distribution rather than log-normal. In general, the estimates appear to be similar to the dimension of the corresponding model, particularly at the second stage, when the correct model was used. If two models of the same dimension are compared then on average an incorrect model will result in a slightly larger trace term than that from a correct model. Stage 1 estimates of the trace in general seem to be slightly lower. These results assume correct specification of the treatment model, ensuring consistency of our estimators as per the double-robustness property of G-estimation. We have found (results omitted) that when both treatment and treatment-free models are mis-specified the resulting trace terms can be much larger. In addition, if the distribution of the error term $\epsilon$ is skewed, as in the previously reported results, then the trace term is again much larger. This suggests the possibility of comparing the trace term with candidate model dimension to investigate the validity of the treatment and treatment-free models, following use of residual plots (or similar techniques) to assess the normality of the residuals.

We note finally that the use of bootstrap procedures to estimate the trace (or penalty) term $\text{tr}\{\bm{I}(\widehat{\bm{\psi}})\widehat{\bm{V}}(\widehat{\bm{\psi}})\}$ have been recommended in preference to the use of termwise plug-in versions of matrices $\bm{I}(.)$ and $\bm{J}(.)$ (see Burnham and Anderson (2002)). We have not investigated this possibility in our analysis as the plug-in procedure appears to work well in our examples, and due to the additional computational burden. Bootstrap procedures are valid for inference in the regular G-estimation setting (see for example Shortreed and Moodie (2012)), and are straightforward, although computationally expensive, to implement.

5.4 Sequenced Treatment Alternatives to Relieve Depression study

We now illustrate application of our proposed IRLS and QICG approaches to real data from the Sequenced Treatment Alternatives to Relieve Depression (STARD) study. STARD was a multi-stage randomized control trial designed to compare different treatment regimes for patients with major depressive disorder (Fava et al. 2003; Rush et al. 2004). The study was split into 4 levels (one of which was itself split into two sub-levels), with patients receiving a different treatment or combination of treatments within each level.

At study entry (level 1) patients were prescribed citalopram and followed up at regularly scheduled clinic visits. Those whose depression did not enter remission – defined as a Quick Inventory of Depressive Symptomatology (QIDS) score less than or equal to 5 – could proceed to a second level of treatment where seven treatment options were available. The second level of treatment was characterized by ‘switching’ from citalopram to one of four new treatments, or ‘augmenting’ the current treatment by receiving citalopram alongside one of three other treatment options. Patients who received cognitive therapy at level 2 (either alone or combined with citalopram) were eligible to enter the sublevel 2A where they received one of the treatments available at level 2. All patients without remission could then proceed to level 3 (and, if their depression persisted, a further level 4) where again their previous treatment was either switched to or augmented with a number of options. Full details of the study design and treatment options, are described elsewhere (Rush et al. 2004).

An important aspect of the study is that patients were asked for their treatment preference and would then be randomized to one of the treatment options consistent with their preference. This is typically characterized as patients choosing to ‘switch’ from their current treatment to, or ‘augment’ it with, a different one, although the reality was slightly more complex. In particular, when moving from level 1 to level 2, patients were asked about their preference to switch or augment their current treatment with cognitive therapy separately to their preference to switch or augmenting with a pharmacological (i.e., drug-based) treatment.

We conduct an analysis following Chakraborty et al. (2013) who investigated the dichotomy between treatments that were, or included, a selective serotonin reuptake inhibitor (SSRI) and those that did not. We restrict attention to two stages of the study and consider level 2 (including level 2A) and level 3 as our first and second stages of treatment, respectively. Treatment at each stage was coded as 1 if an SSRI was received, either alone or in combination, with level 2A treatments (both of which were non-SSRI) combined with level 2 treatments for this purpose. Treatment was coded as 0 if no SSRI was received throughout a stage. Of 1,027 total patients, only 273 entered level 3. Our outcome is defined as negative QIDS score at end of treatment (i.e., at the end of stage 2 if a patient entered level 3, and at the end of stage 1 otherwise). By taking the negative, larger values are preferred, and we therefore seek a DTR that maximizes this outcome. We pursue an analysis analogous to those undertaken by previous authors, viewing QIDS score as a continuous outcome, and then use our new IRLS-based approach and the associated QICG to apply a log-linear model.

We consider the following tailoring variables: QIDS score measured at the start of the corresponding level (denoted $q_{j}$ for stage $j$), the change in QIDS score divided by time across the previous level (QIDS slope, denoted $s_{j}$), and patient preference prior to receiving treatment ($p_{j}$). Patient preference is binary and coded as 1 if the patient rejected all treatments consistent with switching to a different pharmacological treatment, and 0 otherwise. The assumed treatment model at each stage was fit by logistic regression of observed treatment on preference only, while the treatment-free models are specified as:

• •

stage 1: $E[{G}_{1}(\bm{\underline{\psi}}_{1})|\bm{h}_{\beta 1};\bm{\beta}_{1}]=\beta_{10}+q_{1}\beta_{11}+s_{1}\beta_{12}+p_{1}\beta_{13}$; and

• •

stage 2: $E[{G}_{2}(\bm{\underline{\psi}}_{2})|\bm{h}_{\beta 2};\bm{\beta}_{2}]=\beta_{20}+q_{2}\beta_{21}+s_{2}\beta_{22}+p_{2}\beta_{23}+a_{1}\beta_{24}$.

We consider the ‘full’ blip models

• •

stage 1: $\bm{\gamma}_{1}(\bm{h}_{\psi 1},{a}_{1};\bm{\psi}_{1})=a_{1}(\psi_{10}+q_{1}\psi_{11}+s_{1}\psi_{12}+p_{1}\psi_{13})$; and

• •

stage 2: $\bm{\gamma}_{2}(\bm{h}_{\psi 2},{a}_{2};\bm{\psi}_{2})=a_{2}(\psi_{20}+q_{2}\psi_{21}+s_{2}\psi_{22})$,

and investigate all sub-models (every covariate combination from the full models, including intercept-only models). There are eight candidate models at stage 1, and four at stage 2.

Rather than proceed directly to a stepwise procedure, we instead fit all possible models as dimensionality was low. G-estimation was therefore first performed four times: once for each of the candidate stage 2 blip models, with an intercept-only blip model specified at stage 1 (this choice not affecting stage 2 analysis). The lowest value of $\text{QIC}_{G}$ was found when the intercept-only model was used, suggesting it is the best choice of stage 2 blip model (based on either a forward or backward selection procedure). This result was reinforced by application of both Wald-type approaches of the previous section.

We then repeated the analysis for each of the eight candidate stage 1 blip models using the recommended stage 2 blip model, and found that the model containing stage 1 preference only returned the lowest QICG (and would be recommended by either a forward or backward procedure). The Wald-type approaches also both recommended the preference-only model.

Our log-linear analyses were based on a similar setup with treatment-free and blip models exponentiated, adding 27 to all outcome measures to ensure they were positive. Applying G-estimation via IRLS and computing QICG for every blip model we found the same models recommended at each stage. We observe that in all analyses the same model was recommended at stage 1 regardless of which stage 2 model was chosen.

Overall, these results are broadly consistent with the analysis of Chakraborty et al. who, despite using the somewhat different approach of Q-learning, found that no stage 2 blip covariates were statistically significant, while stage 1 preference alone was significant in the stage 1 blip model. It is encouraging that QICG for both modeling setups indicated the preference term should be included. We also find that at both stages both estimation processes predict the same optimal treatment for every patient.

Finally, we can compare observed outcomes among patients based on how consistent their observed treatments were with the optimal ones as recommended by our models. Among patients who entered stage 2, mean improvement in QIDS score among those who received optimal treatment at both stages, one stage, or no stages was, respectively, 4.67 (sd = 3.92, $n$ = 15), 3.80 (sd = 6.08, $n$ = 75), and 2.56 (sd = 5.19, $n$ = 183). Among patients who did not enter stage 2, mean improvement was 6.13 (sd = 4.72, $n$ = 145) for those who received optimal treatment, and 5.49 (sd = 4.84, $n$ = 609) for those who did not.

6 Discussion

Personalized medicine and the development of dynamic treatment regimes is an important frontier in biostatistical research. This has been reflected in a rapidly expanding literature focusing on DTR estimation techniques, but more practical concerns have received comparatively little attention. In this paper we have presented two extensions to the G-estimation framework. First, we have demonstrated how G-estimation may be applied for log-linear models via iteratively-reweighted least squares; this provides a relatively straightforward route to G-estimation use in a greater variety of contexts. Further, we have presented an approach to model selection for SNMMs within the DTR framework. By demonstrating how G-estimation in its typical application (for continuous treatments) may be reduced to a relatively simple form, we have derived a quasi-likelihood for each stage of a multi-stage, recursive analysis process. We have then extended the work of Pan (2001) and Taguri et al. (2014) to derive a general quasi-likelihood information criterion for DTR estimation using G-estimation. Furthermore, while we have focused on the binary treatment setting, the theory extends to the case of a continuous treatment, dramatically increasing its applicability.

Our simulations involving log-linear SNMMs indicate how G-estimation may be implemented for discrete outcomes with relative ease, especially if extant IRLS routines in standard software packages can be used for this purpose. We note that there may be further room to improve on the approximation used to handle the zero-outcomes, which at present are set to a very small number only when the non-optimal treatment was received. This concern extends to the binary outcome setting, where estimation is even more problematic: the ‘blip’ cannot be separated from the treatment-free component of the mean using a logit transformation, and the use of a log-linear model only provides reasonable (unbiased) estimators in a small range of settings. The question of how to construct pseudo-outcomes that better address this issue is an avenue for further research. Note, finally, that if the (counterfactual) outcomes are binary, then modeling the expected counterfactual on the usual logistic scale is problematic for G-estimation, as it is not possible to separate the blip component from the treatment-free expected counterfactual component, as would be the case on the linear or log-linear scale. This issue warrants further attention.

Through simulation studies we have shown that our quasi-likelihood information criterion performs as well as or better than simpler Wald-type approaches for continuous outcomes, particularly when sample or effect sizes are small, or there is correlation between candidate covariates. In addition, we found greater agreement between the forward and backward stepwise approaches when using QICG than the Wald-type approach, a potentially attractive feature in practice. We note, however, that QICG does seem to overfit, and as such slight modifications may lead to more balanced results. We have experimented with ad hoc corrections inspired by the Bayesian Information Criterion and the corrected AIC, which have yielded promising results. Moreover while one may often argue that for the purposes of an explanatory analysis, compared with underfitting, overfitting is the less serious error, this is perhaps more justified in a multi-stage setting where mis-specification by underfitting can have a more severe knock-on effect in our analysis. When the outcome was generated using normal errors, rather than the skewed log-normal errors discussed here, the Wald-type approaches become more competitive, but still badly underfit (full results are included in the Supplementary Material). As in any simulation-based analysis we appreciate that the results presented here cannot possibly be comprehensive, and so would encourage further experimentation (and analysis) to assess the properties of these various criteria. Of particular interest is a more extensive investigation of the QICG for discrete outcomes, with our preliminary analyses providing encouraging results.

The trace term, $K$, warrants additional research. In our simulations we observed that when errors were normally distributed, and the blip model correctly specified, $K$ was approximately equal to the dimension of the model (as is the case in the likelihood-based setting). This was not generally the case, however, especially when both the treatment and treatment-free models were badly mis-specified or errors were non-normal. From a practical perspective, the quasi-likelihood criteria presented may tend to underfit when models are mis-specified (as the corresponding penalty term is large), but if the researcher has doubts about the legitimacy of parameter estimators (due to mis-specification of both the treatment and treatment-free models), then model selection is a secondary concern. In our analysis of the STAR*D data we found that for the various candidate models, the trace term was somewhat larger than the dimension of the associated blip function while residual plots from the proposed models were consistent with normal errors. This could indicate an inadequacy in the treatment or treatment-free models which might merit further investigation.

Our QIC formulation has focused exclusively on selection of the blip, or contrast, component of the outcome mean model. However, the approach can easily be adapted to select the entire mean model (the contrast and treatment-free models simultaneously), allowing the use of our information criterion for G-estimation of static treatment sequences or G-estimation for mediation, as well as in binary outcome settings.

7 Supplementary Material

Appendix:

Appendix containing additional theory, and simulation results.

R code:

R code for the simulation studies of section 5 and the Appendix.

Appendix: Supplementary Material for Manuscript titled “Generalized G-estimation and Model Selection”

8 Quasi-likelihood for continuous treatments

In the binary treatment setting we showed that our estimating equations ${U}(\bm{\psi})$ may be reduced to

[TABLE]

We now extend to the case of a continuous treatment and include a quadratic term in our blip such that

[TABLE]

where we have compartmentalized our blip parameters $\bm{\psi}=(\bm{\psi}_{1},\bm{\psi}_{2})$ and history design matrix $\bm{h}_{\psi}=(\bm{h}_{\psi_{1}},\bm{h}_{\psi_{2}})$ depending on whether they are associated with the linear or quadratic term of ${a}$ in our blip. Writing $\bm{D}_{1}$ and $\bm{D}_{2}$ for the diagonal matrices with $(i,i)^{th}$ entry $a_{i}-{E}[A_{i}|H_{i}]$ and $a_{i}^{2}-{E}[A_{i}^{2}|H_{i}]$, respectively, our estimating equations become (Rich et al. 2014)

[TABLE]

, which yields a quasi-likelihood of the form

[TABLE]

where $\bm{m}_{c}$ and $\bm{M}_{c}$ may be thought of as the continuous treatment analogs to $\bm{m}$ and $\bm{M}$ derived in the main paper for the binary case.

9 Proof of discrepancy theorem

Recall the theorem from section 2.4 of the associated paper:

Theorem: Suppose that $Q(\bm{\psi})$ is twice continuously differentiable with bounded expectation of its second derivative in a neighbourhood $\mathcal{N}$ of $\bm{\psi}_{(m,*)}$. Then, under the stable unit treatment value and no unmeasured confounding assumptions (detailed in section 2.1), the expected divergence $\Delta(m)$ can be approximated

[TABLE]

which is consistently estimated by

[TABLE]

where $\bm{I}(.)$ and $\bm{J}(.)$ are the observed (empirical) versions of $\mathcal{I}$ and $\mathcal{J}$. Thus, the model selection procedure that chooses a model by minimizing $QIC_{G}(m)$ across $\mathcal{M}(m)$ identifies the model that minimizes $\Delta(m)$ with probability 1 as $n\longrightarrow\infty$.

Proof: Following Takeuchi (1976) – see also Burnham and Anderson (2002) – our quasi-likelihood information criterion (QIC) is based on an estimate of $\Delta(m)$, a function of the Kullback-Leibler discrepancy between the data generating model and model $m$. Consider a decomposition of $\Delta(m)=E[\delta(\widehat{\bm{\psi}}_{(m)})]$, where

[TABLE]

given by

[TABLE]

We consider first an expansion of $Q(.)$ around the “true” blip parameter for model $m$, ${\bm{\psi}}_{(m,*)}$: we have that

[TABLE]

so taking expectations with respect to the data generating model with $\widehat{\bm{\psi}}_{(m)}$ fixed, we have

[TABLE]

where, recall,

[TABLE]

and hence

[TABLE]

Using the conventional theory of misspecified models, we have that

[TABLE]

where, by construction of the quasi-likelihood function, we have

[TABLE]

Hence

[TABLE]

as $n\longrightarrow\infty$, where

[TABLE]

and

[TABLE]

By standard convergence results

[TABLE]

and by a standard result for quadratic forms

[TABLE]

Under standard regularity conditions on $Q(.)$, we have that $\widehat{\bm{\psi}}_{(m)}\stackrel{{\scriptstyle p}}{{\longrightarrow}}{\bm{\psi}}_{(m,*)}$ as $n\longrightarrow\infty$, and hence for large $n$ $E[-2Q(\widehat{\bm{\psi}}_{(m)})]=E[-2Q({\bm{\psi}}_{(m,*)})]+o(1)$, so

[TABLE]

This completes the proof.

The asymptotic variance of estimator $\widehat{\bm{\psi}}_{(m)}$ may itself be estimated by

[TABLE]

10 Investigating the trace term

We examine in simulation whether the trace term is a good approximation to the true dimension of the underlying model when that model is fitted, as described in section 5.3 of the main paper. Figure 1 displays estimates of stage 1 (top) and stage 2 (bottom) trace term $K$ from 1,000 simulations with $n=100$ for different candidate models (y-axis). The panels correspond to true blip models containing an intercept term along with stage $j$ covariates $x_{j1}$ (left), $x_{j1},x_{j2}$ (middle), $x_{j1},x_{j2},x_{j3}$ (right), and gray boxes correspond to simulations where the true blip model was fit. In expectation, the trace term matches the dimension of the data generating model even in this relatively small sample setting.

11 Simulation results

11.1 Discrete Outcome

Next, we demonstrate the use of QICG in the DTR framework for a discrete outcome case. We generate data as follows:

• •

stage 1 patient information: $X_{11}\sim N(1,1)$, $X_{12}\sim N(-1,1)$, $X_{13}\sim N(1,1)$;

• •

stage 1 treatment: $a_{1}\in\{0,1\}$, ${P}(A_{1}=1|\bm{h}_{1})=\text{expit}(x_{11})$;

• •

stage 2 patient information: $X_{21}\sim N(a_{1},1)$, $X_{22}\sim N(-1,1)$, $X_{23}\sim N(1,1)$;

• •

stage 2 treatment: $a_{2}\in\{0,1\}$, ${P}(A_{2}=1|\bm{h}_{2})=\text{expit}(x_{21})$;

• •

stage $j$ blip: $\gamma_{j}(a_{j},\bm{h}_{j})=a_{j}(0.5+\psi_{j1}x_{j1}+\psi_{j2}x_{j2}+\psi_{j3}x_{j3})$

such that $a_{j}^{opt}=1$ if $0.5+\psi_{j1}x_{j1}+\psi_{j2}x_{j2}+\psi_{j3}x_{j3}>0$ and 0 otherwise;

• •

outcome: $P(Y=k)=\lambda^{k}e^{-k}/k!$,

with $\lambda=\exp\left[\beta_{0}-\sum_{j=1}^{2}[\gamma_{j}(a_{j}^{opt},\bm{h}_{j})-\gamma_{j}(a_{j},\bm{h}_{j})]\right]$,

where we vary $\beta_{0}$ so that for the various $\psi_{jk}$ we consider, $P(Y=0)=0.1$. We set the blip parameters to $(\psi_{j1},\psi_{j2},\psi_{j3})=(0.5,0,0)$, $(0.5,0.5,0)$ or $(0.5,0.5,0.5)$ giving a range of models including one, two, or all three variables at each stage.

Given the computational requirements for the IRLS algorithm, our simulations are somewhat more limited. We restrict ourselves to sample sizes of $n=200$, and only consider the four blip models containing no covariates, $x_{j1}$ only, $x_{j1}$ and $x_{j2}$ only, or all three covariates, and choose whichever model resulted in the lowest $\text{QIC}_{G}$. As with our simulations in the main paper (Table 1), we correctly specify our treatment model, and mis-specify the treatment-free model, supposing it is linear in $x_{11}$ at stage 1 and linear in $x_{21}$ at stage 2. Results are summarized in Table 4, where the stage 1 results are based on analyses where the stage 2 blip models were correctly specified. The IRLS algorithm was implemented with an iteration limit of 1,000 and a tolerance limit (between successive log-mean function estimates) of 0.001.

The presented results summarize over the simulation runs where all four IRLS algorithms converged. This corresponds to all 1,000 simulated datasets for the simplest models (with only one covariate in the true blip model at each stage), 996 (stage 2) and 994 (stage 1) datasets when two covariates were included, and 904 (stage 2) and 863 (stage 1) datasets where all three were included. If we instead presume model selection is based on the lowest QIC only for those candidate blip models where IRLS converged (as might occur in a real-life analysis), results were near-identical with the slight exception of model selection for the most complex blip model at both stages. In this case, the correct model selection rate drops from 0.921 to 0.833 for stage 2, and from 0.919 to 0.877 for stage 1, with the largest candidate model most likely to result in failure.

11.2 Continuous Outcome

Additional simulation results referred to in the main paper. Table 5 contains results for varying effect sizes, Table 6 contains results for varying correlation strength between covariates, and Table 7 contains non-aggregated results. Tables 8-12 contain analogous results to Tables 1-2 in the main paper and Tables 5-7 in this appendix, but with standard normal errors in the generation of the outcome $Y$, rather than log-normal errors. In all tables bold indicates the most successful method for each scenario.