R version 3.1.0 (2014-04-10) -- "Spring Dance"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.1.0 (64-bit)
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)
Attaching package: 'survey'
The following object is masked from 'package:graphics':
dotchart
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
>
>
> ## regression estimator of total, three ways
> pop<-data.frame(enroll=sum(apipop$enroll, na.rm=TRUE))
> npop <- sum(!is.na(apipop$enroll))
>
> api.reg <- svyglm(api.stu~enroll, design=dstrat)
> a <- predict(api.reg, newdata=pop, total=npop)
> b <- svytotal(~api.stu, calibrate(dstrat, ~enroll, pop=c(npop, pop$enroll)))
>
> all.equal(as.vector(coef(a)),as.vector(coef(b)))
[1] TRUE
> all.equal(as.vector(SE(a)), as.vector(SE(b)))
[1] TRUE
> if(!is.null(getOption("DEBUG"))){ ## uses 6194x6194 matrix
+ d <- predict(api.reg, newdata=na.omit(apipop[,"enroll",drop=FALSE]))
+ all.equal(as.vector(coef(a)), as.vector(sum(coef(d))))
+ all.equal(as.vector(SE(a)), as.vector(sqrt(sum(vcov(d)))))
+ }
>
> ## classical ratio estimator, four ways.
> api.reg2 <- svyglm(api.stu~enroll-1, design=dstrat,
+ family=quasi(link="identity", var="mu"))
>
> a <- predict(api.reg2, newdata=pop, total=npop)
> b <- svytotal(~api.stu,
+ calibrate(dstrat, ~enroll-1, pop= pop$enroll, variance=2))
> e <- predict(svyratio(~api.stu, ~enroll, dstrat),total=pop$enroll)
>
> all.equal(as.vector(coef(a)),as.vector(coef(b)))
[1] TRUE
> all.equal(as.vector(SE(a)), as.vector(SE(b)))
[1] TRUE
> all.equal(as.vector(coef(a)),as.vector(e$total))
[1] TRUE
> all.equal(as.vector(SE(a)), as.vector(e$se))
[1] TRUE
> if(!is.null(getOption("DEBUG"))){## uses 6194x6194 matrix
+ d <- predict(api.reg2, newdata=na.omit(apipop[,"enroll",drop=FALSE]))
+ all.equal(as.vector(coef(a)), as.vector(sum(coef(d))))
+ all.equal(as.vector(SE(a)), as.vector(sqrt(sum(vcov(d)))))
+ }
>
> proc.time()
user system elapsed
0.239 0.026 0.272
Generated by dwww version 1.15 on Thu May 23 17:25:42 CEST 2024.