R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

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 a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("survey-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('survey')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "SE"
> 
> ### * SE
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SE
> ### Title: Extract standard errors
> ### Aliases: SE SE.default SE.svrepstat
> ### Keywords: models
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "api"
> 
> ### * api
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: api
> ### Title: Student performance in California schools
> ### Aliases: api apipop apiclus1 apiclus2 apistrat
> ### Keywords: datasets
> 
> ### ** Examples
> 
> library(survey)
> data(api)
> mean(apipop$api00)
[1] 664.7126
> sum(apipop$enroll, na.rm=TRUE)
[1] 3811472
> 
> #stratified sample
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> summary(dstrat)
Stratified Independent Sampling design
svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02262 0.02262 0.03587 0.04014 0.05339 0.06623 
Stratum Sizes: 
             E  H  M
obs        100 50 50
design.PSU 100 50 50
actual.PSU 100 50 50
Population stratum sizes (PSUs): 
   E    M    H 
4421 1018  755 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc"     
> svymean(~api00, dstrat)
        mean     SE
api00 662.29 9.4089
> svytotal(~enroll, dstrat, na.rm=TRUE)
         total     SE
enroll 3687178 114642
> 
> # one-stage cluster sample
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> summary(dclus1)
1 - level Cluster Sampling design
With (15) clusters.
svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> svymean(~api00, dclus1)
        mean     SE
api00 644.17 23.542
> svytotal(~enroll, dclus1, na.rm=TRUE)
         total     SE
enroll 3404940 932235
> 
> # two-stage cluster sample
> dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
> summary(dclus2)
2 - level Cluster Sampling design
With (40, 126) clusters.
svydesign(id = ~dnum + snum, fpc = ~fpc1 + fpc2, data = apiclus2)
Probabilities:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.003669 0.037740 0.052840 0.042390 0.052840 0.052840 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc1"     "fpc2"    
> svymean(~api00, dclus2)
        mean     SE
api00 670.81 30.099
> svytotal(~enroll, dclus2, na.rm=TRUE)
         total     SE
enroll 2639273 799638
> 
> # two-stage `with replacement'
> dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
> summary(dclus2wr)
2 - level Cluster Sampling design
With (40, 126) clusters.
svydesign(id = ~dnum + snum, weights = ~pw, data = apiclus2)
Probabilities:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.003669 0.037740 0.052840 0.042390 0.052840 0.052840 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc1"     "fpc2"    
> svymean(~api00, dclus2wr)
        mean     SE
api00 670.81 30.712
> svytotal(~enroll, dclus2wr, na.rm=TRUE)
         total     SE
enroll 2639273 820261
> 
> # convert to replicate weights
> rclus1<-as.svrepdesign(dclus1)
> summary(rclus1)
Call: as.svrepdesign(dclus1)
Unstratified cluster jacknife (JK1) with 15 replicates.
Variables: 
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> svymean(~api00, rclus1)
        mean     SE
api00 644.17 26.329
> svytotal(~enroll, rclus1, na.rm=TRUE)
          mean     SE
enroll 3404940 932235
> 
> # post-stratify on school type
> pop.types<-xtabs(~stype, data=apipop)
> 
> rclus1p<-postStratify(rclus1, ~stype, pop.types)
> dclus1p<-postStratify(dclus1, ~stype, pop.types)
> summary(dclus1p)
1 - level Cluster Sampling design
With (15) clusters.
postStratify(dclus1, ~stype, pop.types)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01854 0.03257 0.03257 0.03040 0.03257 0.03257 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> summary(rclus1p)
Call: postStratify(rclus1, ~stype, pop.types)
Unstratified cluster jacknife (JK1) with 15 replicates.
Variables: 
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> 
> svymean(~api00, dclus1p)
        mean     SE
api00 642.31 23.921
> svytotal(~enroll, dclus1p, na.rm=TRUE)
         total     SE
enroll 3680893 406293
> 
> svymean(~api00, rclus1p)
        mean    SE
api00 642.31 26.45
> svytotal(~enroll, rclus1p, na.rm=TRUE)
          mean     SE
enroll 3680893 346014
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "as.fpc"
> 
> ### * as.fpc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.fpc
> ### Title: Package sample and population size data
> ### Aliases: as.fpc
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "as.svrepdesign"
> 
> ### * as.svrepdesign
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.svrepdesign
> ### Title: Convert a survey design to use replicate weights
> ### Aliases: as.svrepdesign
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> scddes<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE, fpc=rep(5,6))
> scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE)
> 
> # convert to BRR replicate weights
> scd2brr <- as.svrepdesign(scdnofpc, type="BRR")
> scd2fay <- as.svrepdesign(scdnofpc, type="Fay",fay.rho=0.3)
> # convert to JKn weights 
> scd2jkn <- as.svrepdesign(scdnofpc, type="JKn")
> 
> # convert to JKn weights with finite population correction
> scd2jknf <- as.svrepdesign(scddes, type="JKn")
> 
> ## with user-supplied hadamard matrix
> scd2brr1 <- as.svrepdesign(scdnofpc, type="BRR", hadamard.matrix=paley(11))
> 
> svyratio(~alive, ~arrests, design=scd2brr)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2brr)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> svyratio(~alive, ~arrests, design=scd2brr1)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2brr1)
Ratios=
        arrests
alive 0.1535064
SEs=
           [,1]
[1,] 0.01001468
> svyratio(~alive, ~arrests, design=scd2fay)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2fay)
Ratios=
        arrests
alive 0.1535064
SEs=
            [,1]
[1,] 0.009525187
> svyratio(~alive, ~arrests, design=scd2jkn)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2jkn)
Ratios=
        arrests
alive 0.1535064
SEs=
            [,1]
[1,] 0.009846457
> svyratio(~alive, ~arrests, design=scd2jknf)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2jknf)
Ratios=
        arrests
alive 0.1535064
SEs=
            [,1]
[1,] 0.007627033
> 
> data(api)
> ## one-stage cluster sample
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> ## convert to JK1 jackknife
> rclus1<-as.svrepdesign(dclus1)
> ## convert to bootstrap
> bclus1<-as.svrepdesign(dclus1,type="bootstrap", replicates=100)
> 
> svymean(~api00, dclus1)
        mean     SE
api00 644.17 23.542
> svytotal(~enroll, dclus1)
         total     SE
enroll 3404940 932235
> 
> svymean(~api00, rclus1)
        mean     SE
api00 644.17 26.329
> svytotal(~enroll, rclus1)
          mean     SE
enroll 3404940 932235
> 
> svymean(~api00, bclus1)
        mean     SE
api00 644.17 24.589
> svytotal(~enroll, bclus1)
          mean     SE
enroll 3404940 896130
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "as.svydesign2"
> 
> ### * as.svydesign2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.svydesign2
> ### Title: Update to the new survey design format
> ### Aliases: as.svydesign2 .svycheck
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bootweights"
> 
> ### * bootweights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bootweights
> ### Title: Compute survey bootstrap weights
> ### Aliases: bootweights bootstratum
> ### Keywords: survey
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "brrweights"
> 
> ### * brrweights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: brrweights
> ### Title: Compute replicate weights
> ### Aliases: jk1weights jknweights brrweights
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE)
> 
> ## convert to BRR replicate weights
> scd2brr <- as.svrepdesign(scdnofpc, type="BRR")
> svymean(~alive, scd2brr)
        mean     SE
alive 46.333 3.5824
> svyratio(~alive, ~arrests, scd2brr)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, scd2brr)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> 
> ## with user-supplied hadamard matrix
> scd2brr1 <- as.svrepdesign(scdnofpc, type="BRR", hadamard.matrix=paley(11))
> svymean(~alive, scd2brr1)
        mean     SE
alive 46.333 3.5824
> svyratio(~alive, ~arrests, scd2brr1)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, scd2brr1)
Ratios=
        arrests
alive 0.1535064
SEs=
           [,1]
[1,] 0.01001468
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "calibrate"
> 
> ### * calibrate
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: calibrate
> ### Title: G-calibration (GREG) estimators
> ### Aliases: calibrate.survey.design2 calibrate.svyrep.design is.calibrated
> ###   calibrate
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> 
> pop.totals<-c(`(Intercept)`=6194, stypeH=755, stypeM=1018)
> 
> ## For a single factor variable this is equivalent to
> ## postStratify
> 
> (dclus1g<-calibrate(dclus1, ~stype, pop.totals))
1 - level Cluster Sampling design
With (15) clusters.
calibrate(dclus1, ~stype, pop.totals)
> 
> svymean(~api00, dclus1g)
        mean     SE
api00 642.31 23.921
> svytotal(~enroll, dclus1g)
         total     SE
enroll 3680893 406293
> svytotal(~stype, dclus1g)
       total        SE
stypeE  4421 3.999e-14
stypeH   755 1.382e-13
stypeM  1018 1.459e-13
> 
> ## Now add sch.wide
> (dclus1g2 <- calibrate(dclus1, ~stype+sch.wide, c(pop.totals, sch.wideYes=5122)))
1 - level Cluster Sampling design
With (15) clusters.
calibrate(dclus1, ~stype + sch.wide, c(pop.totals, sch.wideYes = 5122))
> 
> svymean(~api00, dclus1g2)
      mean     SE
api00  641 23.829
> svytotal(~enroll, dclus1g2)
         total     SE
enroll 3654414 403074
> svytotal(~stype, dclus1g2)
       total        SE
stypeE  4421 1.918e-13
stypeH   755 1.518e-13
stypeM  1018 1.251e-13
> 
> ## Finally, calibrate on 1999 API and school type
> 
> (dclus1g3 <- calibrate(dclus1, ~stype+api99, c(pop.totals, api99=3914069)))
1 - level Cluster Sampling design
With (15) clusters.
calibrate(dclus1, ~stype + api99, c(pop.totals, api99 = 3914069))
> 
> svymean(~api00, dclus1g3)
        mean     SE
