Calibrate against multi-dimensional population totals using the survey package?

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.

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")
> 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.)