1
2R version 3.1.0 (2014-04-10) -- "Spring Dance"
3Copyright (C) 2014 The R Foundation for Statistical Computing
4Platform: x86_64-apple-darwin13.1.0 (64-bit)
5
6R is free software and comes with ABSOLUTELY NO WARRANTY.
7You are welcome to redistribute it under certain conditions.
8Type 'license()' or 'licence()' for distribution details.
9
10R is a collaborative project with many contributors.
11Type 'contributors()' for more information and
12'citation()' on how to cite R or R packages in publications.
13
14Type 'demo()' for some demos, 'help()' for on-line help, or
15'help.start()' for an HTML browser interface to help.
16Type 'q()' to quit R.
17
18> library(survey)
19
20Attaching package: 'survey'
21
22The following object is masked from 'package:graphics':
23
24    dotchart
25
26>
27> data(api)
28> dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
29> rclus1 <- as.svrepdesign(dclus1)
30>
31> ## population marginal totals for each stratum
32> pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
33> pop.schwide <- data.frame(sch.wide=c("No","Yes"), Freq=c(1072,5122))
34>
35> rclus1r <- rake(rclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
36>
37> svymean(~api00, rclus1r)
38        mean     SE
39api00 641.23 26.873
40> svytotal(~enroll, rclus1r)
41         total     SE
42enroll 3647300 463511
43>
44> ff<-~stype+sch.wide
45> poptotals<-colSums(model.matrix(ff,model.frame(ff,apipop)))
46> rclus1g<-calibrate(rclus1, ~stype+sch.wide, poptotals,calfun="raking")
47Loading required package: MASS
48>
49> svymean(~api00,rclus1g)
50        mean     SE
51api00 641.23 26.874
52> svytotal(~enroll,rclus1g)
53         total     SE
54enroll 3647280 463582
55>
56> summary(weights(rclus1g)/weights(rclus1r))
57       V1           V2          V3          V4           V5          V6
58 Min.   :1    Min.   :1   Min.   :1   Min.   :1    Min.   :1   Min.   :1
59 1st Qu.:1    1st Qu.:1   1st Qu.:1   1st Qu.:1    1st Qu.:1   1st Qu.:1
60 Median :1    Median :1   Median :1   Median :1    Median :1   Median :1
61 Mean   :1    Mean   :1   Mean   :1   Mean   :1    Mean   :1   Mean   :1
62 3rd Qu.:1    3rd Qu.:1   3rd Qu.:1   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1
63 Max.   :1    Max.   :1   Max.   :1   Max.   :1    Max.   :1   Max.   :1
64 NA's   :11   NA's   :4   NA's   :2   NA's   :13   NA's   :2   NA's   :4
65       V7          V8           V9         V10          V11
66 Min.   :1   Min.   :1    Min.   :1   Min.   :1    Min.   :1
67 1st Qu.:1   1st Qu.:1    1st Qu.:1   1st Qu.:1    1st Qu.:1
68 Median :1   Median :1    Median :1   Median :1    Median :1
69 Mean   :1   Mean   :1    Mean   :1   Mean   :1    Mean   :1
70 3rd Qu.:1   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1    3rd Qu.:1
71 Max.   :1   Max.   :1    Max.   :1   Max.   :1    Max.   :1
72 NA's   :4   NA's   :16   NA's   :9   NA's   :34   NA's   :21
73      V12              V13          V14         V15
74 Min.   :0.9997   Min.   :1    Min.   :1   Min.   :1
75 1st Qu.:1.0001   1st Qu.:1    1st Qu.:1   1st Qu.:1
76 Median :1.0001   Median :1    Median :1   Median :1
77 Mean   :1.0000   Mean   :1    Mean   :1   Mean   :1
78 3rd Qu.:1.0001   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1
79 Max.   :1.0002   Max.   :1    Max.   :1   Max.   :1
80 NA's   :37       NA's   :13   NA's   :1   NA's   :12
81>
82>
83> ## Do it for a design without replicate weights
84> dclus1r<-rake(dclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
85>
86> svymean(~api00, dclus1r)
87        mean     SE
88api00 641.23 23.704
89> svytotal(~enroll, dclus1r)
90         total     SE
91enroll 3647300 400603
92>
93> dclus1g<-calibrate(dclus1, ~stype+sch.wide, poptotals,calfun="raking")
94>
95> svymean(~api00,dclus1g)
96        mean     SE
97api00 641.23 23.704
98> svytotal(~enroll,dclus1g)
99         total     SE
100enroll 3647280 400603
101>
102> summary(weights(dclus1g)/weights(dclus1r))
103   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
104      1       1       1       1       1       1
105>
106>
107>
108> ## Example of raking with partial joint distributions
109> pop.table <- xtabs(~stype+sch.wide,apipop)
110> pop.imp<-data.frame(comp.imp=c("No","Yes"),Freq=c(1712,4482))
111> dclus1r2<-rake(dclus1, list(~stype+sch.wide, ~comp.imp),
112+                list(pop.table, pop.imp))
113> svymean(~api00, dclus1r2)
114        mean     SE
115api00 642.62 22.732
116>
117> ff1 <-~stype*sch.wide+comp.imp
118>
119> poptotals1<-colSums(model.matrix(ff1,model.frame(ff1,apipop)))
120> dclus1g2<-calibrate(dclus1, ~stype*sch.wide+comp.imp, poptotals1, calfun="raking")
121>
122> svymean(~api00, dclus1g2)
123        mean     SE
124api00 642.61 22.731
125>
126> summary(weights(dclus1r2)/weights(dclus1g2))
127   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
128  0.999   1.000   1.000   1.000   1.000   1.002
129>
130> proc.time()
131   user  system elapsed
132  0.459   0.032   0.499
133