api00 665.31 3.4418
> svytotal(~enroll, dclus1g3)
         total     SE
enroll 3638487 385524
> svytotal(~stype, dclus1g3)
       total        SE
stypeE  4421 1.669e-13
stypeH   755 1.217e-13
stypeM  1018 1.136e-13
> 
> ## Same syntax with replicate weights
> rclus1<-as.svrepdesign(dclus1)
> 
> (rclus1g3 <- calibrate(rclus1, ~stype+api99, c(pop.totals, api99=3914069)))
Call: calibrate(rclus1, ~stype + api99, c(pop.totals, api99 = 3914069))
Unstratified cluster jacknife (JK1) with 15 replicates.
> 
> svymean(~api00, rclus1g3)
        mean     SE
api00 665.31 3.9429
> svytotal(~enroll, rclus1g3)
          mean     SE
enroll 3638487 478749
> svytotal(~stype, rclus1g3)
       mean        SE
stypeE 4421 5.146e-12
stypeH  755 5.327e-13
stypeM 1018 7.455e-13
> 
> ## Ratio estimators
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> rstrat<-as.svrepdesign(dstrat)
> 
> svytotal(~api.stu,dstrat)
          total    SE
api.stu 3086009 99477
> 
> common<-svyratio(~api.stu, ~enroll, dstrat, separate=FALSE)
> sep<-svyratio(~api.stu,~enroll, dstrat,separate=TRUE)
> stratum.totals<-list(E=1877350, H=1013824, M=920298)
> predict(sep, total=stratum.totals)
$total
         enroll
api.stu 3190022

$se
          enroll
api.stu 29756.44

> predict(common, total=do.call("sum",stratum.totals))
$total
         enroll
api.stu 3190038

$se
          enroll
api.stu 29565.98

> 
> pop<-colSums(model.matrix(~stype*enroll-1,model.frame(~stype*enroll,apipop)))
> pop
       stypeE        stypeH        stypeM        enroll stypeH:enroll 
         4397           751          1009       3811472       1013824 
stypeM:enroll 
       920298 
> ## common ratio
> dstratg1<-calibrate(dstrat,~enroll-1, pop[4], lambda=1)
> svytotal(~api.stu, dstratg1)
          total    SE
api.stu 3190038 29566
> rstratg1<-calibrate(rstrat,~enroll-1, pop[4], lambda=1)
> svytotal(~api.stu, rstratg1)
           mean    SE
api.stu 3086121 29625
> 
> ## similar (but not identical) to separate ratio.
> dstratg2<-calibrate(dstrat,~stype*enroll-1, pop,lambda=c(0,0,0,1,0,0))
> svytotal(~api.stu,dstratg2)
          total    SE
api.stu 3189234 30768
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "compressWeights"
> 
> ### * compressWeights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: compressWeights
> ### Title: Compress replicate weight matrix
> ### Aliases: compressWeights compressWeights.default
> ###   compressWeights.repweights_compressed [.repweights_compressed
> ###   dim.repweights_compressed dimnames.repweights_compressed
> ###   as.matrix.repweights_compressed as.matrix.repweights
> ###   as.vector.repweights_compressed compressWeights.svyrep.design
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> rclus1c<-as.svrepdesign(dclus1,compress=TRUE)
> rclus1<-as.svrepdesign(dclus1,compress=FALSE)
> 
> 
> 
> cleanEx(); ..nameEx <- "crowd"
> 
> ### * crowd
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: crowd
> ### Title: Household crowding
> ### Aliases: crowd
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(crowd)
> 
> ## Example 1-1
> i1.1<-as.svrepdesign(svydesign(id=~cluster, weight=~weight,data=crowd))
> i1.1<-update(i1.1, room.ratio=rooms/person,
+ overcrowded=factor(person>rooms))
> svymean(~rooms+person+room.ratio,i1.1)
             mean     SE
rooms      6.0000 0.6831
person     4.0000 1.1832
room.ratio 2.3274 0.6006
> svytotal(~rooms+person+room.ratio,i1.1)
             mean     SE
rooms      36.000 4.0988
person     24.000 7.0993
room.ratio 13.964 3.6037
> svymean(~rooms+person+room.ratio,subset(i1.1,overcrowded==TRUE))
              mean     SE
rooms      5.50000 0.6455
person     7.50000 0.6455
room.ratio 0.73214 0.0231
> svytotal(~rooms+person+room.ratio,subset(i1.1,overcrowded==TRUE))
              mean     SE
rooms      11.0000 7.0000
person     15.0000 9.5184
room.ratio  1.4643 0.9265
> 
> ## Example 1-2
> i1.2<-as.svrepdesign(svydesign(id=~cluster,weight=~weight,strata=~stratum, data=crowd))
> svymean(~rooms+person,i1.2)
       mean     SE
rooms     6 0.2357
person    4 0.4082
> svytotal(~rooms+person,i1.2)
       mean     SE
rooms    36 1.4142
person   24 2.4495
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "fpc"
> 
> ### * fpc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fpc
> ### Title: Small survey example
> ### Aliases: fpc
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(fpc)
> fpc
  stratid psuid weight nh Nh   x
1       1     1      3  5 15 2.8
2       1     2      3  5 15 4.1
3       1     3      3  5 15 6.8
4       1     4      3  5 15 6.8
5       1     5      3  5 15 9.2
6       2     1      4  3 12 3.7
7       2     2      4  3 12 6.6
8       2     3      4  3 12 4.2
> 
> withoutfpc<-svydesign(weights=~weight, ids=~psuid, strata=~stratid, variables=~x, data=fpc, nest=TRUE)
> 
> withoutfpc
Stratified Independent Sampling design
svydesign(weights = ~weight, ids = ~psuid, strata = ~stratid, 
    variables = ~x, data = fpc, nest = TRUE)
> svymean(~x, withoutfpc)
    mean     SE
x 5.4481 0.7413
> 
> withfpc<-svydesign(weights=~weight, ids=~psuid, strata=~stratid,
+ fpc=~Nh, variables=~x, data=fpc, nest=TRUE)
> 
> withfpc
Stratified Independent Sampling design
svydesign(weights = ~weight, ids = ~psuid, strata = ~stratid, 
    fpc = ~Nh, variables = ~x, data = fpc, nest = TRUE)
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> ## Other equivalent forms 
> withfpc<-svydesign(prob=~I(1/weight), ids=~psuid, strata=~stratid,
+ fpc=~Nh, variables=~x, data=fpc, nest=TRUE)
> 
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> withfpc<-svydesign(weights=~weight, ids=~psuid, strata=~stratid,
+ fpc=~I(nh/Nh), variables=~x, data=fpc, nest=TRUE)
> 
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> withfpc<-svydesign(weights=~weight, ids=~interaction(stratid,psuid),
+ strata=~stratid, fpc=~I(nh/Nh), variables=~x, data=fpc)
> 
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> withfpc<-svydesign(ids=~psuid, strata=~stratid, fpc=~Nh,
+  variables=~x,data=fpc,nest=TRUE)
> 
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> withfpc<-svydesign(ids=~psuid, strata=~stratid,
+ fpc=~I(nh/Nh), variables=~x, data=fpc, nest=TRUE)
> 
> svymean(~x, withfpc)
    mean    SE
x 5.4481 0.616
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "ftable.svystat"
> 
> ### * ftable.svystat
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ftable.svystat
> ### Title: Lay out tables of survey statistics
> ### Aliases: ftable.svystat ftable.svrepstat ftable.svyby
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> 
> a<-svymean(~interaction(stype,comp.imp), design=dclus1)
> b<-ftable(a, rownames=list(stype=c("E","H","M"),comp.imp=c("No","Yes")))
> b
              stype          E          H          M
comp.imp                                            
No       mean       0.17486339 0.03825137 0.06010929
         SE         0.02599552 0.01607602 0.02457177
Yes      mean       0.61202186 0.03825137 0.07650273
         SE         0.04167572 0.01605469 0.02167084
> round(100*b,1)
              stype    E    H    M
comp.imp                          
No       mean       17.5  3.8  6.0
         SE          2.6  1.6  2.5
Yes      mean       61.2  3.8  7.7
         SE          4.2  1.6  2.2
> 
> a<-svymean(~interaction(stype,comp.imp), design=dclus1, deff=TRUE)
> b<-ftable(a, rownames=list(stype=c("E","H","M"),comp.imp=c("No","Yes")))
> b
              stype          E          H          M
comp.imp                                            
No       mean       0.17486339 0.03825137 0.06010929
         SE         0.02599552 0.01607602 0.02457177
         Deff       0.88317737 1.32472009 2.01525343
Yes      mean       0.61202186 0.03825137 0.07650273
         SE         0.04167572 0.01605469 0.02167084
         Deff       1.37932843 1.32120655 1.25347214
> round(100*b,1)
              stype     E     H     M
comp.imp                             
No       mean        17.5   3.8   6.0
         SE           2.6   1.6   2.5
         Deff        88.3 132.5 201.5
Yes      mean        61.2   3.8   7.7
         SE           4.2   1.6   2.2
         Deff       137.9 132.1 125.3
> 
> rclus1<-as.svrepdesign(dclus1)
> a<-svytotal(~interaction(stype,comp.imp), design=rclus1)
> b<-ftable(a, rownames=list(stype=c("E","H","M"),comp.imp=c("No","Yes")))
> b
              stype          E          H          M
comp.imp                                            
No       mean       1083.10388  236.92897  372.31696
         SE          373.15133  108.21477  114.69172
Yes      mean       3790.86359  236.92897  473.85795
         SE          994.30863   83.05377  166.10754
