Computing inverse probability weights -- conditional (multivariate) density estimation?

Question

The general version:

I need to estimate $f(A | X)$ where $A$ and $X$ are continuous and multivariate. I'd rather do it nonparametrically because I don't have a good functional form in mind and $\hat{f}(A | X)$ needs to be something like unbiased. I wanted to use a conditional kernel density estimator, but I realized I would need to quantize $X$ first. Then I had the idea to estimate $\hat{f}(A , X)$ and $\hat{f}(X)$ from the data and use that to compute $\hat{f}(A | X)$, or maybe I read it somewhere and don't remember where.

Is there a reason this procedure wouldn't be valid? Is there a better or more honest approach than kernel density? Also, is there an issue with estimating population densities from sample densities nonparametrically? The data is survey data, and I have survey weights; should I incorporate them somehow?

---

The case-specific version:

It's probably worth mentioning that I'm going to use these estimates for inverse probability-of-treatment weights in a marginal structural model, as in Robins (2000) (ungated PDF). I observe a sequence of "treatments" $\{a_t\}_{t=0}^{4}$ and a sequence of time-varying confounders $\{x_t\}_{t=0}^{4}$ with respect to some outcome $\tilde{y}$ that occurs at $t=T+1$. I've hypothesized a simple parametric causal relationship, $E[\tilde{Y} | \vec{a}]=\beta'\vec{a}$, but since there's a time-varying confounder $\beta$ is a biased estimate of the "average treatment effect," and the confounder can't be added as a regressor because it's on the causal path and that will also bias $\beta$. Fortunately Doc Robins figured out that I can get unbiased/unconfounded and reasonably efficient estimates if I re-weight my observations by: $$ w_i = \prod_{s=0}^{4} \frac{ f(a_s | a_{s<t}) }{ f(a_s | a_{s<t},x_{s<t}) } $$

My question: That sequence of weights is really what I need an estimate for. Robins recommends logistic regression. But $a_t$ lies in $[0,\infty)^7$, is measured on $\{0,\dots\}^{7}$, and for all practical purposes lies in a finite subset thereof. $x_t$ lies in a closed interval, but only because it's really the mean of a few variables, each measured on a finite subset of $\{0,\dots,12\}$.

So I had a few ideas:

1. Estimate $f(a_t, a_{s<t}, x_{s<t})$ and $f(x,a_{s<t})$ nonparametrically
2. Estimate $f(a_t | a_{s<t}, x_{s<t})$ with beta regression and $f( x_{s<t}, a_{s<t})$ nonparametrically
3. Estimate $f(x_{t-1}|a_t,a_{s<t},x_{s<(t-1)})$ with beta regression, and estimate $f(a_t, a_{s<t},x_{s<(t-1)})$) by "chaining" beta regressions back through time to express the whole thing as a conditional.
4. Something actually coherent and honest in propagating uncertainty, that I obviously haven't thought of.
5. Bayes? I do know Stan and JAGS, but MCMC will probably explode my computer (I don't want to deal with EC2).

I haven't found any hints in the literature, since multivariate treatments are rare in causal modeling. What should I do?

Bonus points: how do you feel about the notation $a_{s<t}$ to represent $\{a_s\}_{s=0}^{t}$ instead of something like $\vec{a}_{t-1}$?

Answer 1

The basic idea

As per Chen, Linton, and Robinson (2001), the "default" technique for conditional univariate kernel density estimation is to find, for bandwidths $a,b,c$, $$\frac{\hat{f}_{ab}(y,z)}{\hat{f}_c(z)}=\hat{f}_{abc}(y|z)$$

Then, with numerator bandwidth $(a,b)$ and denominator bandwidth $c$ and $a=b=c$, the following central limit result holds under certain independence and consistency assumptions (which are only really restrictive when $y=x_t,z=x_{t-1}$): $$ \sqrt{na^2}\left(\hat{f}_{abc=aaa}(y|z)-f(y|z)\right)\xrightarrow{d}N(0,V) $$

where $$\begin{align} \hat{V}&=\left(\int K(u)^2du\right)^2\cdot\frac{\hat{f}_{aaa}(y|z)}{\hat{f}_a(z)}\\&=\left(\int K(u)^2du\right)^2\cdot\hat{f}_{aa }(y,z) \end{align}$$

Although I've never see a frequentist weighted model (even intro-stats WLS) make an attempt to account for the variance of the estimated weights. For now I'm going to follow that convention but if I get results here I'll see if I can work it into a fully Bayesian model that will propagate uncertainty more honestly. So yes, estimating the conditional density by estimating the joint and marginal densities is standard procedure.

Applicability to my case

It's not explicitly clear from that paper is how this generalizes to the case when when $y=x_t$ and $z=\left(x_s\right)_{s=1}^{t-1}$, and $x_s=\left(\begin{smallmatrix} x_{s,1}\\ \ddots \\x_{s,D} \end{smallmatrix}\right)$. But I guess this is really just the same thing as one big long sequence $x=\left(\left(x_{s,d}\right)_{d=1}^{D}\right)_{s=1}^{t-1}$ which seems perfectly manageable according to Robinson (1983) (cited in Chen, et al). Again, using Bayes' rule to estimate the conditional density seems perfectly acceptable.

Bandwidth

The final issue is bandwidth selection. Bandwidth is now a block matrix of the form $$ B=\left(\begin{matrix} B^{numerator}&0\\0&B^{denominator} \end{matrix}\right)=\left(\begin{matrix} \left(\begin{matrix}a_{1,1}&&B^{num}_1\\&\ddots&\\B^{num}_2&&a_{t,D}\end{matrix}\right)&0\\0&\left(\begin{matrix}c_{1,1}&&B^{denom}_1\\&\ddots&\\B^{denom}_2&&c_{t-1,D}\end{matrix}\right) \end{matrix}\right) $$

which is a mess. When bandwidth $H=hH_0$ such that $|H_0|=1$, then $b^*\sim\sqrt[4+D]{N}$, but this result would apply separately to $B^{num}$ and $B^{denom}$ rather than to $B$ as a whole (source, someone's lecture notes).

Chen et al find an optimal bandwidth $a=b=c$ (in their 2-d case) for a given level of $z$ that looks like it generalizes to the case when $y$ and $z$ are multivariate. They suggest setting $z=\mu$ where $\mu$ is the theoretical mean that would be induced under joint normality, and they derive $\hat{a}(\mu)$.

A more general version of the same result is in another section of those lecture notes, called "rule-of-thumb" bandwidth. They also derive optimal bandwidth as a function of a general cross-validation procedure.

Computation

I have a 7-dimensional treatment over 3 time periods, so I have up to a 21-dimensional density to estimate. And I forgot about baseline covariates. I have something like 30 baseline covariates, so I would end up trying to estimate a 51-dimensional distribution, a 44-dimensional distribution, and a 37-dimensional distribution. And that's not to mention that the extreme dimensionality will require an impossibly large sample. Scott & Wand (1991) report that a sample size of 50 in one dimension is equivalent to well over 1 million in 8 dimensions... no mention of 30. No amount of these can express how I feel right now.

Conclusion

So I just wasted a week of my life on this. Oh well. Instead, I'm going to use MCMC to fit parametric treatment and outcome models simultaneously, so that the IPT weights end up being a function of the posterior predictive densities from the treatment model. Then I'll step through linear, quadratic, and cubic forms for the treatment model and see which one fits the best.