Use raking to impute multivariate distribution

Question

The Problem

I have a detailed dataset of 6 variabes, but for all but one year, I only have marginal distributions of 5 variables, the rest is missing. From that, I would like to obtain a full dataset.

To illustrate my problem imagine the case if I had a dataset of 3 variables (Age, Level, Year) and had only 2 margins for all but one year.

2011 (the “full year")

       Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7
1-20   2758  3531  2759  1472   963  1025   790   13298
21-40  2421  3530  2371  1583  1119  1376   720   13120
41-60  3497  3322  1818  1276   710   562   240   11425
61-80  4122  5554  2754  1786  1180   496   282   16174
80+    2740  5116  3050  2333  1889   656   334   16118
Sum   15538 21053 12752  8450  5861  4115  2366   70135

2010 (only 2 marginals)

       Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7
1-20                                              13197
21-40                                             12922
41-60                                             11369
61-80                                             16353
80+                                               15774
Sum   15151 21643 12611  8273  5586  4026  2325   69615

2009 (only 2 marginals)

       Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7
1-20                                              12684
21-40                                             12477
41-60                                             10844
61-80                                             16022
80+                                               14911
Sum   14367 20999 12201  8032  5235  3841  2263   66938

In this example, the challenge would be to fill in the missing values of the 2010 and 2009 data, using the 2011 distribution to make all necessary assumptions on multivariate marginal distributions.

Note that my problem is somewhat more complex. I have 5 dimensions (Age, Level, Region, Gender, Year) and three 4-dimensional marginal distributions tables (Gender, Level, Region, Year), (Age, Gender, Region, Year) and (Gender, Age, Level, Year) – but the problem is basically the same (I think … ?).

How can that be done?

Here’s my theory. Survey statistics has done a lot of reaserach on how to make any given (sample) data agree with population margins, I could use their methods. They do so by reweighting observations. Raking can be used to obtain a set of such weights. The rake function in the survey package allows raking with multidimensional joint distributions.

My idea was to use the full data as “sample”, and let the rake function compute weights that make the data agree with the marginal distributions.

Question 1: Is that a good idea? Do you have a better idea?

As you can imagine, there are several other problems to be solved. First one is that survey methods always assume that the “sample” data is individual data, I have counts, and no other characteristics of the individuals other than the stratification variables.

What I tried to do is specifying the counts as pre-raking survey weights. I tried that but the raking results in unplausible weights, that is, most weights are 0 and some are Inf.

Another possibility would be to duplicate every entry of the “population” data.frame depending on the number of observations, but that would bloat the data.frame significantly (it would have approximately 18 000 000 rows instead of currently 24 480 (the product of all dimension sizes) and cause memory problems.

If I specify no weights (all probabilities are 1), I get the following distribution of weights after raking:

> mean(lpg.rake$prob)
[1] 0.1289419
> median(lpg.rake$prob)
[1] 0.01700375
> quantile(lpg.rake$prob,c(.075,.1,.5,.9,.95,.99,.999))
        7.5%          10%          50%          90%          95%          99%        99.9% 
 0.001274434  0.001697273  0.017003749  0.155742479  0.322445467  1.597854446 11.711632910

Given that the mean of frequencies in the one year for which I have a full dataset is around 700, I would expect that to be the average of the weights.

Question 2: How should I deal with the fact that survey methods assume that observations are individuals, whereas I have one row for every combination of stratification variables plus a “count” variable?

A working example

Since my problem is very complex, and I am not familiar with the survey package, it is difficult for me to produce a minimal example that reproduces my problem. In some way, the complexity is the main problem. There is the example data in the survey package, you could take a look at that (link).

However, since my data is not secret, I don’t see why I could not give you the full dataset and what I have tried so far. Here it is.

The password for the zip file is: xB2KzbuK

Source: Ministry for Social Affairs (1995-2012): Pflegevorsorgebericht / Bericht des Arbeitskreises für Pflegevorsorge. Vienna, Austria AND own calculations.

Consequently, question 3 is: Is there an error in my code? What could I do in order to get meaningful results?

Answer 1

Generally speaking, you are on the right track. You do need to supply the frequencies in the full 5-way table as raking input weights, and then you should be able to produce the counts for a given year. 5-way raking isn't something very unusual: in most of the surveys that I weight, it's about five margins that I use, anyway (usually something like gender, age, race-ethnicity, education, geography), and in many surveys, there's more. So I am somewhat surprised that rake cannot fully handle this.

In your toy 2010 example, this is what I am getting with my own raking program written in Stata:

         Lvl1       Lvl2       Lvl3       Lvl4       Lvl5       Lvl6        Lvl7  
1-20   2686.67    3631.06    2731.71    1443.76     920.35    1005.43     778.02  
21-40  2341.20    3603.58    2330.44    1541.32    1061.65    1339.89     703.92  
41-60  3415.14    3424.75    1804.55    1254.68     680.27     552.66     236.96  
61-80  4079.18    5802.12    2770.08    1779.57    1145.66     494.26     282.14  
80+    2628.81    5181.50    2974.21    2253.68    1778.07     633.76     323.97  

If you are getting something like that in R, you are definitely on the right track syntax-wise, and if you don't get convergence with the full problem, then that's indeed a computational difficulty issue (some cells that are small in the full data year have to be adjusted a lot to take care of all the changes).

For reference, my Stata code is:

    clear       

    * get the raking package off the web
    net from http://www.stata-journal.com/software/sj14-1
    net install st0323 , force

    * toy data
    input age Count1  Count2  Count3  Count4  Count5  Count6  Count7
    1   2758  3531  2759  1472   963  1025   790   
    21  2421  3530  2371  1583  1119  1376   720   
    41  3497  3322  1818  1276   710   562   240   
    61  4122  5554  2754  1786  1180   496   282   
    80  2740  5116  3050  2333  1889   656   334   
    end

    * make it one record per line
    reshape long Count , i(age) j(level) 

    * minimal checks
    sum Count
    assert r(sum) == 70135

    gen byte _one = 1

    * define the margins
    mat age2010 = (13197, 12922, 11369, 16353, 15774)
    mat rownames age2010 = age
    mat colnames age2010 = 1 21 41 61 80
    mat coleq    age2010 = _one

    mat level2010 = (15151, 21643, 12611,  8273,  5586,  4026,  2325)
    mat rownames level2010 = level
    mat colnames level2010 = 1 2 3 4 5 6 7
    mat coleq    level2010 = _one

    * raking
    ipfraking [pw=Count], ctotal( age2010 level2010 ) gen( Count2010 )

    * results
    list age level Count Count2010, sepby(age)

Answer 2

"In this example, the challenge would be to fill in the missing values of the 2010 and 2009 data, using the 2011 distribution to make all necessary assumptions on multivariate marginal distributions."

Can you assume that in 2010 and 2009 the joint distribution is similar to that in 2011?

Not an answer to your question, but knowledge of the marginals give some information on individual cells: you can obtain upper and lower bounds for counts in each cell (Fréchet bounds). You might want to look up the book by Marcello D'Orazio et al. Statistical Matching, which as I recall deals with this topic.