Influence Re-weighted G-Estimation

2.1 G-estimation of SNMMs for DTRs

A dynamic treatment regimen is a sequence of functions d=(d1,…,dK) where K is the number of treatment intervals and for each interval j=1,…,K the function dj:Hj ↦ Aj maps the set of possible observed histories at interval j to the set of possible treatments at interval j. Letting Aj denote treatment received at interval j, Lj denote observations made subsequent to receiving the treatment but prior to receiving the next treatment, and L0 baseline covariates, an observed history is given by Hj=(L0,A1,L1,A2,…,Lj−1). In this longitudinal setting, the presence of qualitative interactions dictates that the best treatment strategy involves individualization based on the observed patient history up to the point at which the decision is made. To estimate such a strategy, Robins [9] proposed the SNMM. This model is specified via a blip function, which models the interactions between treatment and covariates in the following way. Let Y=g(O) be the final outcome of interest, where g(⋅) is a suitably defined reward function (so that greater values of Y are more desirable). The blip function is defined in terms of potential outcomes (or counterfactuals); let Y(d) denote the potential outcome under DTR d. Let dhj,a,opt denote the DTR consisting of a1,a2,…,aj−1 (components of hj) for the first (j−1) treatment intervals, then a at the jth interval, then the optimal DTR for the remaining (K−j) intervals (it is understood that either (j−1) or (K−j) may be zero, in which case either the component before or after a becomes unnecessary). Let 0 denote a reference treatment, such as no treatment or placebo. The blip function for the jth treatment interval is given by

Two key features of the blip function are (1) the blip function evaluated at a=0 is zero (by definition); (2) the blip function determines the optimal DTR through dopt(hj)=argmaxaγj(a,hj). Under standard identifiability and regularity assumptions given in Robins [9], the blip function can be estimated. A SNMM sets γj(a,hj) equal to some function of a parameter ψj and covariates xj=xj(hj), for instance γj(a,hj;ψj)=ψj⊤⋅axj. (Typically, the blip function is assumed correctly specified, but if not, then the parameter psij would correspond to a projection of the truth onto the specified model, resulting in a DTR that is optimal within the class defined by the model.) Let

where ψ_j+1=(ψj+1⊤,…,ψK⊤)⊤. Doubly robust estimators of the interval-specific parameters ψj are obtained (under regularity conditions) by recursively solving the g-estimating equations

where πj(hj;αj) is a posited or working model for E[Aj|Hj=hj]=Pr[Aj=1|Hj=hj], ηj(hj;ςj) is a posited or working linear model for E[Gj(ψj,ψ_j+1)|Hj,Aj]=E[Y(dhj,0,opt)|Hj,Aj], and ςˆj(ψj,ψ_j+1) is the ordinary least-squares estimator of ςj, determined up to a parameter.

We shall now consider a partially misspecified setting, namely we shall suppose that while the treatment or propensity model πj(hj;αj) is correctly specified, the treatment-free outcome model ηj(hj;ςj) is not. As noted previously, the resulting g-estimator ψˆj is CAN but may be susceptible to high variation due to influential observations. This is the point we aim to address. Model diagnostics can help. In a previous paper [14], we proposed residual-based diagnostics for g-estimation of optimal dynamic treatment regimens which can serve as a guide to obtaining a better fitting outcome nuisance model, and also showed in simulations that improved nuisance model fit resulted in improved estimator efficiency for the parameter of interest. In this paper, we will pursue a different approach for obtaining more robust (in the sense of Hampel [15]) g-estimators in the same context through re-weighting. Here we will assume that the model specification is fixed, i.e. cannot be improved through diagnostics, and focus on the properties of the g-estimator for this fixed model under a varying set of weights.

2.2 Weighting in estimating equations: some general results

Suppose that we have observed i.i.d. data O1,…,On from a population, and let O denote a generic observation for a single individual from the same population, and o a realization. Let β be a parameter of interest and let u(o,β) be a consistent estimating function, i.e. we have E[u(O,β0)]=0, where β0 denotes the truth (we will suppose that β0 is the unique solution to this equation, which in the present context is true if models γj(a,hj;ψj) and ηj(hj;ςj) are linear in parameters).