> round(b)
              stype    E    H    M
comp.imp                          
No       mean       1083  237  372
         SE          373  108  115
Yes      mean       3791  237  474
         SE          994   83  166
> 
> a<-svyby(~api99 + api00, ~stype + sch.wide, rclus1, svymean, keep.var=TRUE)
> ftable(a)
              sch.wide        No                 Yes          
                           api99     api00     api99     api00
stype                                                         
E     svymean          601.66667 596.33333 577.63636 607.45455
      SE                70.04669  64.50553  57.38815  53.97142
H     svymean          608.34848 653.64394 611.37500 606.37500
      SE                23.67277  22.37296  48.19716  48.27853
M     svymean          662.00000 659.33333 607.29412 643.23529
      SE                40.92204  37.80385  49.49574  49.34813
> print(ftable(a),digits=2)
              sch.wide    No         Yes      
                       api99 api00 api99 api00
stype                                         
E     svymean            602   596   578   607
      SE                  70    65    57    54
H     svymean            608   654   611   606
      SE                  24    22    48    48
M     svymean            662   659   607   643
      SE                  41    38    49    49
> 
> b<-svyby(~api99 + api00, ~stype + sch.wide, rclus1, svymean, keep.var=TRUE, deff=TRUE)
> ftable(b)
              sch.wide         No                   Yes           
                            api99      api00      api99      api00
stype                                                             
E     svymean          601.666667 596.333333 577.636364 607.454545
      SE                70.046693  64.505534  57.388148  53.971422
      DEff               6.887863   7.169080   3.317109   3.329886
H     svymean          608.348485 653.643939 611.375000 606.375000
      SE                23.672765  22.372956  48.197160  48.278530
      DEff               5.692758   6.033725   2.033992   2.062942
M     svymean          662.000000 659.333333 607.294118 643.235294
      SE                40.922039  37.803855  49.495744  49.348127
      DEff               2.121054   2.121054   3.332404   3.562371
> print(ftable(b),digits=2)
              sch.wide    No         Yes      
                       api99 api00 api99 api00
stype                                         
E     svymean          601.7 596.3 577.6 607.5
      SE                70.0  64.5  57.4  54.0
      DEff               6.9   7.2   3.3   3.3
H     svymean          608.3 653.6 611.4 606.4
      SE                23.7  22.4  48.2  48.3
      DEff               5.7   6.0   2.0   2.1
M     svymean          662.0 659.3 607.3 643.2
      SE                40.9  37.8  49.5  49.3
      DEff               2.1   2.1   3.3   3.6
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "hadamard"
> 
> ### * hadamard
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hadamard
> ### Title: Hadamard matrices
> ### Aliases: hadamard
> ### Keywords: survey
> 
> ### ** Examples
> 
> 
> par(mfrow=c(2,2))
> ## Sylvester-type
> image(hadamard(63),main="Sylvester: 64")
> ## Paley-type
> image(hadamard(57),main="Paley: 60")
> ## hybrid
> image(hadamard(23),main="Stored: 24")
> image(hadamard(90),main="Constructed: 96=4x24")
> 
> par(mfrow=c(1,1))
> plot(2:150,sapply(2:150,function(i) ncol(hadamard(i))),type="S",
+      ylab="Matrix size",xlab="n",xlim=c(1,150),ylim=c(1,150))
> abline(0,1,lty=3)
> lines(2:150, 2:150-(2:150 %% 4)+4,col="purple",type="S",lty=2)
> legend(c(x=10,y=140),legend=c("Actual size","Minimum possible size"),
+      col=c("black","purple"),bty="n",lty=c(1,2))
> 
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "hospital"
> 
> ### * hospital
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hospital
> ### Title: Sample of obstetric hospitals
> ### Aliases: hospital
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(hospital)
> hospdes<-svydesign(strata=~oblevel, id=~hospno, weights=~weighta,
+ fpc=~tothosp, data=hospital)
> hosprep<-as.svrepdesign(hospdes)
> 
> svytotal(~births, design=hospdes)
        total    SE
births 183983 34014
> svytotal(~births, design=hosprep)
         mean    SE
births 183983 34014
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "mu284"
> 
> ### * mu284
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mu284
> ### Title: Two-stage sample from MU284
> ### Aliases: mu284
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(mu284)
> (dmu284<-svydesign(id=~id1+id2,fpc=~n1+n2, data=mu284))
2 - level Cluster Sampling design
With (5, 15) clusters.
svydesign(id = ~id1 + id2, fpc = ~n1 + n2, data = mu284)
> (ytotal<-svytotal(~y1, dmu284))
   total     SE
y1 15080 2274.3
> vcov(ytotal)
        [,1]
[1,] 5172234
> 
> 
> 
> cleanEx(); ..nameEx <- "nonresponse"
> 
> ### * nonresponse
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nonresponse
> ### Title: Experimental: Construct non-response weights
> ### Aliases: nonresponse sparseCells neighbours joinCells
> ###   weights.nonresponse print.nonresponse print.nonresponseSubset
> ###   [.nonresponse
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> ## pretend the sampling was stratified on three variables
> poptable<-xtabs(~sch.wide+comp.imp+stype,data=apipop)
> sample.count<-xtabs(~sch.wide+comp.imp+stype,data=apiclus1)
> sample.weight<-xtabs(pw~sch.wide+comp.imp+stype, data=apiclus1)
> 
> ## create a nonresponse object
> nr<-nonresponse(sample.weight,sample.count, poptable)
> 
> ## sparse cells
> sparseCells(nr)
sparseCells(nr)
Cells: 3 5 7 11
Indices:
   sch.wide comp.imp stype
3  "No"     "Yes"    "E"  
5  "No"     "No"     "H"  
7  "No"     "Yes"    "H"  
11 "No"     "Yes"    "M"  
Summary:
   NRwt  wt n
3   Inf Inf 0
5   3.2 108 3
7   Inf Inf 0
11  Inf Inf 0
> 
> ## Look at neighbours
> neighbours(3,nr)
"[.nonresponse"(object, nbour.index)
Cells: 4 7 1
Indices:
  sch.wide comp.imp stype
4 "Yes"    "Yes"    "E"  
7 "No"     "Yes"    "H"  
1 "No"     "No"     "E"  
Summary:
  NRwt   wt   n
4 0.92 31.1 112
7  Inf  Inf   0
1 1.04 35.2  12
> neighbours(11,nr)
"[.nonresponse"(object, nbour.index)
Cells: 12 9 7
Indices:
   sch.wide comp.imp stype
12 "Yes"    "Yes"    "M"  
9  "No"     "No"     "M"  
7  "No"     "Yes"    "H"  
Summary:
    NRwt   wt  n
12 1.290 43.6 14
9  0.916 31.0  8
7    Inf  Inf  0
> 
> ## Collapse some contiguous cells
> nr1<-joinCells(nr,3,5,7)
> 
> ## sparse cells now
> sparseCells(nr1)
sparseCells(nr1)
Cells: 3 11
Indices:
   sch.wide comp.imp stype
3  "No"     "Yes"    "E"  
11 "No"     "Yes"    "M"  
Summary:
   NRwt  wt n
3  3.78 128 3
11  Inf Inf 0
> nr2<-joinCells(nr1,3,11,8)
> 
> nr2
Call: nonresponse(sample.weight, sample.count, poptable)
12 original cells,  8 distinct cells remaining
Joins:
3 5 7 
3 5 7 8 11 
     counts         NRweights         totalwts    
 Min.   :  3.00   Min.   :0.6840   Min.   :23.15  
 1st Qu.:  7.00   1st Qu.:0.8956   1st Qu.:30.31  
 Median : 11.00   Median :0.9793   Median :33.15  
 Mean   : 22.88   Mean   :1.1461   Mean   :38.79  
 3rd Qu.: 15.50   3rd Qu.:1.3142   3rd Qu.:44.48  
 Max.   :112.00   Max.   :2.0977   Max.   :71.00  
> 
> ## one relatively sparse cell
> sparseCells(nr2)
sparseCells(nr2)
Cells: 3
Indices:
  sch.wide comp.imp stype
3 "No"     "Yes"    "E"  
Summary:
  NRwt wt  n
3  2.1 71 10
> ## but nothing suitable to join it to
> neighbours(3,nr2)
"[.nonresponse"(object, nbour.index)
Cells: 4 1 6 9 12
Indices:
   sch.wide comp.imp stype
4  "Yes"    "Yes"    "E"  
1  "No"     "No"     "E"  
6  "Yes"    "No"     "H"  
9  "No"     "No"     "M"  
12 "Yes"    "Yes"    "M"  
Summary:
    NRwt   wt   n
4  0.920 31.1 112
1  1.040 35.2  12
6  0.835 28.3   4
9  0.916 31.0   8
12 1.290 43.6  14
> 
> ## extract the weights
> weights(nr2)
, , stype = E

        comp.imp
sch.wide        No       Yes
     No  1.0389893 2.0976751
     Yes 0.6839602 0.9195794

, , stype = H

        comp.imp
sch.wide        No      Yes
     No  2.0976751 2.097675
     Yes 0.8346383 2.097675

, , stype = M

        comp.imp
sch.wide        No      Yes
     No  0.9158863 2.097675
     Yes 1.3886018 1.289416

