# Computing inverse probability weights — conditional (multivariate) density estimation?

## 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}$?

## 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.

• "So I just wasted a week of my life on this." That's called learning and research. Being a master student you should embrace that because there is more like this coming. There usually are no shortcuts in research because often no-one knows the way! – Momo Jun 26 '14 at 20:39