R Under development (unstable) (2024-03-08 r86070 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(survey)
Loading required package: grid
Loading required package: Matrix
Loading required package: survival

Attaching package: 'survey'

The following object is masked from 'package:graphics':

    dotchart

> 
> data(api)
> dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> rclus1 <- as.svrepdesign(dclus1)
> 
> ## population marginal totals for each stratum
> pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
> pop.schwide <- data.frame(sch.wide=c("No","Yes"), Freq=c(1072,5122))
> 
> rclus1r <- rake(rclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
> 
> svymean(~api00, rclus1r)
        mean     SE
api00 641.23 26.873
> svytotal(~enroll, rclus1r)
         total     SE
enroll 3647300 463511
> 
> ff<-~stype+sch.wide
> poptotals<-colSums(model.matrix(ff,model.frame(ff,apipop)))
> rclus1g<-calibrate(rclus1, ~stype+sch.wide, poptotals,calfun="raking")
> 
> svymean(~api00,rclus1g)
        mean     SE
api00 641.23 26.874
> svytotal(~enroll,rclus1g)
         total     SE
enroll 3647280 463582
> 
> summary(weights(rclus1g)/weights(rclus1r))
       V1           V2          V3          V4           V5          V6   
 Min.   :1    Min.   :1   Min.   :1   Min.   :1    Min.   :1   Min.   :1  
 1st Qu.:1    1st Qu.:1   1st Qu.:1   1st Qu.:1    1st Qu.:1   1st Qu.:1  
 Median :1    Median :1   Median :1   Median :1    Median :1   Median :1  
 Mean   :1    Mean   :1   Mean   :1   Mean   :1    Mean   :1   Mean   :1  
 3rd Qu.:1    3rd Qu.:1   3rd Qu.:1   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1  
 Max.   :1    Max.   :1   Max.   :1   Max.   :1    Max.   :1   Max.   :1  
 NA's   :11   NA's   :4   NA's   :2   NA's   :13   NA's   :2   NA's   :4  
       V7          V8           V9         V10          V11    
 Min.   :1   Min.   :1    Min.   :1   Min.   :1    Min.   :1   
 1st Qu.:1   1st Qu.:1    1st Qu.:1   1st Qu.:1    1st Qu.:1   
 Median :1   Median :1    Median :1   Median :1    Median :1   
 Mean   :1   Mean   :1    Mean   :1   Mean   :1    Mean   :1   
 3rd Qu.:1   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1    3rd Qu.:1   
 Max.   :1   Max.   :1    Max.   :1   Max.   :1    Max.   :1   
 NA's   :4   NA's   :16   NA's   :9   NA's   :34   NA's   :21  
      V12              V13          V14         V15    
 Min.   :0.9997   Min.   :1    Min.   :1   Min.   :1   
 1st Qu.:1.0001   1st Qu.:1    1st Qu.:1   1st Qu.:1   
 Median :1.0001   Median :1    Median :1   Median :1   
 Mean   :1.0000   Mean   :1    Mean   :1   Mean   :1   
 3rd Qu.:1.0001   3rd Qu.:1    3rd Qu.:1   3rd Qu.:1   
 Max.   :1.0002   Max.   :1    Max.   :1   Max.   :1   
 NA's   :37       NA's   :13   NA's   :1   NA's   :12  
> 
> 
> ## Do it for a design without replicate weights
> dclus1r<-rake(dclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))
> 
> svymean(~api00, dclus1r)
        mean     SE
api00 641.23 23.704
> svytotal(~enroll, dclus1r)
         total     SE
enroll 3647300 400603
> 
> dclus1g<-calibrate(dclus1, ~stype+sch.wide, poptotals,calfun="raking")
> 
> svymean(~api00,dclus1g)
        mean     SE
api00 641.23 23.704
> svytotal(~enroll,dclus1g)
         total     SE
enroll 3647280 400603
> 
> summary(weights(dclus1g)/weights(dclus1r))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       1       1       1       1       1 
> 
> 
> 
> ## Example of raking with partial joint distributions
> pop.table <- xtabs(~stype+sch.wide,apipop)
> pop.imp<-data.frame(comp.imp=c("No","Yes"),Freq=c(1712,4482))
> dclus1r2<-rake(dclus1, list(~stype+sch.wide, ~comp.imp),
+                list(pop.table, pop.imp))
> svymean(~api00, dclus1r2)
        mean     SE
api00 642.62 22.732
> 
> ff1 <-~stype*sch.wide+comp.imp
> 
> poptotals1<-colSums(model.matrix(ff1,model.frame(ff1,apipop)))
> dclus1g2<-calibrate(dclus1, ~stype*sch.wide+comp.imp, poptotals1, calfun="raking")
> 
> svymean(~api00, dclus1g2)
        mean     SE
api00 642.61 22.731
> 
> summary(weights(dclus1r2)/weights(dclus1g2))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.999   1.000   1.000   1.000   1.000   1.002 
> 
> proc.time()
   user  system elapsed 
   1.48    0.15    1.59