> 
> 
> 
> cleanEx(); ..nameEx <- "paley"
> 
> ### * paley
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: paley
> ### Title: Paley-type Hadamard matrices
> ### Aliases: paley
> ### Keywords: survey algebra
> 
> ### ** Examples
> 
> 
> p<-paley(11)
> 
> p1<-2*p-1
> ## HH^T is diagonal for any Hadamard matrix
> p1%*%t(p1)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]   12    0    0    0    0    0    0    0    0     0     0     0
 [2,]    0   12    0    0    0    0    0    0    0     0     0     0
 [3,]    0    0   12    0    0    0    0    0    0     0     0     0
 [4,]    0    0    0   12    0    0    0    0    0     0     0     0
 [5,]    0    0    0    0   12    0    0    0    0     0     0     0
 [6,]    0    0    0    0    0   12    0    0    0     0     0     0
 [7,]    0    0    0    0    0    0   12    0    0     0     0     0
 [8,]    0    0    0    0    0    0    0   12    0     0     0     0
 [9,]    0    0    0    0    0    0    0    0   12     0     0     0
[10,]    0    0    0    0    0    0    0    0    0    12     0     0
[11,]    0    0    0    0    0    0    0    0    0     0    12     0
[12,]    0    0    0    0    0    0    0    0    0     0     0    12
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "postStratify"
> 
> ### * postStratify
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: postStratify
> ### Title: Post-stratify a survey
> ### Aliases: postStratify postStratify.svyrep.design
> ###   postStratify.survey.design
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> rclus1<-as.svrepdesign(dclus1)
> 
> svymean(~api00, rclus1)
        mean     SE
api00 644.17 26.329
> svytotal(~enroll, rclus1)
          mean     SE
enroll 3404940 932235
> 
> # post-stratify on school type
> pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
> #or: pop.types <- xtabs(~stype, data=apipop)
> #or: pop.types <- table(stype=apipop$stype)
> 
> rclus1p<-postStratify(rclus1, ~stype, pop.types)
> summary(rclus1p)
Call: postStratify(rclus1, ~stype, pop.types)
Unstratified cluster jacknife (JK1) with 15 replicates.
Variables: 
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> svymean(~api00, rclus1p)
        mean    SE
api00 642.31 26.45
> svytotal(~enroll, rclus1p)
          mean     SE
enroll 3680893 346014
> 
> ## and for svydesign objects
> dclus1p<-postStratify(dclus1, ~stype, pop.types)
> summary(dclus1p)
1 - level Cluster Sampling design
With (15) clusters.
postStratify(dclus1, ~stype, pop.types)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01854 0.03257 0.03257 0.03040 0.03257 0.03257 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> svymean(~api00, dclus1p)
        mean     SE
api00 642.31 23.921
> svytotal(~enroll, dclus1p)
         total     SE
enroll 3680893 406293
> 
> 
> 
> cleanEx(); ..nameEx <- "rake"
> 
> ### * rake
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rake
> ### Title: Raking of replicate weight design
> ### Aliases: rake
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> rclus1 <- as.svrepdesign(dclus1)
> 
> svymean(~api00, rclus1)
        mean     SE
api00 644.17 26.329
> svytotal(~enroll, rclus1)
          mean     SE
enroll 3404940 932235
> 
> ## 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.871
> svytotal(~enroll, rclus1r)
          mean     SE
enroll 3647300 463253
> 
> ## marginal totals correspond to population
> xtabs(~stype, apipop)
stype
   E    H    M 
4421  755 1018 
> svytable(~stype, rclus1r, round=TRUE)
stype
   E    H    M 
4421  755 1018 
> xtabs(~sch.wide, apipop)
sch.wide
  No  Yes 
1072 5122 
> svytable(~sch.wide, rclus1r, round=TRUE)
sch.wide
  No  Yes 
1072 5122 
> 
> ## joint totals don't correspond 
> xtabs(~stype+sch.wide, apipop)
     sch.wide
stype   No  Yes
    E  472 3949
    H  334  421
    M  266  752
> svytable(~stype+sch.wide, rclus1r, round=TRUE)
     sch.wide
stype   No  Yes
    E  478 3943
    H  201  554
    M  393  625
> 
> ## 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.739
> svytotal(~enroll, dclus1r)
         total     SE
enroll 3647300 399094
> 
> ## compare to joint post-stratification
> ## (only possible if joint population table is known)
> ##
> pop.table <- xtabs(~stype+sch.wide,apipop)
> rclus1ps <- postStratify(rclus1, ~stype+sch.wide, pop.table)
> svytable(~stype+sch.wide, rclus1ps, round=TRUE)
     sch.wide
stype   No  Yes
    E  472 3949
    H  334  421
    M  266  752
> 
> svymean(~api00, rclus1ps)
        mean     SE
api00 643.15 26.661
> svytotal(~enroll, rclus1ps)
          mean     SE
enroll 3610219 346688
> 
> ## Example of raking with partial joint distributions
> 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.776
> 
> 
> 
> cleanEx(); ..nameEx <- "regTermTest"
> 
> ### * regTermTest
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: regTermTest
> ### Title: Wald test for a term in a regression model
> ### Aliases: regTermTest print.regTermTest
> ### Keywords: regression
> 
> ### ** Examples
> 
>  data(esoph)
>  model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * 
+      alcgp, data = esoph, family = binomial())
>  anova(model1)
Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(ncases, ncontrols)

Terms added sequentially (first to last)


            Df Deviance Resid. Df Resid. Dev
NULL                           87    227.241
agegp        5   88.128        82    139.112
tobgp        3   19.085        79    120.028
alcgp        3   66.054        76     53.973
tobgp:alcgp  9    6.489        67     47.484
> 
>  regTermTest(model1,"tobgp")
Wald test for tobgp
 in glm(formula = cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, 
    family = binomial(), data = esoph)
Chisq =  11.88584  on  3  df: p= 0.0077846 
>  regTermTest(model1,"tobgp:alcgp")
Wald test for tobgp:alcgp
 in glm(formula = cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, 
    family = binomial(), data = esoph)
Chisq =  6.205544  on  9  df: p= 0.71918 
>  regTermTest(model1, ~alcgp+tobgp:alcgp)
Wald test for alcgp alcgp:tobgp
 in glm(formula = cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, 
    family = binomial(), data = esoph)
Chisq =  60.2626  on  12  df: p= 2.0217e-08 
> 
> 
> 
> cleanEx(); ..nameEx <- "scd"
> 
> ### * scd
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: scd
> ### Title: Survival in cardiac arrest
> ### Aliases: scd
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(scd)
> 
> ## survey design objects
> scddes<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE, fpc=rep(5,6))
> scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE)
> 
> # convert to BRR replicate weights
> scd2brr <- as.svrepdesign(scdnofpc, type="BRR")
> 
> # use BRR replicate weights from Levy and Lemeshow
> repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+ c(0,1,0,1,1,0))
> scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
> 
> # ratio estimates
> svyratio(~alive, ~arrests, design=scddes)
Ratio estimator: svyratio.survey.design2(~alive, ~arrests, design = scddes)
Ratios=
        arrests
alive 0.1535064
SEs=
          arrests
alive 0.007596705
> svyratio(~alive, ~arrests, design=scdnofpc)
Ratio estimator: svyratio.survey.design2(~alive, ~arrests, design = scdnofpc)
Ratios=
        arrests
alive 0.1535064
SEs=
          arrests
alive 0.009807304
> svyratio(~alive, ~arrests, design=scd2brr)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2brr)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> svyratio(~alive, ~arrests, design=scdrep)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scdrep)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> 
> # or a logistic regression
> summary(svyglm(cbind(alive,arrests-alive)~1, family=quasibinomial, design=scdnofpc))

Call:
svyglm.survey.design(formula = cbind(alive, arrests - alive) ~ 
    1, design = scdnofpc, family = quasibinomial)

Survey design:
svydesign(data = scd, prob = ~1, id = ~ambulance, strata = ~ESA, 
    nest = TRUE)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.70736    0.07547  -22.62 3.14e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for quasibinomial family taken to be 4.704437)

Number of Fisher Scoring iterations: 4

> summary(svyglm(cbind(alive,arrests-alive)~1, family=quasibinomial, design=scdrep))

Call:
svyglm.svyrep.design(formula = cbind(alive, arrests - alive) ~ 
    1, design = scdrep, family = quasibinomial)

Survey design:
svrepdesign(data = scd, type = "BRR", repweights = repweights)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.7074     0.1241  -13.75 3.65e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for quasibinomial family taken to be 5.273208)

Number of Fisher Scoring iterations: 4

> 
> 
> 
> 
> cleanEx(); ..nameEx <- "subset.survey.design"
> 
> ### * subset.survey.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: subset.survey.design
> ### Title: Subset of survey
> ### Aliases: subset.survey.design subset.svyrep.design [.survey.design
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(fpc)
> dfpc<-svydesign(id=~psuid,strat=~stratid,weight=~weight,data=fpc,nest=TRUE)
> dsub<-subset(dfpc,x>4)
> summary(dsub)
Stratified Independent Sampling design
subset(dfpc, x > 4)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2500  0.2708  0.3333  0.3056  0.3333  0.3333 
Stratum Sizes: 
           1 2
obs        4 2
design.PSU 5 3
actual.PSU 4 2
Data variables:
[1] "stratid" "psuid"   "weight"  "nh"      "Nh"      "x"      
> svymean(~x,design=dsub)
   mean     SE
