1library(survey)
2
3data(api)
4dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
5rclus1 <- as.svrepdesign(dclus1)
6
7## population marginal totals for each stratum
8pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
9pop.schwide <- data.frame(sch.wide=c("No","Yes"), Freq=c(1072,5122))
10
11rclus1r <- rake(rclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
12
13svymean(~api00, rclus1r)
14svytotal(~enroll, rclus1r)
15
16ff<-~stype+sch.wide
17poptotals<-colSums(model.matrix(ff,model.frame(ff,apipop)))
18rclus1g<-calibrate(rclus1, ~stype+sch.wide, poptotals,calfun="raking")
19
20svymean(~api00,rclus1g)
21svytotal(~enroll,rclus1g)
22
23summary(weights(rclus1g)/weights(rclus1r))
24
25
26## Do it for a design without replicate weights
27dclus1r<-rake(dclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
28
29svymean(~api00, dclus1r)
30svytotal(~enroll, dclus1r)
31
32dclus1g<-calibrate(dclus1, ~stype+sch.wide, poptotals,calfun="raking")
33
34svymean(~api00,dclus1g)
35svytotal(~enroll,dclus1g)
36
37summary(weights(dclus1g)/weights(dclus1r))
38
39
40
41## Example of raking with partial joint distributions
42pop.table <- xtabs(~stype+sch.wide,apipop)
43pop.imp<-data.frame(comp.imp=c("No","Yes"),Freq=c(1712,4482))
44dclus1r2<-rake(dclus1, list(~stype+sch.wide, ~comp.imp),
45               list(pop.table, pop.imp))
46svymean(~api00, dclus1r2)
47
48ff1 <-~stype*sch.wide+comp.imp
49
50poptotals1<-colSums(model.matrix(ff1,model.frame(ff1,apipop)))
51dclus1g2<-calibrate(dclus1, ~stype*sch.wide+comp.imp, poptotals1, calfun="raking")
52
53svymean(~api00, dclus1g2)
54
55summary(weights(dclus1r2)/weights(dclus1g2))
56