1
$\begingroup$

Is there a way to calibrate against multidimensional control totals using the survey package? I have been able to estimate poststratification weights, but now I'm stuck with the calibration. (I have a clustered design and need the weights to stay the same within each cluster.) Do I have to recode all combinations into a new attribute?

> df <- data.frame(id=c(1,1,1,2,2,3,4,5,5),
+                  h=factor(c(1,1,1,0,0,1,0,1,1)),
+                  p=factor(c(1,0,1,0,1,0,0,0,1)))
> 
> df
  id h p
1  1 1 1
2  1 1 0
3  1 1 1
4  2 0 0
5  2 0 1
6  3 1 0
7  4 0 0
8  5 1 0
9  5 1 1
> des <- svydesign(id=~id, data=df)
Warning message:
In svydesign.default(id = ~id, data = df) :
  No weights or probabilities supplied, assuming equal probability
> pop <- xtabs(Freq~., data.frame(h=factor(c(0, 1, 0, 1)),
+                                 p=factor(c(0, 0, 1, 1)),
+                                 Freq=c(40, 35, 25, 10)))
> pop
   p
h    0  1
  0 40 25
  1 35 10
> 
> ps <- postStratify(des, strata=~h*p, population=pop)
> weights(ps)
[1]  3.333333 11.666667  3.333333 20.000000 25.000000 
[6]  11.666667 20.000000 11.666667  3.333333

The documentation and Lumley's "survey book" both don't show how to do this and if it's possible at all.

$\endgroup$
2
$\begingroup$

The documentation contains a small hint on how the names of the population vector for calibrate() should be formed: They are compared with the results of model.matrix(). The question remains how to derive a proper population vector from a table of totals like the one passed to postStratify().

I use a special call to model.matrix() to convert a totals table to a population vector. The main idea is to compute a model.matrix() from the totals table plus a Freq column, to scale the matrix contents by the Freq column and then to sum up column-wise. This results exactly in the population vector. Please see the complete example below for details.

> library(survey)

Attaching package: ‘survey’

The following object(s) are masked from ‘package:graphics’:

    dotchart

> df <- data.frame(id=c(1,1,1,2,2,3,4,5,5,6,6,6,7,7,7)
+                  , h=factor(c(1,1,1,0,0,1,0,1,1,0,0,0,0,0,0))
+                  , p=factor(c(1,0,1,0,1,0,0,0,1,1,0,1,1,0,0))
+                  , q=factor(c(0,0,1,1,1,0,1,1,1,0,0,1,1,1,0))
+                  )
> df$hp <- do.call(interaction, df[, c('h', 'p')])
    > df$hq <- do.call(interaction, df[, c('h', 'q')])
> df$pq <- do.call(interaction, df[, c('p', 'q')])
    > df$hpq <- do.call(interaction, df[, c('h', 'p', 'q')])
> df
   id h p q  hp  hq  pq   hpq
1   1 1 1 0 1.1 1.0 1.0 1.1.0
2   1 1 0 0 1.0 1.0 0.0 1.0.0
3   1 1 1 1 1.1 1.1 1.1 1.1.1
4   2 0 0 1 0.0 0.1 0.1 0.0.1
5   2 0 1 1 0.1 0.1 1.1 0.1.1
6   3 1 0 0 1.0 1.0 0.0 1.0.0
7   4 0 0 1 0.0 0.1 0.1 0.0.1
8   5 1 0 1 1.0 1.1 0.1 1.0.1
9   5 1 1 1 1.1 1.1 1.1 1.1.1
10  6 0 1 0 0.1 0.0 1.0 0.1.0
11  6 0 0 0 0.0 0.0 0.0 0.0.0
12  6 0 1 1 0.1 0.1 1.1 0.1.1
13  7 0 1 1 0.1 0.1 1.1 0.1.1
14  7 0 0 1 0.0 0.1 0.1 0.0.1
15  7 0 0 0 0.0 0.0 0.0 0.0.0
> 
> des <- svydesign(id=~id, data=df)
Warning message:
In svydesign.default(id = ~id, data = df) :
  No weights or probabilities supplied, assuming equal probability