Standard M-estimator theory Tsiatis [16] tells us that we can construct an (unweighted) estimator of β that is regular asymptotically linear (RAL) and hence consistent asymptotically normal as the solution to 0=∑i=1nu(Oi,βˆn). Under regularity conditions, βˆn is consistent, has influence function given by:

and asymptotic variance given by MVar(u(O,β0))MT where

Now, let w(Oi) be a strictly positive weight function that does not depend on β. We will use the more concise notation Wi=w(Oi), and Ui(⋅)=u(O,⋅). Note that the weights Wi are data-dependent; this generalized form of M-estimation does require additional technical conditions, which we shall discuss following Lemma 1 below. Then, a weighted estimator βˆnw is given by the solution to the estimating equation 0=∑i=1nWiUi(βˆnw). We will explore the properties of weighted estimators in terms of bias and variance.

First, we observe that the usual Taylor expansion gives

where β∗ lies between β0 and βˆnw. So

Where Z has the standard normal distribution.

If βˆnw happens to be n-consistent, i.e. if n(βˆnw−β0)→dN(0,Var(βˆnw)), then the same M-estimator theory applies to βˆnw as to βˆn, and consequently, under regularity conditions, the asymptotic variance of βˆnw is given by MwVar(WiUi)MwT where

We now give a sufficient condition for the consistency of βˆnw:

Proof. If the weights are not data-dependent, then by iterated expectation:

so the usual M-estimator theory also applies to βˆnw. As the weights are, in general, data-dependent, the proof can be made under the limiting values of Wi following the consistency requirements for the estimators used to define Wi using the arguments for adaptive weighted estimators outlined in Section 4.4 of Newey [17].

3.1 Influence-based weighting schemes

Under correct specification of the treatment-free outcome model, Robins [9] gives an expression for a weight function for g-estimation that is optimal in an efficiency sense (under some further assumptions). The weight is proportional to the inverse of the conditional variance Var(Gj(ψj_)|Hj). Under a misspecified treatment-free outcome model, however, this weight function would be unavailable. Joffe and Brensinger [18] discuss weighting in g-estimation, but in the context of an instrumental-variables analysis for estimating treatments effects in randomized trials with non-compliance. Focussing on the single-interval case, their proposal of weighting by the compliance score assigns larger weights to compliers. Some optimality results are given, but again only under correct model specification.

We propose a different approach to constructing weights which seeks to stabilize the standard estimator in a manner which does not depend on correct model specification and can improve its performance as measured by mean squared error in certain cases (as we will see in a simulation study of Section 4).

Data points that are candidates to receive reduced weights in the estimating equation are points with high sample influence. Hinkley [19] and Krasker and Welsch [20] discuss the use of weights constructed as functions of the sample influence in the context of robust regression estimators. Here, the proposal is to identify data points that are highly influential in the estimation of the parameter of interest, and reducing the influence of those individuals through subsequent re-weighting, will stabilize the estimation and make it more robust to misspecification of the nuisance treatment-free outcome model. The sample influence of a data point is, however, constructed from all the available data. Thus, while it is straightforward to describe an algorithm in which weights are constructed such that the weight for individual I is inversely proportional to the sample influence of individual i on the estimate of ψj (as we now proceed to do), the asymptotic properties of the resulting estimator are difficult to assess the terms of the estimating equation are not i.i.d. since the weights depend on all of the available data and the treatment mechanism and outcome regression models depend on all data (although the weight functions converge to functions which only rely on the individual’s history). As the terms are not independent, the standard weighting theory described above does not apply and the resulting estimator will not in general be consistent for the parameter of interest. We begin, for simplicity, by ignoring these issues, then address them with a simple adaptation of the influence re-weighting g-estimation procedure.

3.2 Influence measures in g-estimation

