Proposer of the vote of thanks to Evans and Didelez and contribution to the Discussion of ‘Parameterizing and simulating from causal models’
Shaun R. Seaman · Journal of the Royal Statistical Society Series B (Statistical Methodology) · 2024
Simulation studies are a key tool for assessing performance of statistical methods. It is common in such studies to generate data in a way that makes the model of interest correctly specified. This is not always straightforward when the model of interest is a causal/structural model, i.e. model for potential outcomes. Evans and Didelez (ED) outline a general solution to this problem. They first simulate data that would arise if exposure were randomized, and then use rejection sampling to obtain ‘observational’ data, i.e. data where exposure is not randomized and there is confounding. This rejection sampling step can be viewed as the opposite of the inverse propensity score weighting (IPW) commonly used when analysing observational data: rejection sampling creates the very propensity score weighting that the analyst uses IPW to eliminate. Superscripts will denote potential random variables under an intervention, e.g. Yx is the variable Y when we set X=x, and FV and fV will denote distribution function and probability density/mass function of generic variable(s) V. Fundamental to ED’s approach is the factorization of the joint distribution. In the scenario of Figure 1a, joint distribution FZ,X,Yx is factorized as the product of FZ,Yx and FX∣Z, and FZ,Yx is factorized in terms of marginal distributions FZ and FYx and some specification of the association between Z and Yx. Note that (a) Z=Zx, because Z is causally prior to X and (b) FYx is at least partly defined by the causal model, e.g. marginal structural model E(Yx)=β0+β1x implies FYx must satisfy ∫yfYx(y)dy=β0+β1x. One way to specify the association between Z and Yx is by specifying a copula. When Z and Y are continuous, the variables UZ=FZ(Z) and UY=FYx(Yx) are both marginally Uniform(0,1) and their joint distribution FUZ,UY is called a copula (Aas et al., 2009). Rather than simulating data where X is randomized and then using rejection sampling to obtain observational data, a more computationally efficient method is as follows. Suppose Z and Y are continuous. Sample Z from FZ and then X from FX∣Z. Denote this sampled X value as x. Calculate UZ=FZ(Z) and sample a variable UY from FUY∣UZ, the conditional distribution of UY given UZ implied by the joint distribution FUZ,UY. This ensures UY is marginally Uniform(0,1). Hence, if we set Yx=FYx−1(UY), then Yx has marginal distribution FYx (thus satisfying the causal model) and Yx is correlated with Z (so there is confounding). By the consistency assumption that Y=Yx when X=x, this sampled Yx is the observed outcome Y. This method also works when Y is discrete. Seaman and Keogh’s (2023) proposal for simulating data for marginal structural survival models uses this approach, and it became apparent during the RSS Discussion Meeting that ED have also been using it. If Z is discrete, instead of setting UZ=FZ(Z), draw UZ∣Z∼Uniform(limz→Z−FZ(z),FZ(Z)). This ensures UZ is marginally Uniform(0,1). When Z is a random vector, ED suggest using vine copulas, which involves multiple (bivariate) copulas. Might it be easier to choose a scalar-valued function g and use a single (bivariate) copula to describe the association between random variables g(Z) and Yx? Seaman and Keogh (2023) call g a ‘risk score function’. Any function g could be chosen, allowing considerable flexibility in the choice of association. Evans and Didelez also provide a basis for maximum likelihood and Bayesian analysis of causal models. Example R1 with continuous L and Y serves to illustrate how this would work (or to reveal that I have misunderstood Section 5.1!) Here, Yab and La denote Y and L when we intervene to set A=a and B=b (note that Lab=La, because L is causally prior to B). We would specify models for FA, FLa, FB∣A,L, and FYab, with parameters θ, α, γ, and β, respectively (β are the parameters of interest). We could specify the association between La and Yab via a copula with parameter ρ, i.e. specify the joint distribution FUL,UY of UL=FLa(La) and UY=FYab(Yab) (ρ could depend on a and b). Figure 2 implies FYab∣A=a,L,B=FYab∣La. Hence, So, assuming (θ,γ) and (α,β,ρ) are variation independent, the likelihood for (α,β,ρ) is fUL,UY(FLa(l;α),FYab(y;β);ρ)fLa(l;α)fYab(y;β) and there is no need to specify models for FA or FB∣A,L. (If we wanted to simulate data from this model, we could sample a value a of A from FA, then L=La from FLa, then a value b of B from FB∣A=a,L, then calculate UL=FLa(L), sample UY from FUY∣UL, and calculate Y=Yab=FYab−1(UY)). Evans and Didelez write ‘Although we can fit models via maximum likelihood, […] double robust approaches may […] be more useful in practice’. In the Discussion Meeting, ED suggested a maximum likelihood analysis might be useful as a benchmark against which to compare statistical efficiency of another method, e.g. IPW. However, in Example R1, would not the efficiency of maximum likelihood depend on how flexible were the models for FLa and the association between La and Yab? Might the likelihood function be more useful for Bayesian analysis in contexts where prior information about β is available? I congratulate ED on a fascinating, thought-provoking article. I have learned much from it. As well as enabling me to simulate data from marginal structural survival models with fewer restrictions than previous simulation methods impose, it has helped me to simulate from the structural nested cumulative survival time model of Seaman et al. (2020). I enthusiastically propose the vote of thanks.