x 6.195 0.7555
> 
> ## These should give the same domain estimates and standard errors
> svyby(~x,~I(x>4),design=dfpc, svymean)
Warning in onestrat(x[index, , drop = FALSE], clusters[index], nPSU[index][1],  : 
	 Stratum (1) has only one PSU at stage 1
Warning in onestrat(x[index, , drop = FALSE], clusters[index], nPSU[index][1],  : 
	 Stratum (2) has only one PSU at stage 1
      I(x > 4) statistic.x        SE
FALSE    FALSE    3.314286 0.3117042
TRUE      TRUE    6.195000 0.7555129
> summary(svyglm(x~I(x>4)+0,design=dfpc))

Call:
svyglm.survey.design(formula = x ~ I(x > 4) + 0, design = dfpc)

Survey design:
svydesign(id = ~psuid, strat = ~stratid, weight = ~weight, data = fpc, 
    nest = TRUE)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
I(x > 4)FALSE   3.3143     0.3117   10.63 4.08e-05 ***
I(x > 4)TRUE    6.1950     0.7555    8.20 0.000177 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 2.237706)

Number of Fisher Scoring iterations: 2

> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> rclus1<-as.svrepdesign(dclus1)
> svymean(~enroll, subset(dclus1, sch.wide=="Yes" & comp.imp=="Yes"))
         mean     SE
enroll 534.56 36.248
> svymean(~enroll, subset(rclus1, sch.wide=="Yes" & comp.imp=="Yes"))
         mean     SE
enroll 534.56 40.398
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "surveysummary"
> 
> ### * surveysummary
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: surveysummary
> ### Title: Summary statistics for sample surveys
> ### Aliases: svymean svymean.survey.design svymean.survey.design2
> ###   svymean.svyrep.design svrepmean svytotal svytotal.survey.design
> ###   svytotal.survey.design2 svytotal.svyrep.design print.svystat
> ###   print.svrepstat svreptotal svyvar svyvar.survey.design
> ###   svyvar.svyrep.design coef.svystat vcov.svystat coef.svrepstat
> ###   vcov.svrepstat cv.svyratio cv.svrepratio cv.svrepstat cv.svystat
> ###   cv.default cv deff deff.default svrepvar
> ### Keywords: univar survey
> 
> ### ** Examples
> 
> 
>   data(api)
>   ## population
>   mean(apipop$api00)
[1] 664.7126
>   quantile(apipop$api00,c(.25,.5,.75))
25% 50% 75% 
565 667 761 
>   var(apipop$api00)
[1] 16446.56
>   sum(apipop$enroll)
[1] NA
>   sum(apipop$api.stu)/sum(apipop$enroll)
[1] NA
> 
>   ## one-stage cluster sample
>   dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
>   summary(dclus1)
1 - level Cluster Sampling design
With (15) clusters.
svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
>   svymean(~api00, dclus1, deff=TRUE)
         mean      SE   DEff
api00 644.169  23.542 9.3972
>   svymean(~factor(stype),dclus1)
                   mean     SE
factor(stype)E 0.786885 0.0463
factor(stype)H 0.076503 0.0268
factor(stype)M 0.136612 0.0296
>   svymean(~interaction(stype, comp.imp), dclus1)
                                      mean     SE
interaction(stype, comp.imp)E.No  0.174863 0.0260
interaction(stype, comp.imp)H.No  0.038251 0.0161
interaction(stype, comp.imp)M.No  0.060109 0.0246
interaction(stype, comp.imp)E.Yes 0.612022 0.0417
interaction(stype, comp.imp)H.Yes 0.038251 0.0161
interaction(stype, comp.imp)M.Yes 0.076503 0.0217
>   svyquantile(~api00, dclus1, c(.25,.5,.75))
        0.25 0.5  0.75
api00 551.75 652 717.5
>   svyvar(~api00, dclus1)
      variance     SE
api00    11122 1378.8
>   svytotal(~enroll, dclus1, deff=TRUE)
         total      SE   DEff
enroll 3404940  932235 30.552
>   svyratio(~api.stu, ~enroll, dclus1)
Ratio estimator: svyratio.survey.design2(~api.stu, ~enroll, dclus1)
Ratios=
           enroll
api.stu 0.8497087
SEs=
             enroll
api.stu 0.008386297
> 
>   #stratified sample
>   dstrat<-svydesign(id=~1, strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
>   summary(dstrat)
Stratified Independent Sampling design
svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02262 0.02262 0.03587 0.04014 0.05339 0.06623 
Stratum Sizes: 
             E  H  M
obs        100 50 50
design.PSU 100 50 50
actual.PSU 100 50 50
Population stratum sizes (PSUs): 
   E    M    H 
4421 1018  755 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc"     
>   svymean(~api00, dstrat)
        mean     SE
api00 662.29 9.4089
>   svyquantile(~api00, dstrat, c(.25,.5,.75))
          0.25      0.5     0.75
api00 562.2056 667.2358 755.1226
>   svyvar(~api00, dstrat)
      variance     SE
api00    15115 1249.5
>   svytotal(~enroll, dstrat)
         total     SE
enroll 3687178 114642
>   svyratio(~api.stu, ~enroll, dstrat)
Ratio estimator: svyratio.survey.design2(~api.stu, ~enroll, dstrat)
Ratios=
           enroll
api.stu 0.8369569
SEs=
             enroll
api.stu 0.007757103
>   
>   # replicate weights - jackknife (this is slow)
>   jkstrat<-as.svrepdesign(dstrat)
>   summary(jkstrat)
Call: as.svrepdesign(dstrat)
Stratified cluster jackknife (JKn) with 200 replicates.
Variables: 
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc"     
>   svymean(~api00, jkstrat)
        mean     SE
api00 662.29 9.4089
>   svymean(~factor(stype),jkstrat)
                  mean        SE
factor(stype)E 0.71376 2.331e-14
factor(stype)H 0.12189 4.761e-15
factor(stype)M 0.16435 1.972e-15
>   svyvar(~api00,jkstrat)
      variance     SE
api00    15115 1257.5
>   svyquantile(~api00, jkstrat, c(.25,.5,.75))
Statistic:
       api00
q0.25 562.21
q0.5  667.24
q0.75 755.12
SE:
         api00
q0.25 12.29485
q0.5  11.63687
q0.75 15.70636
>   svytotal(~enroll, jkstrat)
          mean     SE
enroll 3687178 114642
>   svyratio(~api.stu, ~enroll, jkstrat)
Ratio estimator: svyratio.svyrep.design(~api.stu, ~enroll, jkstrat)
Ratios=
           enroll
api.stu 0.8369569
SEs=
            [,1]
[1,] 0.007772509
> 
>   # coefficients of variation
>   cv(svytotal(~enroll,dstrat))
  enroll 
0.031092 
>   cv(svyratio(~api.stu, ~enroll, jkstrat))
            enroll
api.stu 0.00928663
> 
>   # extracting statistic and variance
>   coef(svytotal(~enroll,dstrat))
 enroll 
3687178 
>   vcov(svymean(~api00+api99,jkstrat))
         [,1]     [,2]
[1,] 88.52817 91.80068
[2,] 91.80068 99.28025
> 
>   # Design effect
>   svymean(~api00, dstrat, deff=TRUE)
          mean       SE   DEff
api00 662.2874   9.4089 1.2105
>   svymean(~api00, dstrat, deff="replace")
          mean       SE   DEff
api00 662.2874   9.4089 1.1714
>   svymean(~api00, jkstrat, deff=TRUE)
          mean       SE   DEff
api00 662.2874   9.4089 1.2105
>   svymean(~api00, jkstrat, deff="replace")
          mean       SE   DEff
api00 662.2874   9.4089 1.1714
>  (a<-svytotal(~enroll, dclus1, deff=TRUE))
         total      SE   DEff
enroll 3404940  932235 30.552
>   deff(a)
[1] 30.55255
> 
>   # BRR method
>   data(scd)
>   repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+               c(0,1,0,1,1,0))
>   scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
>   svymean(~arrests+alive, design=scdrep)
           mean      SE
arrests 301.833 25.3930
alive    46.333  3.5824
> 
>  
> 
> 
> cleanEx(); ..nameEx <- "svrVar"
> 
> ### * svrVar
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svrVar
> ### Title: Compute variance from replicates
> ### Aliases: svrVar
> ### Keywords: survey
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svrepdesign"
> 
> ### * svrepdesign
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svrepdesign
> ### Title: Specify survey design with replicate weights
> ### Aliases: svrepdesign [.svyrep.design image.svyrep.design
> ###   print.svyrep.design model.frame.svyrep.design summary.svyrep.design
> ###   print.summary.svyrep.design
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> # use BRR replicate weights from Levy and Lemeshow
> repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+ c(0,1,0,1,1,0))
> scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
> svyratio(~alive, ~arrests, scdrep)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, scdrep)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> 
> 
> 
> cleanEx(); ..nameEx <- "svy.varcoef"
> 
> ### * svy.varcoef
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svy.varcoef
> ### Title: Sandwich variance estimator for glms
> ### Aliases: svy.varcoef
> ### Keywords: regression survey
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyCprod"
> 
> ### * svyCprod
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyCprod
> ### Title: Computations for survey variances
> ### Aliases: svyCprod onestage onestrat
> ### Keywords: utilities survey
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyby"
> 
> ### * svyby
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyby
> ### Title: Survey statistics on subsets
> ### Aliases: svyby
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> 
> svyby(~api99, ~stype, dclus1, svymean)
  stype statistic.api99       SE
E     E        607.7917 22.81660
H     H        595.7143 41.76400
M     M        608.6000 32.56064
> svyby(~api99, ~stype, dclus1, svyquantile, quantiles=0.5,ci=TRUE)
  stype statistic.quantiles      statistic.CIs       SE
E     E                 615 525.6174, 674.1479 37.89113
H     H                 593 428.4810, 701.0065 69.52309
M     M                 611 527.5797, 675.2395 37.66903
> ## without ci=TRUE svyquantile does not compute standard errors
> svyby(~api99, ~stype, dclus1, svyquantile, quantiles=0.5, keep.var=FALSE)
  stype statistic
E     E       615
H     H       593
M     M       611
> svyby(~api99, list(school.type=apiclus1$stype), dclus1, svymean)
  school.type statistic.api99       SE
E           E        607.7917 22.81660
H           H        595.7143 41.76400
M           M        608.6000 32.56064
> svyby(~api99+api00, ~stype, dclus1, svymean, deff=TRUE)
  stype statistic.api99 statistic.api00      SE1      SE2    DEff1    DEff2
