Thinking out loud about missing data

[Ben Goodrich]

I've been reading about so-called "k Nearest Neighbor Density Estimators". This short paper

http://scholar.google.com/scholar?cluster=16043802846384639781&hl=en&as_sdt=0,33

compares three of these things. Basically, they are scalar functions of a vector --- based on Euclidean distances to other rows --- that are not exactly proper density functions (because they integrate to infinity) but converge pointwise to a genuine density function as the number of observations increases. In particular, the function being pushed in the above paper converges to a gamma density. So, someone could potentially utilize this function if it were implemented in Stan and live with the fact that it is not quite a proper density or utilize a gamma density and live with the fact that it is not quite the right function. My gut says that it wouldn't matter much either way relative to the magnitude of the problem I am contemplating, which is that the datamatrix could have an arbitrary missingness pattern.

If all the variables were continuous but partially missing, then a full Bayesian approach would be pretty straightforward: Make parameters for the missing cells, put priors on them, index-fiddle to put these parameters into the right locations in the datamatrix, use something from that article as your "likelihood", let Stan run until convergence, and get a whole bunch of completed datamatrices to use in a subsequent analysis. Basically, it is just a joint multiple imputation algorithm with relaxed parametric assumptions.

Finding the k nearest neighbors is computationally difficult. But among rows that are fully observed, you only have to calculate the Euclidean distances once. And among rows that are incompletely observed, you can calculate a partial sum of squared differences once based on the columns where both rows are observed once and then you just have to tack on the remaining sum of squared differences for the columns where one or both rows are filled with parameters. And even still, we could do the whole thing in doubles to find the k nearest neighbors and then handle the derivatives manually. In other words, we don't need to mess with the derivatives for the (k + m)th distances because they never enter the "likelihood".

If the datamatrix includes discrete variables with missingness, we are not completely stuck by the fact that Stan can't do discrete parameters yet. Explode discrete variables into a full set of dummy variables and if missing, define an underlying parameter as a simplex of probabilities with some reasonable prior. The contribution to the total sum of squared differences for those columns is just 2.0 times the probability that the discrete variable differs between the two rows. To see this, say you have a variable like education with four levels (less than high school, HS or GED, some college, college degree). Then you might have two fully observed people that were both category two

row i: [...,0,1,0,0,...]
row j: [...,0,1,0,0,...]

in which case the contribution to the sum of squared differences is zero. Or two observed people with different levels of eduction

row i: [...,0,1,0,0,...]
row j: [...,0,0,0,1,...]

in which case the contribution to the sum of squared differences is 2. So, if row j were instead missing, the contribution to the sum of squares would be 2 times the probability that j is not in the HS / GED category. If both i and j were missing, then the contribution to the sum of squares is more complicated to calculate but has the same basic form. When Stan is done, you can draw realizations from a multinomial distribution with the corresponding posterior probability vector as your imputations.

Of course, that assumes it is reasonable to say that the "distance" between i and j is the same and has a weird cardinal value no matter how i and j differ in their education levels. E.g. if i drops out of high school and j graduates from college, they have the same distance as if i graduates from high school and j drops out of college. But at least the contribution to the distance is zero for two people with the same education level, and the fact that you have missing data on discrete variables is a much bigger problem.

I can't think of anything clever to deal with variables that are counts, but it is probably fine to treat them as continuous, impute missing elements with real parameters, and perhaps round the imputed values to the nearest integer once Stan is finished.

Thoughts?
Ben

[Ben Goodrich]

On Friday, February 8, 2013 2:38:50 AM UTC-5, Ben Goodrich wrote:

I've been reading about so-called "k Nearest Neighbor Density Estimators". This short paper

http://scholar.google.com/scholar?cluster=16043802846384639781&hl=en&as_sdt=0,33

compares three of these things. Basically, they are scalar functions of a vector --- based on Euclidean distances to other rows --- that are not exactly proper density functions (because they integrate to infinity) but converge pointwise to a genuine density function as the number of observations increases.

