Estimating Distributions of Weighted Data

Question

I'm trying to build a bivariate copula-based model of income and wealth in Italy and I'm having trouble handling weighted data. I have access to micro data, a survey of about 10,000 households that includes the corresponding sample weights.

When calculating basic statistics (like mean and median) and even when performing linear regressions it is pretty easy to account for weights, besides there are useful packages for that (e. g. survey). But what do I do when I want to fit a parametric model of the distribution to weighted data? Or to estimate its kernel density?

I have a few ideas, but they seem to be pretty crude. For one, I could inflate my sample to the size of the universe. That is, I could multiply all weights by 100 (which would turn them into integers) and then create a vector that repeats each value of income and wealth a given number of times. But that would lead to a very large sample (which I believe still wouldn't be a perfect representation of the population) and will certainly put some extra strain on my computer.

I could also just round the weights off instead of multiplying them by 100, but this would still make the sample noticeably bigger and will inevitably skew the real proportions.

Another approach I came up with would be to normalize the weights (so that they sum up to one) and then randomly sample with repetitions from my initial sample with the corresponding vector of probability weights. R doesn't allow to draw the samples that are larger in size than the one that they are being drawn from. But I think that drawing the sample of the same size as the initial one will lead to some loss of information about the observed proportions. So I could draw the samples of the initial size as described above several times (how would I know how many is though?) and then combine them into one sample. And again, I will have a larger sample with some of the information lost along the way.

So I was wondering if there is a better way to handle weighted data. In some cases I think I could technically introduce the weights into the formula for computing the maximum likelihood for fitting a particular model, although I certainly wouldn't like to code that from the ground up. I will have to fit a lot of models as part of my project, both univariate (e. g. Singh-Mandala) for income and wealth and bivariate for copulas. I don't think the built in functions in any of the copula-related packages that I'm aware of allow one to account for weights. So any advice would help!

Answer 1

Actually estimating kernel density of weighted data is pretty straightforward. Recall that standard kernel density estimator is

$$ \hat{f_h}(x) = \frac{1}{n} \sum_{i=1}^n K_h(x-x_i) $$

where $K_h(x) = K(x/h)/h$ for kernel $K$ and bandwidth $h$, then it may be thought as weighted kernel density with all weights equal to $w_i = 1/n$,

$$ \hat{f_h}(x) = \sum_{i=1}^n w_i K_h(x-x_i) $$

this easily translates to the unequal weights scenario, where $\sum_{i=1}^n w_i = 1$. Moreover you don't have to look far for software estimating it as it is often available, e.g. R's density has weights parameter for such case.

As about fitting traditional distributions to weighted data, you can introduce weights in the likelihood function. Notice that if you re-sampled your data, so that some points would get duplicated, then you could obtain sample such as $x_1,x_1,x_1,x_2,x_2,\dots,x_k,x_k,x_k,x_k$, then the likelihood function becomes

$$ L(\theta\mid x_1,\dots,x_k) = f(x_1;\theta)\times f(x_1;\theta)\times f(x_1;\theta)\times f(x_2;\theta)\times f(x_2;\theta)\times \dots \\ \times f(x_k;\theta)\times f(x_k;\theta)\times f(x_k;\theta)\times f(x_k;\theta) = \\ f(x_1;\theta)^{3}\times f(x_2;\theta)^{2}\times\dots\times f(x_k;\theta)^{4} = \prod_i f(x_i;\theta)^{n_i} $$

where $f$ are density functions evaluated at $x_i$ points and parametrized by $\theta$. So if weights are non-negative and proportional to numbers of samples, then you can easily, by analogy, introduce weighted likelihood function of the same form

$$ L^w(\theta\mid x_1,\dots,x_k) = \prod_i f(x_i;\theta)^{w_i} $$

that can be maximized to find the optimal value of $\theta$. To lower the risk of numerical precision issues, it'd be better to work with it's logarithm,

$$ \log L^w(\theta\mid x_1,\dots,x_k) = \sum_i w_i \log f(x_i;\theta) $$

So you can use if for fitting the distribution by maximizing the weighted likelihood function.

Answer 2

I'm thinking through a similar problem right now - my initial thought is to generate a histogram using the weights (eg. numpy's histogram takes a weights vector) and then fit a parametric distribution to the resulting (normalized) frequency vector. Curve fitting is pretty low headache to implement. I'm happy to go into more detail if it would be helpful.