E     E        607.7917        648.8681 22.81660 22.36241 5.936963 6.629713
H     H        595.7143        618.5714 41.76400 38.02025 2.382009 2.399664
M     M        608.6000        631.4400 32.56064 31.60947 2.319781 2.254062
> svyby(~api99+api00, ~stype+sch.wide, dclus1, svymean, keep.var=FALSE)
      stype sch.wide statistic.api99 statistic.api00
E.No      E       No        601.6667        596.3333
E.Yes     E      Yes        608.3485        653.6439
H.No      H       No        662.0000        659.3333
H.Yes     H      Yes        577.6364        607.4545
M.No      M       No        611.3750        606.3750
M.Yes     M      Yes        607.2941        643.2353
> 
> rclus1<-as.svrepdesign(dclus1)
> 
> svyby(~api99, ~stype, rclus1, svymean)
  stype statistic.api99       SE
E     E        607.7917 25.83542
H     H        595.7143 50.75106
M     M        608.6000 34.82521
> svyby(~api99, ~stype, rclus1, svyquantile, quantiles=0.5)
  stype statistic       SE
E     E       615 43.11213
H     H       593 74.97454
M     M       611 42.42604
> svyby(~api99, list(school.type=apiclus1$stype), rclus1, svymean)
  school.type statistic.api99       SE
E           E        607.7917 25.83542
H           H        595.7143 50.75106
M           M        608.6000 34.82521
> svyby(~enroll,~stype, rclus1,svytotal, deff=TRUE)
  stype statistic.enroll       SE       DEff
E     E        2109717.1 631349.4 2991144956
H     H         535594.9 226716.6    1123422
M     M         759628.1 213635.5    9707132
> svyby(~api99+api00, ~stype+sch.wide, rclus1, svymean, keep.var=FALSE)
      stype sch.wide statistic.api99 statistic.api00
E.No      E       No        601.6667        596.3333
E.Yes     E      Yes        608.3485        653.6439
H.No      H       No        662.0000        659.3333
H.Yes     H      Yes        577.6364        607.4545
M.No      M       No        611.3750        606.3750
M.Yes     M      Yes        607.2941        643.2353
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svychisq"
> 
> ### * svychisq
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svytable
> ### Title: Contingency tables for survey data
> ### Aliases: svreptable svytable svytable.svyrep.design
> ###   svytable.survey.design svychisq svychisq.survey.design
> ###   summary.svytable print.summary.svytable summary.svreptable
> ### Keywords: survey category htest
> 
> ### ** Examples
> 
>   data(api)
>   xtabs(~sch.wide+stype, data=apipop)
        stype
sch.wide    E    H    M
     No   472  334  266
     Yes 3949  421  752
> 
>   dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
>   summary(dclus1)
1 - level Cluster Sampling design
With (15) clusters.
svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
Probabilities:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "fpc"      "pw"      
> 
>   svytable(~sch.wide+stype, dclus1)
        stype
sch.wide         E         H         M
     No   406.1640  101.5410  270.7760
     Yes 4467.8035  372.3170  575.3989
>   svychisq(~sch.wide+stype, dclus1)

	Pearson's X^2: Rao & Scott adjustment

data:  svychisq.survey.design(~sch.wide + stype, dclus1) 
F = 5.1934, ndf = 1.495, ddf = 20.925, p-value = 0.02175

>   svychisq(~sch.wide+stype, dclus1, statistic="Chisq")

	Pearson's X^2: Rao & Scott adjustment

data:  svychisq.survey.design(~sch.wide + stype, dclus1, statistic = "Chisq") 
X-squared = 11.9409, df = 2, p-value = 0.005553

>   svychisq(~sch.wide+stype, dclus1, statistic="adjWald")

	Design-based Wald test of association

data:  svychisq.survey.design(~sch.wide + stype, dclus1, statistic = "adjWald") 
F = 2.2296, ndf = 2, ddf = 13, p-value = 0.1471

> 
>   rclus1 <- as.svrepdesign(dclus1)
>   svytable(~sch.wide+stype, rclus1, round=TRUE)
        stype
sch.wide    E    H    M
     No   406  102  271
     Yes 4468  372  575
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svycoxph"
> 
> ### * svycoxph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svycoxph
> ### Title: Survey-weighted Cox models.
> ### Aliases: svycoxph svycoxph.survey.design2 svycoxph.survey.design
> ###   svycoxph.svyrep.design print.svycoxph model.frame.svycoxph
> ###   summary.svycoxph anova.svycoxph extractAIC.svycoxph survfit.svycoxph
> ### Keywords: regression survival survey
> 
> ### ** Examples
> 
> ## Somewhat unrealistic example of nonresponse bias.
> data(pbc, package="survival")
> 
> biasmodel<-glm(I(trt>0)~age*edema,data=pbc)
> pbc$randprob<-fitted(biasmodel)
> 
> dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,trt>0))
> rpbc<-as.svrepdesign(dpbc)
> 
> svycoxph(Surv(time,status)~log(bili)+protime+alb,design=dpbc)
Loading required package: survival
Loading required package: splines

Attaching package: 'survival'


	The following object(s) are masked _by_ .GlobalEnv :

	 pbc 

Stratified Independent Sampling design
svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
    trt > 0))
Call:
svycoxph.survey.design(formula = Surv(time, status) ~ log(bili) + 
    protime + alb, design = dpbc)


            coef exp(coef) se(coef)     z       p
log(bili)  0.882      2.42   0.0931  9.48 0.0e+00
protime    0.324      1.38   0.0827  3.91 9.1e-05
alb       -1.172      0.31   0.2222 -5.28 1.3e-07

Likelihood ratio test=NA  on 3 df, p=NA  n= 312 
> 
> svycoxph(Surv(time,status)~log(bili)+protime+alb,design=rpbc)
Call: as.svrepdesign(dpbc)
Stratified cluster jackknife (JKn) with 312 replicates.
Call:
svycoxph.svyrep.design(formula = Surv(time, status) ~ log(bili) + 
    protime + alb, design = rpbc)


            coef exp(coef) se(coef)     z       p
log(bili)  0.882      2.42    0.102  8.66 0.0e+00
protime    0.324      1.38    0.102  3.16 1.6e-03
alb       -1.172      0.31    0.241 -4.85 1.2e-06

Likelihood ratio test=NA  on 3 df, p=NA  n= 312 
> 
> 
> 
> cleanEx(); ..nameEx <- "svydesign"
> 
> ### * svydesign
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svydesign
> ### Title: Survey sample analysis.
> ### Aliases: svydesign oldsvydesign summary.survey.design
> ###   summary.survey.design2 print.summary.survey.design
> ###   print.summary.survey.design2 print.survey.design print.survey.design2
> ###   [<-.survey.design [.survey.design2 model.frame.survey.design
> ###   na.omit.survey.design na.exclude.survey.design na.fail.survey.design
> ###   dim.survey.design
> ### Keywords: survey univar manip
> 
> ### ** Examples
> 
>   data(api)
> # stratified sample
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> # one-stage cluster sample
> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
> # two-stage cluster sample: weights computed from population sizes.
> dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
> 
> ## multistage sampling has no effect when fpc is not given, so
> ## these are equivalent.
> dclus2wr<-svydesign(id=~dnum+snum, weights=weights(dclus2), data=apiclus2)
> dclus2wr2<-svydesign(id=~dnum, weights=weights(dclus2), data=apiclus2)
> 
> ## syntax for stratified cluster sample
> ##(though the data weren't really sampled this way)
> svydesign(id=~dnum, strata=~stype, weights=~pw, data=apistrat, nest=TRUE)
Stratified 1 - level Cluster Sampling design
With (162) clusters.
svydesign(id = ~dnum, strata = ~stype, weights = ~pw, data = apistrat, 
    nest = TRUE)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyglm"
> 
> ### * svyglm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyglm
> ### Title: Survey-weighted generalised linear models.
> ### Aliases: svyglm svyglm.survey.design svyglm.svyrep.design svrepglm
> ###   print.svyglm summary.svyglm summary.svrepglm print.summary.svyglm
> ###   vcov.svyglm residuals.svyglm residuals.svrepglm coef.svyglm
> ###   extractAIC.svyglm extractAIC.svrepglm logLik.svyglm logLik.svrepglm
> ### Keywords: regression survey
> 
> ### ** Examples
> 
> 
>   data(api)
> 
>   glm(api00~ell+meals+mobility, data=apipop)

Call:  glm(formula = api00 ~ ell + meals + mobility, data = apipop) 

Coefficients:
(Intercept)          ell        meals     mobility  
   833.0947      -0.9915      -2.9180      -0.3233  

Degrees of Freedom: 6189 Total (i.e. Null);  6186 Residual
Null Deviance:	    101800000 
Residual Deviance: 30820000 	AIC: 70270 
> 
>   dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
>   dclus2<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
>   rstrat<-as.svrepdesign(dstrat)
>   rclus2<-as.svrepdesign(dclus2)
> 
>   summary(svyglm(api00~ell+meals+mobility, design=dstrat))

Call:
svyglm.survey.design(formula = api00 ~ ell + meals + mobility, 
    design = dstrat)

Survey design:
svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 820.8873    10.0777  81.456   <2e-16 ***
ell          -0.4806     0.3920  -1.226    0.222    
meals        -3.1415     0.2839 -11.064   <2e-16 ***
mobility      0.2257     0.3932   0.574    0.567    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 5146.106)

Number of Fisher Scoring iterations: 2

>   summary(svyglm(api00~ell+meals+mobility, design=dclus2))

Call:
svyglm.survey.design(formula = api00 ~ ell + meals + mobility, 
    design = dclus2)

