When and when not to use optimal model averaging

Abstract

Traditionally model averaging has been viewed as an alternative to model selection with the ultimate goal to incorporate the uncertainty associated with the model selection process in standard errors and confidence intervals by using a weighted combination of candidate models. In recent years, a new class of model averaging estimators has emerged in the literature, suggesting to combine models such that the squared risk, or other risk functions, are minimized. We argue that, contrary to popular belief, these estimators do not necessarily address the challenges induced by model selection uncertainty, but should be regarded as attractive complements for the machine learning and forecasting literature, as well as tools to identify causal parameters. We illustrate our point by means of several targeted simulation studies.

Appendices

Appendix A: Software

We have implemented Mallow’s model averaging, Jackknife model averaging, and Lasso averaging in the R-package MAMI (Schomaker 2017a), available at http://mami.r-forge.r-project.org/. In addition to this, we have implemented several wrappers that make optimal model averaging easily useable for super learning (Polley et al. 2017), and in conjunction with causal inference packages such as tmle (Gruber and van der Laan 2012) and ltmle (Lendle et al. 2017). Available wrappers are explained by calling listSLWrappers(), and examples are given in the documentation (Schomaker 2017b).

Appendix B: Notation and background on the sequential g-formula

Consider a sample of size n of which measurements are available both at baseline ($t=0$) and during a series of follow-up times $t=1,\ldots ,T$. At each point in time we measure the outcome $Y_t$, the intervention $A_t$, time-dependent covariates $\mathbf {L}_t=\{L^1_t,\ldots ,L^q_t\}$, and a censoring indicator $C_t$. $\mathbf {{L}}_{t}$ may include baseline variables $\mathbf {V}=\{L^1_0,\ldots ,L^{q_V}_0\}$ and can potentially contain variables which refer to the outcome variable before time t, for instance ${Y}_{t-1}$. The treatment and covariate history of an individual i up to and including time t is represented as $\bar{A}_{t,i}=(A_{0,i},\ldots ,A_{t,i})$ and $\bar{L}^s_{t,i}=(L^s_{0,i},\ldots ,L^s_{t,i})$ respectively. $C_t$ equals 1 if a subject gets censored in the interval $(t-1,t]$, and 0 otherwise. Therefore, $\bar{C}_t=0$ is the event that an individual remains uncensored until time t.

The counterfactual outcome $Y_{t,i}^{\bar{a}_{t}}$ refers to the hypothetical outcome that would have been observed at time t if subject i had received, possibly contrary to the fact, the treatment history $\bar{A}_{t,i}=\bar{a}_t$. Similarly, $\mathbf {L}_{t,i}^{\bar{a}_{t}}$ are the counterfactual covariates related to the intervention $\bar{A}_{t,i}=\bar{a}_t$. The above notation refers to static treatment rules; a treatment rule may, however, depend on covariates, and in this case it is called dynamic. A dynamic rule ${d}_t({\bar{\mathbf {L}}}_{t})$ assigns treatment ${A}_{t,i} \in \{0,1\}$ as a function of the covariate history ${\bar{\mathbf {L}}}_{t,i}$. The vector of decisions $d_t, t=0,\ldots ,T$, is denoted as $\bar{d}_T=\bar{d}$. The notation $\bar{A}_t = \bar{d}$ refers to the treatment history according to the rule $\bar{d}$. The counterfactual outcome related to a dynamic rule $\bar{d}$ is $Y_{t,i}^{\bar{d}}$, and the counterfactual covariates are $\mathbf {L}_{t,i}^{\bar{d}}$.

In Sect. 3.3 we consider the expected value of Y at time 6, under no censoring, for a given treatment rule $\bar{d}$ to be the main quantity of interest, that is $\psi =\mathbb {E}(Y_6^{\bar{d}})$.

The sequential g-formula can estimate this target quantity by sequentially marginalizing the distribution with respect to $\mathbf {L}$ given the intervention rule of interest. It holds that

$$\begin{aligned} \mathbb {E}(Y_T^{\bar{d}})= & {} \mathbb {E}(\,\mathbb {E}(\,\ldots \mathbb {E}(\,\mathbb {E}(Y_T|\bar{A}_T=\bar{d}_{T}, {\bar{\mathbf{L}}}_T) | \bar{A}_{T-1}=\bar{d}_{T-1}, {\bar{\mathbf{L}}}_{T-1}\,)\ldots |\bar{A}_{0}={d}_{0}, \mathbf {{L}}_{0}\,)|\mathbf {{L}}_{0}\,), \end{aligned}$$

