R Under development (unstable) (2024-03-07 r86063 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)
> 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 
    1.1     0.1     1.2