Or ... is this potentially a convergence diagnostic (and a steak knife, etc.) My thinking here is that as the number of iterations increases a Markov chain converges to its stationary distribution and the non-parametric density estimator converges to the true density. Since we have the true log-posterior (ignoring normalizing constants) from Stan and can approximate the log-posterior non-parametrically, then the difference between them should become constant (i.e. zero if Stan didn't ignore normalizing constants).

This would at least make sense for a Gibbs sampler that moves every iteration. For a Markov chain with a Metropolis acceptance test, then sometimes the chain won't move from one iteration to the next. I presume one should exclude those when searching for the k-th nearest neighbors, but that may well affect the asymptotics and / or the optimal choice of k.

Here is an R function for the non-parametric density estimator pushed in the above paper:

kNN_DE <- function(x, k = floor(2 * nrow(x)^(1/ncol(x)))) {

  dist_x <- as.matrix(dist(x))

  dist_x[dist_x <= .Machine$double.eps] <- NA_real_

  N <- nrow(x)

  half_D <- ncol(x) / 2

  c_D <- pi^(half_D) / gamma(half_D + 1)

  P <- k / (N * c_D * apply(dist_x, 2, FUN = function(y) sort(y)[k])^ncol(x))

  return(P)

}

and a test with 8 schools

par(mfrow = c(2,1))

library(rstan)

fit <- stan_demo(model = "eight_schools", chains = 1, warmup = 0, iter = 4000, 

                 seed = 123456, init = 0)

mat <- as.matrix(fit)

lp <- log(kNN_DE(mat[,-ncol(mat)]))

smoothScatter(1:nrow(mat), mat[,"lp__"] - lp, main = "Not Converged", las = 1,

              ylab = "Difference", xlab = "Post-warmup iteration")

fit <- stan_demo(model = "eight_schools", chains = 1, warmup = 1000, iter = 5000, 

                 seed = 123456, init = 0)

mat <- as.matrix(fit)

lp <- log(kNN_DE(mat[,-ncol(mat)]))

smoothScatter(1:nrow(mat), mat[,"lp__"] - lp, main = "Converged", las = 1,

              ylab = "Difference", xlab = "Post-warmup iteration")

which yields

That at least looks pretty compelling. One limitation of non-parametric density estimators is that they do a worse job in the tails than near the mode, which may be what is happening at the iterations just before 3000.

Thoughts?

Ben

[Bob Carpenter]

Is the idea that the difference should go to zero
in the limit as you get better non-parametric density
estimates?

How do you monitor convergence of the density
estimator?

- Bob

> <https://lh3.googleusercontent.com/-ovITzBIC1ts/USm6gGyVF9I/AAAAAAAAAHk/wtlST1pReLk/s1600/kNN_DE.png>

>
> That at least looks pretty compelling. One limitation of non-parametric density estimators is that they do a worse job
> in the tails than near the mode, which may be what is happening at the iterations just before 3000.
>
> Thoughts?
> Ben
>

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Ben Goodrich]

On Sunday, February 24, 2013 3:52:21 AM UTC-5, Bob Carpenter wrote:

Is the idea that the difference should go to zero
in the limit as you get better non-parametric density
estimates?

I think it would be zero asymptotically if Stan's lp__ reflected the normalizing constants. So, the difference will be non-zero but shouldn't be changing much as the number of draws from the stationary distribution increases.

 

How do you monitor convergence of the density
estimator?

I don't think anyone has derived great theoretical results. You can draw a finite number of times from a known joint distribution like multivariate normal and then see how long it takes for the non-parametric estimator to get close. 

Ben

[Bob Carpenter]

I'd think in high dimensions, it'd be terrible, because
no points are close to each other on average, so kNN tends
to break down. At that point, what I see people in the
machine learning literature do is try to find submanifolds
of lower dimension on which the relevant values lie.

But we're not really doing multivariate convergence analysis,
just a bunch of marginals. (I believe Brooks and Gelman extended
R-hat to a multivariate setting.)

Stan's lp__ should be OK up to a proportion -- I wonder how
much of a problem we get by introducing an arbitrary rescaling?

- Bob

[Ben Goodrich]

On Sunday, February 24, 2013 5:11:03 PM UTC-5, Bob Carpenter wrote:

I'd think in high dimensions, it'd be terrible, because
no points are close to each other on average, so kNN tends
to break down.  

That is possible. The authors do recommend that k be a function of the dimensionality, namely 2 * N^(1/d). I did it with 10000 draws from a multivariate normal distribution with 24 highly correlated variables and it looked okay. I haven't been able to do it with a truly large number of dimensions yet because I haven't optimized the RAM enough to do much more than a 10000x10000 matrix of distances.

 

At that point, what I see people in the
machine learning literature do is try to find submanifolds
of lower dimension on which the relevant values lie.

Possible. I don't know offhand how that would affect the convergence to the true density.

 

But we're not really doing multivariate convergence analysis,
just a bunch of marginals.  (I believe Brooks and Gelman extended
R-hat to a multivariate setting.)

That is what I think is a problem with looking at convergence marginally. But I believe Gelman does not like the Brooks and Gelman statistic.

 

Stan's lp__ should be OK up to a proportion -- I wonder how
much of a problem we get by introducing an arbitrary rescaling?

I think it is not a problem at all, as long as you are just expecting the difference between the parametric and non-parametric to be a constant, rather than constantly zero.

Ben

[Andrew Gelman]

The challenge is that lp__ gives us the multivariate density, but density estimation works well in one dimension. One approach is to monitor one dimension at a time, integrating out all the others. This can actually be done using path sampling. Fan, Brooks, and I proposed such a convergence diagnostic in Section 5 of this paper:
http://www.stat.columbia.edu/~gelman/research/published/fanbrooksgelman.pdf
It was a cool idea but didn't seem so practical at the time. With Stan, though, it could be done automatically, so maybe it is actually worth a try.
Andrew