see for example Bang and Robins (2005). Equation (14) is valid under several assumptions: sequential conditional exchangeability, consistency, positivity and the time ordering $\mathbf {L}_t \rightarrow A_t \rightarrow C_t \rightarrow Y_t$. These assumptions essentially mean that all confounders need to be measured, that the intervention is well-defined and that individuals have a positive probability of continuing to receive treatment according to the assigned treatment rule, given that they have done so thus far and irrespective of the covariate history; see Daniel et al. (2013) and Robins and Hernan (2009) for more details and interpretations. Note that the two interventions defined in Sect. 3 also assign $C_t=0$ meaning that we are interested in the effect estimate under no censoring.

To estimate $\psi$ one needs to the the following for $t=T,\ldots ,0$: (i) use an appropriate model to estimate $\mathbb {E}(Y_T|{\bar{A}}_{T-1}={\bar{d}}_{T-1}, {\bar{\mathbf{L}}}_T)$. The model is fit on all subjects that are uncensored (until $T-1$). Note that the outcome refers to the measured outcome for $t=T$ and to the prediction (of the conditional outcome) from step (ii) if $t<T$. Then, (ii) plug in $\bar{A}_t=\bar{d}_{t}$ to predict $Y_t$ at time t; (iii) For $t=0$ the estimate $\hat{\psi }$ is obtained by calculating the arithmetic mean of the predicted outcome from the second step.

Appendix C: Data generating process in the simulation study

Both baseline data ($t=0$) and follow-up data ($t=1,\ldots ,12$) were created using structural equations using the R-package simcausal (Sofrygin et al. 2017). The below listed distributions, listed in temporal order, describe the data-generating process motivated by the analysis from Schomaker et al. (2016). Baseline data refers to region, sex, age, CD4 count, CD4%, WAZ and HAZ respectively ($V^1$, $V^2$, $V^3$, $L^1_0$, $L^2_0$, $L^3_0$, $Y_0$), see Schomaker et al. (2016) for a full explanation of variables and motivational question. Follow-up data refers to CD4 count, CD4%, WAZ and HAZ ($L^1_t$, $L^2_t$, $L^3_t$, $Y_t$), as well as a treatment ($A_t$) and censoring ($C_t$) indicator. In addition to Bernoulli (B), uniform (U) and normal (N) distributions, we also use truncated normal distributions which are denoted by $N_{[a,b]}$ where a and b are the truncation levels. Values which are smaller a are replaced by a random draw from a $U(a_1,a_2)$ distribution and values greater than b are drawn from a $U(b_1,b_2)$ distribution. Values for $(a_1, a_2, b_1, b_2)$ are (0, 50, 5000, 10000) for $L^1$, (0.03,0.09,0.7,0.8) for $L^2$, and $(-10,3,3,10)$ for both $L^3$ and Y. The notation $\bar{\mathcal {D}}$ means “conditional on the data that has already been measured (generated) according the the time ordering”.

For $t=0$:

$$\begin{aligned} V^1\sim & {} B(p=4392/5826) \\ V^2|\bar{\mathcal {D}}\sim & {} \left\{ \begin{array}{ll} B(p=2222/4392) &{} \quad \text{ if } \quad V^1 = 1\\ B(p=758/1434) &{} \quad \text{ if } \quad V^1 = 0\\ \end{array} \right. \\ V^3|\bar{\mathcal {D}}\sim & {} U(1,5) \\ L^1_0|\bar{\mathcal {D}}\sim & {} \left\{ \begin{array}{cl} N_{[0,10000]}(650,350) &{} \quad \text{ if } \quad V^1 = 1\\ N_{[0,10000]}(720,400)) &{} \quad \text{ if } \quad V^1 = 0\\ \end{array} \right. \\ \tilde{L}^1_0|\bar{\mathcal {D}}\sim & {} N((L^1_0-671.7468)/(10\cdot 352.2788)+1,0)\\ L^2_0|\bar{\mathcal {D}}\sim & {} N_{[0.06,0.8]}(0.16+0.05\cdot (L^1_0-650)/650,0.07) \\ \tilde{L}^2_0|\bar{\mathcal {D}}\sim & {} N((L^2_0-0.1648594)/(10\cdot 0.06980332)+1,0)\\ L^3_0|\bar{\mathcal {D}}\sim & {} \left\{ \begin{array}{ll} N_{[-5,5]}(- 1.65 + 0.1 \cdot V^3 + 0.05 \cdot (L^1_0 - 650)/650 \\ +\, 0.05 \cdot (L^2_0 - 16)/16,1) &{} \quad \text{ if } \quad V^1 = 1\\ N_{[-5,5]}( -2.05 + 0.1 \cdot V^3 + 0.05 \cdot (L^1_0 - 650)/650 \\ +\, 0.05 \cdot (L^2_0 - 16)/16,1)) &{} \quad \text{ if } \quad V^1 = 0\\ \end{array} \right. \\ A_0|\bar{\mathcal {D}}\sim & {} B(p=0) \\ C_0|\bar{\mathcal {D}}\sim & {} B(p=0) \\ Y_0|\bar{\mathcal {D}}\sim & {} N_{[-5,5]}(-2.6 + 0.1 \cdot I(V^3 > 2) + 0.3 \cdot I(V^1 = 0) + (L^3_0 + 1.45),1.1) \\ \end{aligned}$$

For $t>0$:

$$\begin{aligned} L^1_t|\bar{\mathcal {D}}\sim & {} \left\{ \begin{array}{cl} &{} N_{[0,10000]}(13\cdot \log (t \cdot (1034-662)/8) + L^1_{t-1} + 2 \cdot L^2_{t-1}\\ &{} +\, 2 \cdot L^3_{t-1} + 2.5 \cdot A_{t-1},50) \qquad \text{ if } \quad t \in \{1,2,3,4\}\\ &{} N_{[0,10000]}(4\cdot \log (t \cdot (1034-662)/8) + L^1_{t-1} + 2 \cdot L^2_{t-1} \\ &{}+\, 2 \cdot L^3_{t-1} + 2.5 \cdot A_{t-1},50) \qquad \text{ if } \quad t \in \{5,6\}\\ \end{array} \right. \\ L^2_t|\bar{\mathcal {D}}\sim & {} N_{[0.06,0.8]}(L^2_{t-1} + 0.0003 \cdot (L^1_t - L^1_{t-1}) + 0.0005 \cdot (L^3_{t-1}) + 0.0005 \cdot A_{t-1} \cdot \tilde{L}^1_0,0.02) \\ L^3_t|\bar{\mathcal {D}}\sim & {} N_{-5,5}(L^3_{t-1} + 0.0017 \cdot (L^1_t - L^1_{t-1}) + 0.2 \cdot (L^2_t - L^2_{t-1}) + 0.005 \cdot A_{t-1} \cdot \tilde{L}^2_0,0.5) \\ A_t|\bar{\mathcal {D}}\sim & {} \left\{ \begin{array}{ll} &{}B(p=1) \qquad \text{ if } \quad A_{t-1} = 1\\ &{} B(p=1/(1+\exp (-[-2.4 + 0.015 \cdot (750 - L^1_t) + 5 \cdot (0.2 - L^2_t) \\ &{} -\, 0.8 \cdot L^3_t + 0.8 \cdot t]))) \qquad \text{ if } \quad A_{t-1} = 0\\ \end{array} \right. \\ C_t|\bar{\mathcal {D}}\sim & {} B(p=1/(1+\exp (-[-6+ 0.01 \cdot (750 - L^1_t) + 1 \cdot (0.2 - L^2_t) - 0.65 \cdot L^3_t - A_t]))) \\ Y_t|\bar{\mathcal {D}}\sim & {} N_{[-5,5]}(Y_{t-1} + 0.00005 \cdot (L^1_t - L^1_{t-1}) - 0.000001 \cdot \left( (L^1_t - L^1_{t-1})\cdot \sqrt{\tilde{L}^1_0}\right) ^2\\&+\, 0.01 \cdot (L^2_t - L^2_{t-1})- 0.0001 \cdot \left( (L^2_t - L^2_{t-1})\cdot \sqrt{\tilde{L}^2_0}\right) ^2\\&+\, 0.07 \cdot ((L^3_t-L^3_{t-1})\cdot (L^3_0+1.5135)) - 0.001 \cdot ((L^3_t-L^3_{t-1})\cdot (L^3_0+1.5135))^2 \\&+\, 0.005 \cdot A_t + 0.075 \cdot A_{t-1} + 0.05 \cdot A[t] \cdot A[t-1] ,0.01) \\ \end{aligned}$$