Influence refers to the stability of statistics (e.g. estimators) to perturbations of the data. The mathematical object used to describe the asymptotic stability of an estimator is the influence function (IF) (Hampel et al. [15], Tsiatis [16], Cook and Weisberg [21]). The idea is to consider the distribution function F from which the observed data Z1,…,Zn are independently drawn, a statistic T(F), and form the directional derivative

where δz represents a point-mass density at z. Thus, the IF ϕT,F(z) of T with respect to F describes the effect on T of an infinitesimal “contamination” at z [15].

The IF is a useful theoretical tool, but its use in practice for model diagnostics is limited since it depends on the true data generating distribution, which is only available under correct model misspecification and infinite samples. More useful for diagnostics is the finite sample analogue of the IF, obtained through the jackknife (leave-one-out) approach. Consider Fˆ(i), the empirical distribution function for a data set from which the ith observation has been removed, thus

Then, the sample influence (SI) of observation i is obtained from (2) by setting F=Fˆ(i), z=zi and ϵ=1/n, given

In standard outcome regression such as E[Y|X]=m(X;β), where Z=(X,Y), the usual measures derived from the sample influence (e.g. Seber and Lee [22]) include:

1. DFBETAi [23], defined by:

   βˆ−βˆ(i);

2. DFFITi [23], defined by:

   m(xi;βˆ)−m(xi,βˆ(i));

3. Cook’s distance Di [21], defined by:

   (βˆ−βˆ(i))TΣ−1(βˆ−βˆ(i))

   for a symmetric and positive definite matrix Σ, typically a scaled estimate of the covariance matrix of βˆ.

Consider now the g-estimator ψˆj. Let ψˆj(i) be the estimate of ψj obtained after deletion of individual i. Then, the sample influence of individual i on the individual components of the estimate ψˆj are given by ψˆj−ψˆj(i), which is equivalent to DFBETA above. We will also be interested in the overall or compound influence of individual i, defined in a manner that is analogous to Cook’s Di above, but where we consider the square root:

and Σˆψj is an estimate of the covariance matrix of ψˆj. There are no computational shortcuts as in the case of linear regression, where elegant mathematics allow the deletion diagnostics to be computed very efficiently [21], here the computation is by brute force. Nonetheless, with closed form g-estimation and fast computers, the computational burden is modest for reasonable sample sizes.

3.3 Influence re-weighted g-estimation

The procedure begins with an initial solution to the weighted g-estimating equations, using some initial set of weights that we denote Wij0 (note that we are free to set Wij0=1 for all i,j). Based on this initial solution, the sample influence of each data point at each interval is computed using the usual leave-one-out approach as described in previous section. Let Iij denote the influence of the observation made on individual i at interval j. For instance, let

where ψˆj(i) denotes the estimate of ψj obtained when the ith individual is deleted from the data set (or equivalently, assigned weight zero), and Σˆψj is an estimate of the covariance matrix of ψˆj. Define a monotone non-increasing function ρ:ℝ+→ℝ+. For instance, we propose the family of functions

where λ=(λ1,λ2)∈ℝ+×ℝ+ is a tuning parameter. The function ρλ(x) incorporates a “scale” parameter (λ1) and a “shape” parameter (λ2) is smooth and bounded above at 1 so that weights cannot be absurdly large. The procedure assigns to individual i at interval j the weight Wij=Wij0×ρ(Iij). The g-estimating equations are re-solved using this set of weights and a new estimate of ψ is obtained. We refer to an estimator constructed in this manner as an influence-based weight (IBW) g-estimator.

In the proposed family ρλ, the tuning parameter λ determines the shape of the “down-weighting curve”, i.e. how quickly the weights of points with high influence are shrunk down to zero. While the choice of the λ should be application-specific and is left to the analyst, some guidance can be provided. The consideration is the trade-off between bias and variance. To gain a sense for the bias and variance of the weighted estimator under different choices of λ, bootstrap re-sampling could be used. A less computationally intensive approach is to evaluate the estimator under different choices on a single data set, as Cole and Hernán [25] do in their discussion of the bias-variance trade-off involved in the selection of a truncation point for weights in IPTW estimation. For example, we could study the impact of choosing λ such that the median weight is equal to 0.5, 0.55 0.6,…, 0.95 when the weights have been normalized such that the largest weight is 1.