Survey design:
svydesign(id = ~dnum + snum, weights = ~pw, data = apiclus2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 811.4907    30.8795  26.279   <2e-16 ***
ell          -2.0592     1.4075  -1.463    0.146    
meals        -1.7772     1.1053  -1.608    0.110    
mobility      0.3253     0.5305   0.613    0.541    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 8296.727)

Number of Fisher Scoring iterations: 2

>   summary(svyglm(api00~ell+meals+mobility, design=rstrat))

Call:
svyglm.svyrep.design(formula = api00 ~ ell + meals + mobility, 
    design = rstrat)

Survey design:
as.svrepdesign(dstrat)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 820.8873    10.5190  78.038   <2e-16 ***
ell          -0.4806     0.4060  -1.184    0.238    
meals        -3.1415     0.2939 -10.691   <2e-16 ***
mobility      0.2257     0.4515   0.500    0.618    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 5146.106)

Number of Fisher Scoring iterations: 2

>   summary(svyglm(api00~ell+meals+mobility, design=rclus2))

Call:
svyglm.svyrep.design(formula = api00 ~ ell + meals + mobility, 
    design = rclus2)

Survey design:
as.svrepdesign(dclus2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 811.4907    37.1474  21.845   <2e-16 ***
ell          -2.0592     1.6177  -1.273    0.205    
meals        -1.7772     1.4901  -1.193    0.235    
mobility      0.3253     0.6307   0.516    0.607    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 8296.727)

Number of Fisher Scoring iterations: 2

> 
>   ## use quasibinomial, quasipoisson to avoid warning messages
>   summary(svyglm(sch.wide~ell+meals+mobility, design=dstrat, family=quasibinomial()))

Call:
svyglm.survey.design(formula = sch.wide ~ ell + meals + mobility, 
    design = dstrat, family = quasibinomial())

Survey design:
svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.835837   0.455551   1.835   0.0681 .
ell         -0.002490   0.013251  -0.188   0.8512  
meals       -0.003152   0.009199  -0.343   0.7322  
mobility     0.060897   0.031925   1.908   0.0579 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for quasibinomial family taken to be 1.103133)

Number of Fisher Scoring iterations: 5

> 
>  
> 
> 
> cleanEx(); ..nameEx <- "svyhist"
> 
> ### * svyhist
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyhist
> ### Title: Histograms
> ### Aliases: svyhist
> ### Keywords: survey hplot
> 
> ### ** Examples
> 
> data(api)
> dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
+     fpc = ~fpc)
> svyhist(~ell, dstrat, main="English language learners",col="peachpuff")
> 
> 
> 
> cleanEx(); ..nameEx <- "svymle"
> 
> ### * svymle
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svymle
> ### Title: Maximum pseudolikelihood estimation in complex surveys
> ### Aliases: svymle print.svymle coef.svymle summary.svymle vcov.svymle
> ### Keywords: survey models optimize
> 
> ### ** Examples
> 
> 
>  data(api)
> 
>  dstrat<-svydesign(id=~1, strata=~stype, weight=~pw, fpc=~fpc, data=apistrat)
> 
>  ## fit with glm
>  m0 <- svyglm(api00~api99+ell,family="gaussian",design=dstrat)
>  ## fit as mle (without gradient)
>  m1 <- svymle(loglike=dnorm,gradient=NULL, design=dstrat, formulas=list(mean=api00~api99+ell, sd=~1),start=list(c(80,1,0),c(20)), log=TRUE)
>  ## with gradient
>  gr<- function(x,mean,sd,log){
+          dm<-2*(x - mean)/(2*sd^2)
+          ds<-(x-mean)^2*(2*(2 * sd))/(2*sd^2)^2 - sqrt(2*pi)/(sd*sqrt(2*pi))
+          cbind(dm,ds)
+       }
>  m2 <- svymle(loglike=dnorm,gradient=gr, design=dstrat, formulas=list(mean=api00~api99+ell, sd=~1),
+                start=list(c(80,1,0),c(20)), log=TRUE, method="BFGS")
> 
>  summary(m0)

Call:
svyglm.survey.design(formula = api00 ~ api99 + ell, design = dstrat, 
    family = "gaussian")

Survey design:
svydesign(id = ~1, strata = ~stype, weight = ~pw, fpc = ~fpc, 
    data = apistrat)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 79.89525   14.61154   5.468 1.37e-07 ***
api99        0.92798    0.01879  49.395  < 2e-16 ***
ell         -0.07312    0.13966  -0.524    0.601    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 718.2254)

Number of Fisher Scoring iterations: 2

>  summary(m1,stderr="model")
Survey-sampled mle: 
svymle(loglike = dnorm, gradient = NULL, design = dstrat, formulas = list(mean = api00 ~ 
    api99 + ell, sd = ~1), start = list(c(80, 1, 0), c(20)), 
    log = TRUE)
                        Coef          SE p.value
mean.(Intercept) 79.86241974 14.29751372  <0.001
mean.api99        0.92800541  0.01949146  <0.001
mean.ell         -0.07270673  0.11779640   0.537
sd.(Intercept)   26.79938270  1.33994303  <0.001
Stratified Independent Sampling design
svydesign(id = ~1, strata = ~stype, weight = ~pw, fpc = ~fpc, 
    data = apistrat)
>  summary(m2)
Survey-sampled mle: 
svymle(loglike = dnorm, gradient = gr, design = dstrat, formulas = list(mean = api00 ~ 
    api99 + ell, sd = ~1), start = list(c(80, 1, 0), c(20)), 
    method = "BFGS", log = TRUE)
                        Coef          SE p.value
mean.(Intercept) 79.99861660 14.61396337  <0.001
mean.api99        0.92783577  0.01878966  <0.001
mean.ell         -0.07381515  0.13967139   0.597
sd.(Intercept)   26.89229311  1.75176987  <0.001
Stratified Independent Sampling design
svydesign(id = ~1, strata = ~stype, weight = ~pw, fpc = ~fpc, 
    data = apistrat)
> 
>  ## More complicated censored data example
>  ## showing that the response variable can be multivariate
> 
>  data(pbc, package="survival")
>  biasmodel<-glm(I(trt>0)~age*edema,data=pbc)
>  pbc$randprob<-fitted(biasmodel)
>  dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,trt>0))
> 
>  lcens<-function(x,mean,sd){
+     ifelse(x[,2]==1,
+            dnorm(log(x[,1]),mean,sd,log=TRUE),
+            pnorm(log(x[,1]),mean,sd,log=TRUE,lower.tail=FALSE)
+            )
+  }
> 
>  gcens<- function(x,mean,sd){
+ 
+         dz<- -dnorm(log(x[,1]),mean,sd)/pnorm(log(x[,1]),mean,sd,lower.tail=FALSE)
+ 
+         dm<-ifelse(x[,2]==1,
+                    2*(log(x[,1]) - mean)/(2*sd^2),
+                    dz*-1/sd)
+         ds<-ifelse(x[,2]==1,
+                    (log(x[,1])-mean)^2*(2*(2 * sd))/(2*sd^2)^2 - sqrt(2*pi)/(sd*sqrt(2*pi)),
+                    ds<- dz*-(log(x[,1])-mean)/(sd*sd))
+         cbind(dm,ds)      
+  }
> 
> svymle(loglike=lcens, gradient=gcens, design=dpbc,
+       formulas=list(mean=I(cbind(time,status))~bili+protime+alb,
+                     sd=~1),
+          start=list(c(10,0,0,0),c(1)))
Survey-sampled mle: 
svymle(loglike = lcens, gradient = gcens, design = dpbc, formulas = list(mean = I(cbind(time, 
    status)) ~ bili + protime + alb, sd = ~1), start = list(c(10, 
    0, 0, 0), c(1)))
Coef:  
mean.(Intercept)        mean.bili     mean.protime         mean.alb 
       8.9803185       -0.0927236       -0.3689851        0.9535612 
  sd.(Intercept) 
       0.9888172 
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyplot"
> 
> ### * svyplot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyplot
> ### Title: Plots for survey data
> ### Aliases: svyplot
> ### Keywords: survey hplot
> 
> ### ** Examples
> 
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> 
> svyplot(api00~api99, design=dstrat, style="bubble")
> ## Not run: 
> ##D ## these two require the hexbin package from Bioconductor
> ##D svyplot(api00~api99, design=dstrat, style="hex", xlab="1999 API",ylab="2000 API")
> ##D svyplot(api00~api99, design=dstrat, style="grayhex",legend=0)
> ## End(Not run)
> ## Subsampling doesn't really make sense for such a small survey
> svyplot(api00~api99, design=dstrat, style="subsample")
> svyplot(api00~stype, design=dstrat, style="subsample")
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyquantile"
> 
> ### * svyquantile
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyquantile
> ### Title: Quantiles for sample surveys
> ### Aliases: svyquantile print.svyquantile SE.svyquantile
> ###   svyquantile.survey.design svyquantile.svyrep.design svrepquantile
> ### Keywords: univar survey
> 
> ### ** Examples
> 
> 
>   data(api)
>   ## population
>   quantile(apipop$api00,c(.25,.5,.75))
25% 50% 75% 
565 667 761 
> 
>   ## one-stage cluster sample
>   dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
>   svyquantile(~api00, dclus1, c(.25,.5,.75),ci=TRUE)
$quantiles
        0.25 0.5  0.75
api00 551.75 652 717.5

$CIs
, , api00

           0.25      0.5     0.75
(lower 493.2835 564.3250 696.0000
upper) 622.6495 710.8375 761.1355


> 
>   dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
>   (qapi<-svyquantile(~api00, dclus1, c(.25,.5,.75),ci=TRUE, interval.type="score"))
$quantiles
        0.25 0.5  0.75
api00 551.75 652 717.5