> 
> pop <- data.frame(h=factor(c(0, 1, 0, 1, 0, 1, 0, 1)),
+                   p=factor(c(0, 0, 1, 1, 0, 0, 1, 1)),
+                   q=factor(c(0, 0, 0, 0, 1, 1, 1, 1)),
+                   Freq=c(40, 35, 25, 35, 45, 30, 30, 30))
> pop$hp <- do.call(interaction, pop[, c('h', 'p')])
    > pop$hq <- do.call(interaction, pop[, c('h', 'q')])
> pop$pq <- do.call(interaction, pop[, c('p', 'q')])
    > pop$hpq <- do.call(interaction, pop[, c('h', 'p', 'q')])
> pop
  h p q Freq  hp  hq  pq   hpq
1 0 0 0   40 0.0 0.0 0.0 0.0.0
2 1 0 0   35 1.0 1.0 0.0 1.0.0
3 0 1 0   25 0.1 0.0 1.0 0.1.0
4 1 1 0   35 1.1 1.0 1.0 1.1.0
5 0 0 1   45 0.0 0.1 0.1 0.0.1
6 1 0 1   30 1.0 1.1 0.1 1.0.1
7 0 1 1   30 0.1 0.1 1.1 0.1.1
8 1 1 1   30 1.1 1.1 1.1 1.1.1
> 
> #fm <- '~hp+hq+pq'
> fm <- '~h:p+h:q+p:q'
> # Both options above seem to produce the same results
> 
> (mm <- as.data.frame(model.matrix(as.formula(paste0(fm, '+Freq')), pop)))
  (Intercept) Freq h0:p0 h1:p0 h0:p1 h1:p1 h0:q1 h1:q1 p1:q1
1           1   40     1     0     0     0     0     0     0
2           1   35     0     1     0     0     0     0     0
3           1   25     0     0     1     0     0     0     0
4           1   35     0     0     0     1     0     0     0
5           1   45     1     0     0     0     1     0     0
6           1   30     0     1     0     0     0     1     0
7           1   30     0     0     1     0     1     0     1
8           1   30     0     0     0     1     0     1     1
> mm <- mm * mm$Freq
    > mm$Freq <- NULL
> (pop <- unlist(lapply(mm, sum)))
(Intercept)       h0:p0       h1:p0       h0:p1       h1:p1       h0:q1 
        270          85          65          55          65          75 
      h1:q1       p1:q1 
         60          60 
> 
> cal <- NULL
> cal <- calibrate(des, formula=as.formula(fm), population=pop,
+                  calfun="raking")
Loading required package: MASS
> weights(cal)
        1         2         3         4         5         6         7         8 
31.902369 19.048815 16.548815 16.032544  8.967456 19.048815 16.032544 26.902369 
        9        10        11        12        13        14        15 
16.548815 28.097631 18.451185  8.967456  8.967456 16.032544 18.451185 
attr(,"eta")
           [,1]
[1,]  2.5320998
[2,]  0.3830288
[3,]  0.4149051
[4,]  0.8035854
[5,]  0.9305805
[6,] -0.1405080
[7,]  0.3452095
[8,] -1.0015752
> svytotal(~h, cal)
   total SE
h0   140  0
h1   130  0
> svytotal(~p, cal)
   total SE
p0   150  0
p1   120  0
> svytotal(~q, cal)
   total SE
q0   135  0
q1   135  0
> svytotal(~hp, cal)
      total SE
hp0.0    85  0
hp1.0    65  0
hp0.1    55  0
hp1.1    65  0
> svytotal(~hq, cal)
      total SE
hq0.0    65  0
hq1.0    70  0
hq0.1    75  0
hq1.1    60  0
> svytotal(~pq, cal)
      total SE
pq0.0    75  0
pq1.0    60  0
pq0.1    75  0
pq1.1    60  0
> svytotal(~hpq, cal) 
          total     SE
hpq0.0.0 36.902 6.4914
hpq1.0.0 38.098 6.4914
hpq0.1.0 28.098 6.4914
hpq1.1.0 31.902 6.4914
hpq0.0.1 48.098 6.4914
hpq1.0.1 26.902 6.4914
hpq0.1.1 26.902 6.4914
hpq1.1.1 33.098 6.4914

(Just for the record: The R function interaction() recodes all combinations of given factors into one single factor.)

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.