3.4 Re-weighting “inside-the-loop” and “outside-the-loop”

When g-estimation is applied recursively, there are two different variants of the algorithm for computing a re-weighted estimator depending on the ordering of the calculations. We have called them “inside-the-loop” and “outside-the-loop”, making reference to a programming loop with a loop index j that goes down from K to 1. The two variants are sketched in pseudo-code in Figure 1. Essentially, the distinction is whether the initial estimate of ψj uses the weighted or unweighted estimate of ψm for intervals j<m. In the “inside-the-loop” variant, the weighted estimates at later intervals are computed first, and used in the calculation of the initial estimate of ψj. In the “outside-the-loop” variant, all initial estimates are computed before any of the weighted estimates. In the simulations, the “inside-the-loop” variant was used because we expect the weighted estimates of ψm for m>j to be slightly more stable than the unweighted estimates, and hence lead to a slightly better initial estimate of ψj.

Figure 1:

Algorithms for re-weighting “inside-the-loop” (a) and “outside-the-loop” (b).

3.5 Consistent re-weighting

As stated previously, the estimation procedure described above may have undesirable properties. In particular, because the terms of the estimating equations are not i.i.d. and the conditions of Lemma 1 do not necessarily hold, consistency of the estimator is not assured. We now propose a simple way to adapt the procedure and obtain a different estimator, similarly motivated, that is not subject to these difficulties.

The idea is to construct a set of weights such that the weight for individual i at interval j is a function of the history Hij only. Consistency then follows from Lemma 1, Equation (3), and regularity conditions on the weights. To do so, a predictive model for the sample influence Iij using Hij as predictors is developed. To compute the new weights, the function ρ is then applied to the predicted sample influence rather than actual sample influence. We will refer to an estimator that uses weights based on a predicted influence as a predicted influence-based weight (PIBW) estimator.

The challenge with this approach is in constructing predictive models for the sample influence at each interval. We approached this in two different ways: direct prediction of the compound influence Iij, and prediction of the individual components of ψˆj(i). For direct prediction of the compound influence Iij, we propose to regress logIj on Hj using some possibly flexible model, and take the exponential of the predicted value as Iˆij(hj). The log transformation is appropriate because Iij is a non-negative quantity. We refer to this approach as PIBW1.

Unlike the compound influence Iij which is non-negative, the individual components of ψˆj(i) can take any real value so no transformation is needed and we suppose that the relationship between covariates and components of ψˆj(i) will be more linear, and that a less flexible model will be required for prediction. However, multiple models are required since each component must be predicted separately. Once these separate predictions have been computed, I^ij(hj) is obtained using (3) but with each component of ψˆj(i) replaced by its prediction. We refer to this approach as PIBW2.

Yet a third approach that we considered was to orthogonalize the individual components before prediction. This allows for better prediction in the case where the individual components of ψˆj(i) are highly correlated. Specifically, let qij=Σˆψj−1/2(ψˆj(i)−ψˆj) where Σψj−1/2 is the inverse of the Choleski decomposition of the estimate of the covariance matrix of ψˆj. Then, regress each component of the vectors qij on Hj to obtain predictions qˆij(Hij), and finally set Iˆij(Hij)=qˆij(Hij)Tqˆij(Hij). We refer to the approach of orthogonalizing before predicting as PIBW3.

Once the prediction Iˆij(Hij) has been obtained (using any of the three versions described above), the weights are computed analogously to the IBW case, namely Wij(Hij)=Wij0×ρ(Iˆij(Hij)). For inference, the variance estimator for ψˆj can be adjusted to account for the estimation of an additional nuisance model in the usual way (i.e. following the same approach used to adjust for the estimation of the parameter αj in the propensity score model).