$CIs
, , api00

           0.25      0.5     0.75
(lower 514.9998 581.9996 669.0003
upper) 632.0005 699.9997 749.0001


>   SE(qapi)
    0.25      0.5     0.75 
29.84766 30.10264 20.40850 
> 
>   #stratified sample
>   dstrat<-svydesign(id=~1, strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
>   svyquantile(~api00, dstrat, c(.25,.5,.75),ci=TRUE)
$quantiles
          0.25      0.5     0.75
api00 562.2056 667.2358 755.1226

$CIs
, , api00

           0.25      0.5     0.75
(lower 534.0000 636.4618 725.1052
upper) 594.2887 681.0000 776.9814


> 
>   #stratified sample, replicate weights
>   # interval="probability" is necessary for jackknife weights
>   rstrat<-as.svrepdesign(dstrat)
>   svyquantile(~api00, rstrat, c(.25,.5,.75), interval="probability")
Statistic:
       api00
q0.25 562.21
q0.5  667.24
q0.75 755.12
SE:
         api00
q0.25 12.29485
q0.5  11.63687
q0.75 15.70636
> 
>   # BRR method
>   data(scd)
>   repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+               c(0,1,0,1,1,0))
>   scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
>   svyquantile(~arrests+alive, design=scdrep, quantile=0.5, interval="quantile")
Statistic:
     arrests alive
q0.5  265.75  43.5
SE:
      arrests    alive
q0.5 77.22896 5.396758
> 
>  
> 
> 
> cleanEx(); ..nameEx <- "svyratio"
> 
> ### * svyratio
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyratio
> ### Title: Ratio estimation
> ### Aliases: svyratio svrepratio print.svyratio print.svyratio_separate
> ###   svyratio.survey.design svyratio.survey.design2 svyratio.svyrep.design
> ###   SE.svyratio predict.svyratio predict.svyratio_separate
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> 
> ## survey design objects
> scddes<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE, fpc=rep(5,6))
> scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+ nest=TRUE)
> 
> # convert to BRR replicate weights
> scd2brr <- as.svrepdesign(scdnofpc, type="BRR")
> 
> # use BRR replicate weights from Levy and Lemeshow
> repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+ c(0,1,0,1,1,0))
> scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
> 
> # ratio estimates
> svyratio(~alive, ~arrests, design=scddes)
Ratio estimator: svyratio.survey.design2(~alive, ~arrests, design = scddes)
Ratios=
        arrests
alive 0.1535064
SEs=
          arrests
alive 0.007596705
> svyratio(~alive, ~arrests, design=scdnofpc)
Ratio estimator: svyratio.survey.design2(~alive, ~arrests, design = scdnofpc)
Ratios=
        arrests
alive 0.1535064
SEs=
          arrests
alive 0.009807304
> svyratio(~alive, ~arrests, design=scd2brr)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scd2brr)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> svyratio(~alive, ~arrests, design=scdrep)
Ratio estimator: svyratio.svyrep.design(~alive, ~arrests, design = scdrep)
Ratios=
        arrests
alive 0.1535064
SEs=
          [,1]
[1,] 0.0094184
> 
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> 
> ## domain means are ratio estimates, but available directly
> svyratio(~I(api.stu*(comp.imp=="Yes")), ~as.numeric(comp.imp=="Yes"), dstrat)
Ratio estimator: svyratio.survey.design2(~I(api.stu * (comp.imp == "Yes")), ~as.numeric(comp.imp == 
    "Yes"), dstrat)
Ratios=
                                 as.numeric(comp.imp == "Yes")
I(api.stu * (comp.imp == "Yes"))                      439.9305
SEs=
                                 as.numeric(comp.imp == "Yes")
I(api.stu * (comp.imp == "Yes"))                      19.24367
> svymean(~api.stu, subset(dstrat, comp.imp=="Yes"))
          mean     SE
api.stu 439.93 19.244
> 
> ## separate and combined ratio estimates of total
> (sep<-svyratio(~api.stu,~enroll, dstrat,separate=TRUE))
Stratified ratio estimate: svyratio.survey.design2(~api.stu, ~enroll, dstrat, separate = TRUE)
Ratio estimator: Stratum == "E"
Ratios=
           enroll
api.stu 0.8518163
SEs=
            enroll
api.stu 0.00703236
Ratio estimator: Stratum == "H"
Ratios=
           enroll
api.stu 0.8105702
SEs=
            enroll
api.stu 0.02047726
Ratio estimator: Stratum == "M"
Ratios=
           enroll
api.stu 0.8356958
SEs=
            enroll
api.stu 0.01818744
> (com<-svyratio(~api.stu, ~enroll, dstrat))
Ratio estimator: svyratio.survey.design2(~api.stu, ~enroll, dstrat)
Ratios=
           enroll
api.stu 0.8369569
SEs=
             enroll
api.stu 0.007757103
> 
> stratum.totals<-list(E=1877350, H=1013824, M=920298)
> 
> predict(sep, total=stratum.totals)
$total
         enroll
api.stu 3190022

$se
          enroll
api.stu 29756.44

> predict(com, total=sum(unlist(stratum.totals)))
$total
         enroll
api.stu 3190038

$se
          enroll
api.stu 29565.98

> 
> 
> 
> 
> cleanEx(); ..nameEx <- "svyrecvar"
> 
> ### * svyrecvar
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: svyrecvar
> ### Title: Variance estimation for multistage surveys
> ### Aliases: svyrecvar multistage
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(mu284)
> dmu284<-svydesign(id=~id1+id2,fpc=~n1+n2, data=mu284)
> svytotal(~y1, dmu284)
   total     SE
y1 15080 2274.3
> 
> data(api)
> # two-stage cluster sample
> dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
> summary(dclus2)
2 - level Cluster Sampling design
With (40, 126) clusters.
svydesign(id = ~dnum + snum, fpc = ~fpc1 + fpc2, data = apiclus2)
Probabilities:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.003669 0.037740 0.052840 0.042390 0.052840 0.052840 
Population size (PSUs): 757 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc1"     "fpc2"    
> svymean(~api00, dclus2)
        mean     SE
api00 670.81 30.099
> svytotal(~enroll, dclus2,na.rm=TRUE)
         total     SE
enroll 2639273 799638
> 
> # two-stage `with replacement'
> dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
> summary(dclus2wr)
2 - level Cluster Sampling design
With (40, 126) clusters.
svydesign(id = ~dnum + snum, weights = ~pw, data = apiclus2)
Probabilities:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.003669 0.037740 0.052840 0.042390 0.052840 0.052840 
Data variables:
 [1] "cds"      "stype"    "name"     "sname"    "snum"     "dname"   
 [7] "dnum"     "cname"    "cnum"     "flag"     "pcttest"  "api00"   
[13] "api99"    "target"   "growth"   "sch.wide" "comp.imp" "both"    
[19] "awards"   "meals"    "ell"      "yr.rnd"   "mobility" "acs.k3"  
[25] "acs.46"   "acs.core" "pct.resp" "not.hsg"  "hsg"      "some.col"
[31] "col.grad" "grad.sch" "avg.ed"   "full"     "emer"     "enroll"  
[37] "api.stu"  "pw"       "fpc1"     "fpc2"    
> svymean(~api00, dclus2wr)
        mean     SE
api00 670.81 30.712
> svytotal(~enroll, dclus2wr,na.rm=TRUE)
         total     SE
enroll 2639273 820261
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "update.survey.design"
> 
> ### * update.survey.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: update.survey.design
> ### Title: Add variables to a survey design
> ### Aliases: update.survey.design update.svyrep.design
> ### Keywords: survey manip
> 
> ### ** Examples
> 
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat,
+ fpc=~fpc)
> dstrat<-update(dstrat, apidiff=api00-api99)
> svymean(~api99+api00+apidiff, dstrat)
           mean     SE
api99   629.395 9.9639
api00   662.287 9.4089
apidiff  32.893 2.0511
> 
> 
> 
> cleanEx(); ..nameEx <- "weights.survey.design"
> 
> ### * weights.survey.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: weights.survey.design
> ### Title: Survey design weights
> ### Aliases: weights.survey.design weights.svyrep.design weights.survey_fpc
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> 
> scddes<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
+                  nest=TRUE, fpc=rep(5,6))
> repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1), c(0,1,0,1,1,0))
> scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
> 
> weights(scdrep)
NULL
> weights(scdrep, type="sampling")
[1] 1 1 1 1 1 1
> weights(scdrep, type="analysis")
     [,1] [,2] [,3] [,4]
[1,]    2    2    0    0
[2,]    0    0    2    2
[3,]    2    0    2    0
[4,]    0    2    0    2
[5,]    2    0    0    2
[6,]    0    2    2    0
> weights(scddes)
1 2 3 4 5 6 
1 1 1 1 1 1 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "withReplicates"
> 
> ### * withReplicates
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: withReplicates
> ### Title: Compute variances by replicate weighting
> ### Aliases: withReplicates
> ### Keywords: survey
> 
> ### ** Examples
> 
> data(scd)
> repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
+ c(0,1,0,1,1,0))
> scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)
Warning in svrepdesign(data = scd, type = "BRR", repweights = repweights) : 
	 No sampling weights provided: equal probability assumed
> 
> a<-svyratio(~alive, ~arrests, design=scdrep)
> print(a$ratio)
        arrests
alive 0.1535064
> print(a$var)
             [,1]
[1,] 8.870627e-05
> withReplicates(scdrep, quote(sum(.weights*alive)/sum(.weights*arrests)))
       theta     SE
[1,] 0.15351 0.0094
> withReplicates(scdrep, function(w,data) sum(w*data$alive)/sum(w*data$arrests))
       theta     SE
[1,] 0.15351 0.0094
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  25.27 1.09 26.51 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')