Running '3stage2phase.R' [3s]
  Running 'DBIcheck.R' [2s]
  Running 'anova-svycoxph.R' [156s]
  Running 'api.R' [172s]
  Running 'badcal.R' [2s]
  Comparing 'badcal.Rout' to 'badcal.Rout.save' ... OK
  Running 'brewer_cpp.R' [173s]
  Running 'bycovmat.R' [159s]
  Running 'caleg.R' [172s]
  Running 'check.R' [2s]
  Comparing 'check.Rout' to 'check.Rout.save' ... OK
  Running 'confintrep.R' [172s]
  Running 'contrast-replicates.R' [176s]
  Running 'coxph-termtest.R' [149s]
  Running 'deff.R' [2s]
  Comparing 'deff.Rout' to 'deff.Rout.save' ... OK
  Running 'defftest.R' [155s]
  Running 'degf-svrepdesign.R' [2s]
  Running 'domain.R' [2s]
  Comparing 'domain.Rout' to 'domain.Rout.save' ... OK
  Running 'fpc.R' [2s]
  Running 'glm-scoping.R' [152s]
  Running 'kalton.R' [145s]
  Running 'logranktest.R' [2s]
  Running 'lonely.psu.R' [161s]
  Running 'mtcars-var.R' [2s]
  Running 'multistage-rcpp.R' [167s]
  Running 'multistage.R' [2s]
  Running 'na_action.R' [174s]
  Running 'newquantile.R' [2s]
  Running 'nwts-cch.R' [3s]
  Running 'nwts.R' [167s]
  Running 'poisson.R' [2s]
  Running 'pps.R' [2s]
  Running 'qrule-swiss.R' [1s]
  Running 'quantile.R' [162s]
  Running 'quantiles-chile.R' [2s]
  Running 'rakecheck.R' [2s]
  Comparing 'rakecheck.Rout' to 'rakecheck.Rout.save' ... OK
  Running 'raowuboot.R' [157s]
  Running 'regTermTest-missing.R' [2s]
  Running 'regpredict.R' [154s]
  Running 'rss_scores.R' [2s]
  Running 'scoping.R' [2s]
  Running 'survcurve.R' [4s]
  Comparing 'survcurve.Rout' to 'survcurve.Rout.save' ... OK
  Running 'svyby-strings.R' [2s]
  Running 'svyby_bug.R' [150s]
  Running 'svyby_se.R' [175s]
  Running 'svycontrast.R' [2s]
  Running 'svyivreg-var.R' [2s]
  Running 'svyivreg.R' [165s]
  Running 'svyolr-rake-subset.R' [2s]
  Running 'svyolr.R' [2s]
  Running 'toy_example_for_postStratify.R' [157s]
  Running 'twophase.R' [167s]
Running the tests in 'tests/anova-svycoxph.R' failed.
Complete output:
  > 
  >  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
  
  > example(svycoxph, ask=FALSE)
  
  svycxp> ## Somewhat unrealistic example of nonresponse bias.
  svycxp> data(pbc, package="survival")
  
  svycxp> pbc$randomized<-with(pbc, !is.na(trt) & trt>0)
  
  svycxp> biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial)
  
  svycxp> pbc$randprob<-fitted(biasmodel)
  
  svycxp> if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0
  
  svycxp> dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
  
  svycxp> rpbc<-as.svrepdesign(dpbc)
  
  svycxp> (model<-svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=dpbc))
  Call:
  svycoxph(formula = Surv(time, status > 0) ~ log(bili) + protime + 
      albumin, design = dpbc)
  
                coef exp(coef) se(coef) robust se      z        p
  log(bili)  0.88592   2.42522  0.09140   0.09048  9.791  < 2e-16
  protime    0.24487   1.27745  0.07825   0.08122  3.015  0.00257
  albumin   -1.04298   0.35240  0.21211   0.20454 -5.099 3.41e-07
  
  Likelihood ratio test=  on 3 df, p=
  n= 312, number of events= 144 
  
  svycxp> svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=rpbc)
  Call:
  svycoxph.svyrep.design(formula = Surv(time, status > 0) ~ log(bili) + 
      protime + albumin, design = rpbc)
  
                coef exp(coef) se(coef)      z        p
  log(bili)  0.88592   2.42522  0.09838  9.005  < 2e-16
  protime    0.24487   1.27745  0.09373  2.612  0.00899
  albumin   -1.04298   0.35240  0.21966 -4.748 2.05e-06
  
  Likelihood ratio test=NA  on 3 df, p=NA
  n= 312, number of events= 144 
  
  svycxp> s<-predict(model,se=TRUE, type="curve",
  svycxp+      newdata=data.frame(bili=c(3,9), protime=c(10,10), albumin=c(3.5,3.5)))
  
  svycxp> plot(s[[1]],ci=TRUE,col="sienna")
  
  svycxp> lines(s[[2]], ci=TRUE,col="royalblue")
  
  svycxp> quantile(s[[1]], ci=TRUE)
  0.75  0.5 0.25 
  1435 2503 3574 
  attr(,"ci")
       0.025 0.975
  0.75  1217  1786
  0.5   2256  3170
  0.25  3222   Inf
  
  svycxp> confint(s[[2]], parm=365*(1:5))
           0.025     0.975
  365  0.8375139 0.9453781
  730  0.7382750 0.8999016
  1095 0.4784105 0.7478460
  1460 0.3192009 0.6206764
  1825 0.2149475 0.5292978
  > m<-update(model, .~.+I(protime^2))
  > a<-anova(m,model)
  > b<-anova(m, model,force=TRUE)
  > stopifnot(isTRUE(all.equal(b[2:6],a[c(3,4,6,7,8)])))
  > 
  > proc.time()
     user  system elapsed 
     2.26    0.25    2.46 
Running the tests in 'tests/api.R' failed.
Complete output:
  > 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
  
  > options(survey.replicates.mse=TRUE)
  > example(api)
  
  api> library(survey)
  
  api> data(api)
  
  api> mean(apipop$api00)
  [1] 664.7126
  
  api> sum(apipop$enroll, na.rm=TRUE)
  [1] 3811472
  
  api> #stratified sample
  api> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
  
  api> summary(dstrat)
  Stratified Independent Sampling design
  dstrat<-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    H    M 
  4421  755 1018 
  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"     
  
  api> svymean(~api00, dstrat)
          mean     SE
  api00 662.29 9.4089
  
  api> svytotal(~enroll, dstrat, na.rm=TRUE)
           total     SE
  enroll 3687178 114642
  
  api> # one-stage cluster sample
  api> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  
  api> summary(dclus1)
  1 - level Cluster Sampling design
  With (15) clusters.
  dclus1<-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"      
  
  api> svymean(~api00, dclus1)
          mean     SE
  api00 644.17 23.542
  
  api> svytotal(~enroll, dclus1, na.rm=TRUE)
           total     SE
  enroll 3404940 932235
  
  api> # two-stage cluster sample
  api> dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
  
  api> summary(dclus2)
  2 - level Cluster Sampling design
  With (40, 126) clusters.
  dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
  Probabilities:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003669 0.037743 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"    
  
  api> svymean(~api00, dclus2)
          mean     SE
  api00 670.81 30.099
  
  api> svytotal(~enroll, dclus2, na.rm=TRUE)
           total     SE
  enroll 2639273 799638
  
  api> # two-stage `with replacement'
  api> dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
  
  api> summary(dclus2wr)
  2 - level Cluster Sampling design (with replacement)
  With (40, 126) clusters.
  dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
  Probabilities:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003669 0.037743 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"    
  
  api> svymean(~api00, dclus2wr)
          mean     SE
  api00 670.81 30.712
  
  api> svytotal(~enroll, dclus2wr, na.rm=TRUE)
           total     SE
  enroll 2639273 820261
  
  api> # convert to replicate weights
  api> rclus1<-as.svrepdesign(dclus1)
  
  api> summary(rclus1)
  Call: as.svrepdesign.default(dclus1)
  Unstratified cluster jacknife (JK1) with 15 replicates and MSE variances.
  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"      
  
  api> svymean(~api00, rclus1)
          mean     SE
  api00 644.17 26.335
  
  api> svytotal(~enroll, rclus1, na.rm=TRUE)
           total     SE
  enroll 3404940 932235
  
  api> # post-stratify on school type
  api> pop.types<-xtabs(~stype, data=apipop)
  
  api> rclus1p<-postStratify(rclus1, ~stype, pop.types)
  
  api> dclus1p<-postStratify(dclus1, ~stype, pop.types)
  
  api> summary(dclus1p)
  1 - level Cluster Sampling design
  With (15) clusters.
  dclus1p<-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"      
  
  api> summary(rclus1p)
  Call: rclus1p<-postStratify(rclus1, ~stype, pop.types)
  Unstratified cluster jacknife (JK1) with 15 replicates and MSE variances.
  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"      
  
  api> svymean(~api00, dclus1p)
          mean     SE
  api00 642.31 23.921
  
  api> svytotal(~enroll, dclus1p, na.rm=TRUE)
           total     SE
  enroll 3680893 406293
  
  api> svymean(~api00, rclus1p)
          mean     SE
  api00 642.31 26.936
  
  api> svytotal(~enroll, rclus1p, na.rm=TRUE)
           total     SE
  enroll 3680893 473434
  > 
  > options(survey.replicates.mse=FALSE)
  > example(api)
  
  api> library(survey)
  
  api> data(api)
  
  api> mean(apipop$api00)
  [1] 664.7126
  
  api> sum(apipop$enroll, na.rm=TRUE)
  [1] 3811472
  
  api> #stratified sample
  api> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
  
  api> summary(dstrat)
  Stratified Independent Sampling design
  dstrat<-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    H    M 
  4421  755 1018 
  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"     
  
  api> svymean(~api00, dstrat)
          mean     SE
  api00 662.29 9.4089
  
  api> svytotal(~enroll, dstrat, na.rm=TRUE)
           total     SE
  enroll 3687178 114642
  
  api> # one-stage cluster sample
  api> dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  
  api> summary(dclus1)
  1 - level Cluster Sampling design
  With (15) clusters.
  dclus1<-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"      
  
  api> svymean(~api00, dclus1)
          mean     SE
  api00 644.17 23.542
  
  api> svytotal(~enroll, dclus1, na.rm=TRUE)
           total     SE
  enroll 3404940 932235
  
  api> # two-stage cluster sample
  api> dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
  
  api> summary(dclus2)
  2 - level Cluster Sampling design
  With (40, 126) clusters.
  dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
  Probabilities:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003669 0.037743 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"    
  
  api> svymean(~api00, dclus2)
          mean     SE
  api00 670.81 30.099
  
  api> svytotal(~enroll, dclus2, na.rm=TRUE)
           total     SE
  enroll 2639273 799638
  
  api> # two-stage `with replacement'
  api> dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
  
  api> summary(dclus2wr)
  2 - level Cluster Sampling design (with replacement)
  With (40, 126) clusters.
  dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
  Probabilities:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003669 0.037743 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"    
  
  api> svymean(~api00, dclus2wr)
          mean     SE
  api00 670.81 30.712
  
  api> svytotal(~enroll, dclus2wr, na.rm=TRUE)
           total     SE
  enroll 2639273 820261
  
  api> # convert to replicate weights
  api> rclus1<-as.svrepdesign(dclus1)
  
  api> summary(rclus1)
  Call: as.svrepdesign.default(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"      
  
  api> svymean(~api00, rclus1)
          mean     SE
  api00 644.17 26.329
  
  api> svytotal(~enroll, rclus1, na.rm=TRUE)
           total     SE
  enroll 3404940 932235
  
  api> # post-stratify on school type
  api> pop.types<-xtabs(~stype, data=apipop)
  
  api> rclus1p<-postStratify(rclus1, ~stype, pop.types)
  
  api> dclus1p<-postStratify(dclus1, ~stype, pop.types)
  
  api> summary(dclus1p)
  1 - level Cluster Sampling design
  With (15) clusters.
  dclus1p<-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"      
  
  api> summary(rclus1p)
  Call: rclus1p<-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"      
  
  api> svymean(~api00, dclus1p)
          mean     SE
  api00 642.31 23.921
  
  api> svytotal(~enroll, dclus1p, na.rm=TRUE)
           total     SE
  enroll 3680893 406293
  
  api> svymean(~api00, rclus1p)
          mean     SE
  api00 642.31 26.934
  
  api> svytotal(~enroll, rclus1p, na.rm=TRUE)
           total     SE
  enroll 3680893 473431
  > 
  > proc.time()
     user  system elapsed 
     1.89    0.25    2.14 
Running the tests in 'tests/brewer_cpp.R' failed.
Complete output:
  > ## pps="brewer" can't use rcpp
  > ## this checks that it doesn't
  > 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(election)
  > 
  > dpps_br<- svydesign(id=~1,  fpc=~p, data=election_pps, pps="brewer")
  > options(survey.use_rcpp=TRUE)
  > a<-svytotal(~Bush+Kerry+Nader, dpps_br)
  > options(survey.use_rcpp=FALSE)
  > b<-svytotal(~Bush+Kerry+Nader, dpps_br)
  > 
  > stopifnot(identical(a,b))
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.20    1.56 
Running the tests in 'tests/bycovmat.R' failed.
Complete output:
  > 
  > 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)
  > options(survey.replicates.mse=TRUE)
  > dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > rclus1<-as.svrepdesign(dclus1)
  > 
  > a<-svyby(~api00+api99, ~comp.imp+sch.wide,design=rclus1,svymean,
  +          covmat=TRUE,drop.empty.groups=FALSE)
  > b<-svyby(~api00+api99, ~comp.imp+sch.wide,design=rclus1,svymean,
  +          covmat=TRUE,drop.empty.groups=TRUE)
  > 
  > stopifnot(all.equal(
  +               as.vector(as.matrix(SE(a))),as.vector(sqrt(diag(vcov(a))))
  + ))
  > stopifnot(all.equal(
  +               as.vector(as.matrix(SE(b))),as.vector(sqrt(diag(vcov(b))))
  +               ))
  > 
  > rat <- svyratio(~ell+mobility, ~mobility+meals, dclus1,covmat=TRUE)
  > all <- svytotal(~ell+mobility+meals, dclus1)
  > 
  > stopifnot(all(abs(vcov(svycontrast(all,
  +                                    list(quote(ell/mobility),
  +                                         quote(mobility/mobility),
  +                                         quote(ell/meals),quote(mobility/meals))))
  +                   -vcov(rat))<1e-10))
  > 
  > stopifnot(all(abs(SE(rat)-sqrt(diag(vcov(rat))))<1e-10))
  > 
  > rat <- svyratio(~ell+mobility, ~mobility+meals, rclus1,covmat=TRUE)
  > all <- svytotal(~ell+mobility+meals, rclus1, return.replicates=TRUE)
  > 
  > con<-svycontrast(all,
  +                  list(quote(ell/mobility),
  +                       quote(mobility/mobility),
  +                       quote(ell/meals),quote(mobility/meals)))
  > 
  > stopifnot(all(abs(survey:::svrVar(con$replicates, rclus1$scale,rclus1$rscales,mse=rclus1$mse, coef=coef(con))-vcov(rat))<1e-10))
  > 
  > options(survey.replicates.mse=FALSE)
  > dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > rclus1<-as.svrepdesign(dclus1)
  > 
  > a<-svyby(~api00+api99, ~comp.imp+sch.wide,design=rclus1,svymean,
  +          covmat=TRUE,drop.empty.groups=FALSE)
  > b<-svyby(~api00+api99, ~comp.imp+sch.wide,design=rclus1,svymean,
  +          covmat=TRUE,drop.empty.groups=TRUE)
  > 
  > stopifnot(all.equal(
  +               as.vector(as.matrix(SE(a))), as.vector(sqrt(diag(vcov(a))))
  + ))
  > stopifnot(all.equal(
  +               as.vector(as.matrix(SE(b))), as.vector(sqrt(diag(vcov(b))))
  +               ))
  > 
  > rat <- svyratio(~ell+mobility, ~mobility+meals, dclus1,covmat=TRUE)
  > all <- svytotal(~ell+mobility+meals, dclus1)
  > 
  > stopifnot(all(abs(vcov(svycontrast(all,
  +                                    list(quote(ell/mobility),
  +                                         quote(mobility/mobility),
  +                                         quote(ell/meals),quote(mobility/meals))))
  +                   -vcov(rat))<1e-10))
  > 
  > stopifnot(all(abs(SE(rat)-sqrt(diag(vcov(rat))))<1e-10))
  > 
  > rat <- svyratio(~ell+mobility, ~mobility+meals, rclus1,covmat=TRUE)
  > all <- svytotal(~ell+mobility+meals, rclus1, return.replicates=TRUE)
  > 
  > con<-svycontrast(all,
  +                  list(quote(ell/mobility),
  +                       quote(mobility/mobility),
  +                       quote(ell/meals),quote(mobility/meals)))
  > 
  > stopifnot(all(abs(survey:::svrVar(con$replicates, rclus1$scale,rclus1$rscales,mse=rclus1$mse, coef=coef(con))-vcov(rat))<1e-10))
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.56    0.25    1.79 
Running the tests in 'tests/caleg.R' failed.
Complete output:
  > ##
  > ## Calibration examples
  > ##
  > 
  > 
  > ## Example of calibration to first-stage clusters
  > 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)
  > 
  > clusters<-table(apiclus2$dnum)
  > clusters<-clusters[clusters>1 & names(clusters)!="639"]
  > apiclus2a<-subset(apiclus2, dnum %in% as.numeric(names(clusters)))
  > 
  > dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2a)
  > 
  > popclusters<-subset(apipop, dnum %in% as.numeric(names(clusters)))
  > 
  > pop<-lapply(as.numeric(names(clusters)), function(cluster) {
  +   colSums(model.matrix(~api99, model.frame(~api99, subset(popclusters, dnum %in% cluster))))})
  > 
  > names(pop)<-names(clusters)
  > 
  > dclus2g<-calibrate(dclus2, ~api99, pop,stage=1)
  > 
  > svymean(~api99, dclus2)
          mean     SE
  api99 642.14 31.434
  > svymean(~api99, dclus2g)
          mean    SE
  api99 654.49 29.82
  > 
  > round(svyby(~api99, ~dnum, design=dclus2, svymean),4)
      dnum    api99      se
  83    83 694.3333  0.0000
  132  132 505.0000  0.0000
  152  152 574.0000  0.0000
  173  173 894.7500  0.0000
  198  198 533.7500  0.0000
  200  200 589.8000  6.8335
  228  228 477.0000  0.0000
  295  295 646.4000  0.0000
  302  302 903.5000  0.0000
  403  403 852.4000  0.0000
  452  452 533.0000  0.0000
  480  480 614.2000  0.0000
  523  523 580.5000  0.0000
  534  534 564.6000  0.0000
  549  549 896.2000  0.0000
  552  552 730.0000  0.0000
  570  570 518.4000  7.5478
  575  575 800.8000  4.2513
  596  596 785.6000  2.4155
  620  620 591.6000 10.5869
  638  638 560.2000  4.0954
  674  674 760.0000  0.0000
  679  679 610.2500  0.0000
  687  687 718.6667  0.0000
  701  701 651.5000  0.0000
  711  711 690.5000  0.0000
  731  731 702.0000  2.1744
  768  768 562.5000  0.0000
  781  781 854.4000  0.7456
  > 
  > round(svyby(~api99, ~dnum, design=dclus2g, svymean),4)
      dnum    api99 se
  83    83 694.3333  0
  132  132 505.0000  0
  152  152 574.0000  0
  173  173 894.7500  0
  198  198 533.7500  0
  200  200 567.5455  0
  228  228 477.0000  0
  295  295 646.4000  0
  302  302 903.5000  0
  403  403 852.4000  0
  452  452 533.0000  0
  480  480 614.2000  0
  523  523 580.5000  0
  534  534 564.6000  0
  549  549 896.2000  0
  552  552 730.0000  0
  570  570 548.9444  0
  575  575 824.5357  0
  596  596 787.5714  0
  620  620 609.3750  0
  638  638 585.6429  0
  674  674 760.0000  0
  679  679 610.2500  0
  687  687 718.6667  0
  701  701 651.5000  0
  711  711 690.5000  0
  731  731 700.6667  0
  768  768 562.5000  0
  781  781 851.0000  0
  > 
  > ## Averaging to first stage
  > 
  > dclus1<- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
  > pop<-colSums(cbind(1,apipop$enroll),na.rm=TRUE)
  > 
  > dclus1g<-calibrate(dclus1, ~enroll, pop, aggregate=1)
  > 
  > svytotal(~enroll,dclus1g)
           total SE
  enroll 3811472  0
  > svytotal(~api.stu,dclus1g)
            total    SE
  api.stu 3242857 38967
  > 
  > #variation within clusters should be zero
  > all.equal(0, max(ave(weights(dclus1g),dclus1g$cluster,FUN=var),na.rm=TRUE))
  [1] TRUE
  > 
  > ##bounded weights
  >  dclus1g<-calibrate(dclus1, ~enroll, pop)
  >  range(weights(dclus1g)/weights(dclus1))
  [1] 0.7906782 1.7891164
  >  dclus1gb<-calibrate(dclus1, ~enroll, pop, bounds=c(.6,1.5))
  >  range(weights(dclus1gb)/weights(dclus1))
  [1] 0.7198751 1.5000000
  > 
  > ## Ratio estimators
  > dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
  > svytotal(~api.stu,dstrat)
            total    SE
  api.stu 3086009 99477
  > common<-svyratio(~api.stu, ~enroll, dstrat, separate=FALSE)
  > total.enroll<-sum(apipop$enroll,na.rm=TRUE)
  > predict(common, total=total.enroll)
  $total
           enroll
  api.stu 3190038
  
  $se
            enroll
  api.stu 29565.98
  
  > dstratg<-calibrate(dstrat,~enroll-1, total.enroll, variance=1)
  > svytotal(~api.stu, dstratg)
            total    SE
  api.stu 3190038 29566
  > 
  > ## postStratify vs calibrate in stratified sample (Ben French)
  > set.seed(17)
  > dat<-data.frame(y=rep(0:1,each=100),x=rnorm(200)+2*rep(0:1,each=100),
  +                 z=rbinom(200,1,.2), fpc=rep(c(100,10000),each=100))
  > dat$w<-ifelse(dat$y,dat$z,1-dat$z)
  > popw<-data.frame(w=c("0","1"), Freq=c(2000,8000))
  >  des<-svydesign(id=~1,fpc=~fpc, data=dat,strata=~y)
  > postStratify(des,~w,popw)->dps
  > dcal<-calibrate(des,~factor(w), pop=c(10000,8000))
  > 
  > all.equal(SE(svymean(~x,dcal)),SE(svymean(~x,dps)))
  [1] TRUE
  > 
  > ## missing data in calibrated design
  > dps$variables$z[1]<-NA
  > summary(svyglm(y~z+x,design=dps,family=quasibinomial))
  
  Call:
  svyglm(formula = y ~ z + x, design = dps, family = quasibinomial)
  
  Survey design:
  postStratify(des, ~w, popw)
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  -0.1203     0.3380  -0.356    0.722    
  z             6.2118     0.6451   9.630   <2e-16 ***
  x             2.2602     0.2514   8.992   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.919987)
  
  Number of Fisher Scoring iterations: 9
  
  > 
  > ## Ratio estimator using the heteroskedasticity parameter (Daniel Oehm)
  > # should match the ratio estmate above
  > dstratgh <- calibrate(dstrat,~enroll-1, total.enroll, variance=apistrat$enroll)
  > svytotal(~api.stu, dstratgh)
            total    SE
  api.stu 3190038 29566
  > 
  > ## individual boundary constraints as multiplicative values (Daniel Oehm)
  > bnds <- list(
  +   lower = c(1, 1, rep(-Inf, nrow(apistrat)-2)), 
  +   upper = c(1, 1, rep(Inf, nrow(apistrat)-2))) # the first two weights will remain unchanged the others are free to move
  > lapply(bnds, head)
  $lower
  [1]    1    1 -Inf -Inf -Inf -Inf
  
  $upper
  [1]   1   1 Inf Inf Inf Inf
  
  > dstratg1<-calibrate(dstrat, ~enroll-1, total.enroll, bounds = bnds, variance=apistrat$enroll)
  > svytotal(~api.stu, dstratg1)
            total    SE
  api.stu 3190133 29561
  > head(weights(dstrat))
      1     2     3     4     5     6 
  44.21 44.21 44.21 44.21 44.21 44.21 
  > head(weights(dstratg1))
         1        2        3        4        5        6 
  44.21000 44.21000 45.72055 45.72055 45.72055 45.72055 
  > all.equal(weights(dstrat)[1:2], weights(dstratg1)[1:2])
  [1] TRUE
  > 
  > ## individual boundary constraints as constant values (Daniel Oehm)
  > bnds <- list(
  +   lower = c(44.21, 44.21, rep(-Inf, nrow(apistrat)-2)), 
  +   upper = c(44.21, 44.21, rep(Inf, nrow(apistrat)-2))) # the first two weights will remain unchanged
  > lapply(bnds, head)
  $lower
  [1] 44.21 44.21  -Inf  -Inf  -Inf  -Inf
  
  $upper
  [1] 44.21 44.21   Inf   Inf   Inf   Inf
  
  > dstratg2<-calibrate(dstrat, ~enroll-1, total.enroll, bounds = bnds, bounds.const = TRUE, variance=apistrat$enroll)
  > svytotal(~api.stu, dstratg2)
            total    SE
  api.stu 3190133 29561
  > head(weights(dstrat))
      1     2     3     4     5     6 
  44.21 44.21 44.21 44.21 44.21 44.21 
  > head(weights(dstratg2))
         1        2        3        4        5        6 
  44.21000 44.21000 45.72055 45.72055 45.72055 45.72055 
  > all.equal(round(weights(dstrat)[1:2], 8), round(weights(dstratg2)[1:2]), 8) # minor rounding error but all good
  [1] TRUE
  > 
  > # sparse matrix support (Daniel Oehm)
  > dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > pop.totals<-c(`(Intercept)`=6194, stypeH=755, stypeM=1018)
  > dclus1g<-calibrate(dclus1, ~stype, pop.totals)
  > svymean(~api00, dclus1g)
          mean     SE
  api00 642.31 23.921
  > svytotal(~enroll, dclus1g)
           total     SE
  enroll 3680893 406293
  > 
  > pop.totals<-c(`(Intercept)`=6194, stypeH=755, stypeM=1018)
  > dclus1g<-calibrate(dclus1, ~stype, pop.totals, sparse = TRUE)
  > svymean(~api00, dclus1g)
          mean     SE
  api00 642.31 23.921
  > svytotal(~enroll, dclus1g)
           total     SE
  enroll 3680893 406293
  > 
  > pop.totals<-c(`(Intercept)`=6194, stypeH=755, stypeM=1018)
  > dclus1g<-calibrate(dclus1, ~stype, pop.totals, sparse = TRUE, calfun = "raking")
  > svymean(~api00, dclus1g)
          mean     SE
  api00 642.31 23.921
  > svytotal(~enroll, dclus1g)
           total     SE
  enroll 3680893 406293
  > 
  > proc.time()
     user  system elapsed 
     2.84    0.21    3.04 
Running the tests in 'tests/confintrep.R' failed.
Complete output:
  > 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)
  > dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
  > rclus2<-as.svrepdesign(dclus2)
  Warning message:
  In as.svrepdesign.default(dclus2) :
    Finite population corrections after first stage have been dropped
  > 
  > m<-svyglm(I(comp.imp=="Yes")~1, design=dclus2, family=quasibinomial)
  > if(anyNA(confint(m, method="likelihood"))) stop("NA in confint")
  > mrep<-svyglm(I(comp.imp=="Yes")~1, design=rclus2, family=quasibinomial)
  > if(anyNA(confint(mrep))) stop("NA in confint")
  > if(anyNA(confint(mrep, method="likelihood"))) stop("NA in confint")
  > 
  > proc.time()
     user  system elapsed 
     1.70    0.20    1.89 
Running the tests in 'tests/contrast-replicates.R' failed.
Complete output:
  > ## test use of replicates in svyby, svycontrast
  > 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)
  > 
  > meanlogs_without<-svyby(~log(enroll),~stype,svymean, design=rclus1,covmat=TRUE)
  > c_without<-svycontrast(meanlogs_without, quote(exp(E-H)))
  > vcov(c_without)
             contrast
  contrast 0.06862947
  > 
  > meanlogs_with<-svyby(~log(enroll),~stype,svymean, design=rclus1,covmat=TRUE,return.replicates=TRUE)
  > 
  > c_with<-svycontrast(meanlogs_with, quote(exp(E-H)))
  > 
  > v_with<- vcov(rclus1, c_with$replicates)
  > 
  > r<- attr(meanlogs_with, "replicates")
  > vr_with<-vcov(rclus1,exp(r[,1]-r[,2]))
  > 
  > stopifnot(all.equal(as.numeric(v_with),as.numeric(vr_with)))
  > stopifnot(all.equal(as.numeric(v_with),as.numeric(vcov(c_with))))
  > 
  > proc.time()
     user  system elapsed 
     1.54    0.20    1.71 
Running the tests in 'tests/coxph-termtest.R' failed.
Complete output:
  > 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
  
  > library(survival)
  > set.seed(2021-6-25)
  > test1 <- list(time=c(4,3,1,1,2,2,3),
  +               status=c(1,1,1,0,1,1,0),
  +               x1=as.factor(rbinom(7, 2, 0.5)),
  +               x=c(0,2,1,1,1,0,0))
  > # Fit a stratified model
  > mod_c <- coxph(Surv(time, status) ~ x1 + x, test1)
  > mod_d <- coxph(Surv(time, status) ~ x + x1, test1)
  > stopifnot(all.equal(regTermTest(mod_c, ~x1, df = Inf)[c("chisq","df","test.terms","p")],
  + regTermTest(mod_d, ~x1, df = Inf)[c("chisq","df","test.terms","p")]))
  > 
  > data(pbc, package="survival")
  > 
  > pbc$randomized<-with(pbc, !is.na(trt) & trt>0)
  > biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial)
  > pbc$randprob<-fitted(biasmodel)
  > if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0
  > 
  > dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
  > library(splines)
  > model<-svycoxph(formula = Surv(time, status > 0) ~ bili + protime + albumin+ns(bili,4)[,1:3], design = dpbc)
  > test<-regTermTest(model, ~ns(bili,4)[,1:3],method="LRT")
  > stopifnot(all.equal(test$chisq, 47.314, tolerance=1e-4))
  > stopifnot(all.equal(test$lambda, c(1.4764260, 1.0109836, 0.6923415),tolerance=1e-4))
  > 
  > proc.time()
     user  system elapsed 
     1.45    0.21    1.65 
Running the tests in 'tests/defftest.R' failed.
Complete output:
  > 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)
  > 
  > ## one-stage cluster sample
  > dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > 
  > # svyglm model
  > mod <- svyglm(api99 ~ enroll + api.stu, design = dclus1, deff = TRUE)
  > 
  > #deffs returned from svyglm model - implausibly high
  > deff(mod)
  (Intercept)      enroll     api.stu 
     1.919945    2.500983    2.686102 
  > #> (Intercept)      enroll     api.stu 
  > #>    351.3500    457.6799    491.5567
  > 
  > # run mod with same data and glm()
  > srs_mod <- glm(api99 ~ enroll + api.stu, data = apiclus1)
  > 
  > # manually calculate deffs
  > 
  > clust_se <- summary(mod)$coefficients[,2]
  > srs_se <- summary(srs_mod)$coefficients[,2]
  > 
  > deffs <- clust_se^2 / srs_se^2
  > stopifnot(all.equal(deffs, deff(mod)))
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.48    0.12    1.59 
Running the tests in 'tests/glm-scoping.R' failed.
Complete output:
  > ## bug report from Thomas Leeper, fixed in version 3.32-3
  > 
  > 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)
  > 
  > # pass `family` directly (WORKS!)
  > svyglm(api00~ell+meals+mobility, design=dstrat, family = gaussian())
  Stratified Independent Sampling design
  svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
      fpc = ~fpc)
  
  Call:  svyglm(formula = api00 ~ ell + meals + mobility, design = dstrat, 
      family = gaussian())
  
  Coefficients:
  (Intercept)          ell        meals     mobility  
     820.8873      -0.4806      -3.1415       0.2257  
  
  Degrees of Freedom: 199 Total (i.e. Null);  194 Residual
  Null Deviance:	    3023000 
  Residual Deviance: 1029000 	AIC: 2308
  > 
  > # passing `family` via ... (WORKS!)
  > myfun1 <- function(formula, design, ...) {
  +     svyglm(formula, design = design, ...)
  + }
  > myfun1(api00~ell+meals+mobility, design=dstrat, family = gaussian())
  Stratified Independent Sampling design
  svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
      fpc = ~fpc)
  
  Call:  svyglm(formula = formula, design = design, family = ..1)
  
  Coefficients:
  (Intercept)          ell        meals     mobility  
     820.8873      -0.4806      -3.1415       0.2257  
  
  Degrees of Freedom: 199 Total (i.e. Null);  194 Residual
  Null Deviance:	    3023000 
  Residual Deviance: 1029000 	AIC: 2308
  > 
  > # passing `family` via default argument (DOES NOT WORK!)
  > myfun2 <- function(formula, design, family = gaussian()) {
  +     svyglm(formula, design = design, family = family)
  + }
  > myfun2(api00~ell+meals+mobility, design=dstrat, family = gaussian())
  Stratified Independent Sampling design
  svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
      fpc = ~fpc)
  
  Call:  svyglm(formula = formula, design = design, family = family)
  
  Coefficients:
  (Intercept)          ell        meals     mobility  
     820.8873      -0.4806      -3.1415       0.2257  
  
  Degrees of Freedom: 199 Total (i.e. Null);  194 Residual
  Null Deviance:	    3023000 
  Residual Deviance: 1029000 	AIC: 2308
  > 
  > proc.time()
     user  system elapsed 
     1.45    0.23    1.67 
Running the tests in 'tests/kalton.R' failed.
Complete output:
  > 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
  
  > 
  > ab<-expand.grid(a=factor(1:4),b=factor(1:3))
  > 
  > kaltonsample<-ab[rep(1:12,c(20,50,100,30,40,140,50,100,40,310,50,70)),]
  > 
  > kaltonpop<-ab[rep(1:12,c(80,60,170,55,40,150,60,165,55,340,200,125)),]
  > 
  > jointpop<-colSums(model.matrix(~a*b,kaltonpop))
  > marginalpop<-colSums(model.matrix(~a+b,kaltonpop))
  > gregpop<-colSums(model.matrix(~as.numeric(a)+as.numeric(b),kaltonpop))
  > 
  > dkalton<-svydesign(id=~1,data=kaltonsample)
  Warning message:
  In svydesign.default(id = ~1, data = kaltonsample) :
    No weights or probabilities supplied, assuming equal probability
  > 
  > dps<-postStratify(dkalton,~a+b,xtabs(~a+b,kaltonpop))
  > 
  > drake<-rake(dkalton, list(~a,~b),list(xtabs(~a,kaltonpop),xtabs(~b,kaltonpop)),control=list(epsilon=0.0001))
  > 
  > dcalps<-calibrate(dkalton, ~a*b, jointpop)
  > dcalrake<-calibrate(dkalton,~a+b, marginalpop, calfun="raking")
  > dlinear<-calibrate(dkalton, ~a+b, marginalpop)
  > 
  > dtrunclinear<-calibrate(dkalton, ~a+b, marginalpop,bounds=c(0.5,2.2))
  > 
  > dlogit<-calibrate(dkalton, ~a+b, marginalpop,bounds=c(0.5,2.2),calfun="logit")
  > 
  > dgreg<-calibrate(dkalton,~as.numeric(a)+as.numeric(b), gregpop)
  > 
  > 
  > #table A
  >  round(svytable(~a+b,dps)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 4.00 1.00 1.38
    2 1.20 1.07 1.10
    3 1.70 1.20 4.00
    4 1.83 1.65 1.79
  >  round(svytable(~a+b,dcalps)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 4.00 1.00 1.37
    2 1.20 1.07 1.10
    3 1.70 1.20 4.00
    4 1.83 1.65 1.79
  > 
  > #table B
  >  round(svytable(~a+b,drake)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.81 1.45 2.02
    2 1.08 0.87 1.21
    3 2.20 1.76 2.45
    4 1.83 1.47 2.04
  >  round(svytable(~a+b,dcalrake)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.81 1.45 2.02
    2 1.08 0.87 1.21
    3 2.20 1.76 2.45
    4 1.83 1.47 2.04
  > 
  > #table C
  > round(svytable(~a+b,dlinear)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.82 1.50 1.97
    2 1.09 0.78 1.24
    3 2.19 1.88 2.34
    4 1.83 1.52 1.98
  > 
  > #table D
  > round(svytable(~a+b,dgreg)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.21 1.17 1.14
    2 1.43 1.40 1.36
    3 1.66 1.62 1.59
    4 1.88 1.85 1.81
  > 
  > #table G
  > round(svytable(~a+b,dlogit)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.87 1.46 1.98
    2 1.08 0.74 1.27
    3 2.17 2.09 2.18
    4 1.89 1.49 1.99
  > 
  > #table G
  > round(svytable(~a+b,dtrunclinear)/xtabs(~a+b,kaltonsample),2)
     b
  a      1    2    3
    1 1.81 1.48 1.99
    2 1.08 0.75 1.26
    3 2.20 2.00 2.20
    4 1.83 1.50 2.00
  > 
  > proc.time()
     user  system elapsed 
     1.65    0.14    1.78 
Running the tests in 'tests/lonely.psu.R' failed.
Complete output:
  > 
  > ## lonely PSUs by design
  > 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)
  > ## not certainty PSUs by fpc
  > ds<-svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = apiclus1)
  > summary(ds)
  Stratified Independent Sampling design (with replacement)
  svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = apiclus1)
  Probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
  Stratum Sizes: 
             61 135 178 197 255 406 413 437 448 510 568 637 716 778 815
  obs        13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  design.PSU 13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  actual.PSU 13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  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"      
  > 
  > options(survey.lonely.psu="fail")
  > try(svymean(~api00,ds))
  Error in onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  : 
    Stratum (413) has only one PSU at stage 1
  > try(svymean(~api00, as.svrepdesign(ds)))
  Error in jknweights(design$strata[, 1], design$cluster[, 1], fpc = fpc,  : 
    Stratum413has only one PSU
  > options(survey.lonely.psu="remove")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.8058
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.8058
  > options(survey.lonely.psu="certainty")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.8058
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.8058
  > options(survey.lonely.psu="adjust")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.8281
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.8267
  > options(survey.lonely.psu="average")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 6.0096
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.8217
  > 
  > ## fpc specified
  > fpc<-ifelse(apiclus1$dnum==413, 1,1000)
  > ds<-svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = apiclus1,fpc=fpc)
  > summary(ds)
  Stratified Independent Sampling design
  svydesign(id = ~1, weights = ~pw, strata = ~dnum, 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 
  Stratum Sizes: 
             61 135 178 197 255 406 413 437 448 510 568 637 716 778 815
  obs        13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  design.PSU 13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  actual.PSU 13  34   4  13  16   2   1   4  12  21   9  11  37   2   4
  Population stratum sizes (PSUs): 
   135  178  197  255  406  413  437  448  510  568   61  637  716  778  815 
  1000 1000 1000 1000 1000    1 1000 1000 1000 1000 1000 1000 1000 1000 1000 
  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"      
  > 
  > options(survey.lonely.psu="fail")
  > try(svymean(~api00,ds))
          mean     SE
  api00 644.17 5.7344
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.7344
  > options(survey.lonely.psu="remove")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.7344
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.7344
  > options(survey.lonely.psu="certainty")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.7344
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.7344
  > options(survey.lonely.psu="adjust")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.7344
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.7344
  > options(survey.lonely.psu="average")
  > svymean(~api00,ds)
          mean     SE
  api00 644.17 5.7344
  > svymean(~api00, as.svrepdesign(ds))
          mean     SE
  api00 644.17 5.7501
  > 
  > rs<-as.svrepdesign(ds)
  > svytotal(~api00,rs)
          total    SE
  api00 3989986 35616
  > SE(svytotal(~api00,subset(rs, dnum==413)))==0
  [1] TRUE
  > 
  > ## lonely PSUs after subsetting
  > ds<-svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = subset(apiclus1,dnum !=413))
  > ds1<-ds[-31,]
  > summary(ds1)
  Stratified Independent Sampling design (with replacement)
  svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = subset(apiclus1, 
      dnum != 413))
  Probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
  Stratum Sizes: 
             61 135 178 197 255 406 437 448 510 568 637 716 778 815
  obs        13  34   4  13  16   1   4  12  21   9  11  37   2   4
  design.PSU 13  34   4  13  16   2   4  12  21   9  11  37   2   4
  actual.PSU 13  34   4  13  16   1   4  12  21   9  11  37   2   4
  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"      
  > 
  > options(survey.lonely.psu="fail")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  > options(survey.lonely.psu="remove")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  > options(survey.lonely.psu="certainty")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  > options(survey.lonely.psu="adjust")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  > options(survey.lonely.psu="average")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  > 
  > ## with adjustment
  > options(survey.adjust.domain.lonely=TRUE)
  > ds<-svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = subset(apiclus1,dnum !=413))
  > ds1<-ds[-31,]
  > summary(ds1)
  Stratified Independent Sampling design (with replacement)
  svydesign(id = ~1, weights = ~pw, strata = ~dnum, data = subset(apiclus1, 
      dnum != 413))
  Probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.02954 0.02954 0.02954 0.02954 0.02954 0.02954 
  Stratum Sizes: 
             61 135 178 197 255 406 437 448 510 568 637 716 778 815
  obs        13  34   4  13  16   1   4  12  21   9  11  37   2   4
  design.PSU 13  34   4  13  16   2   4  12  21   9  11  37   2   4
  actual.PSU 13  34   4  13  16   1   4  12  21   9  11  37   2   4
  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"      
  > 
  > options(survey.lonely.psu="fail")
  > try(svymean(~api00,ds1))
          mean     SE
  api00 645.14 5.8909
  Warning message:
  In onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  :
    Stratum (406) has only one PSU at stage 1
  > options(survey.lonely.psu="remove")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  Warning message:
  In onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  :
    Stratum (406) has only one PSU at stage 1
  > options(survey.lonely.psu="certainty")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.8909
  Warning message:
  In onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  :
    Stratum (406) has only one PSU at stage 1
  > options(survey.lonely.psu="adjust")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 5.9119
  Warning message:
  In onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  :
    Stratum (406) has only one PSU at stage 1
  > options(survey.lonely.psu="average")
  > svymean(~api00,ds1)
          mean     SE
  api00 645.14 6.0914
  Warning message:
  In onestrat(`attr<-`(x[index, , drop = FALSE], "recentering", recentering),  :
    Stratum (406) has only one PSU at stage 1
  > 
  > ## checks for `svytotal()`
  > 
  > df_w_singleton <- data.frame(
  +   Stratum = c(1, 1, 2, 2, 3),
  +   PSU = c(1, 2, 3, 4, 5),
  +   Design_Weight = c(10.5, 10.5, 20, 20, 15),
  +   Sex = c("M", "F", "F", "M",  "M"),
  +   Age = c(30.5, 40.5, 35, 52, 44),
  +   Height = c(6.2, 5.0, 5.3, 5.7, 5.5)
  + )
  > 
  > design_w_singleton <- survey::svydesign(
  +   data = df_w_singleton,
  +   ids = ~ PSU, strata = ~ Stratum,
  +   weights = ~ Design_Weight
  + )
  > 
  > options("survey.lonely.psu" = "remove")
  > 
  > stopifnot(all.equal(
  +   target = 126625,
  +   current = as.numeric(
  +     vcov(svytotal(x = ~ Age, design = design_w_singleton))
  +   )
  + ))
  > 
  > options("survey.lonely.psu" = "certainty")
  > 
  > stopifnot(all.equal(
  +   target = 126625,
  +   current = as.numeric(
  +     vcov(svytotal(x = ~ Age, design = design_w_singleton))
  +   )
  + ))
  > 
  > options("survey.lonely.psu" = "adjust")
  > 
  > stopifnot(all.equal(
  +   target = 127579.8,
  +   current = as.numeric(
  +     vcov(svytotal(x = ~ Age, design = design_w_singleton))
  +   ),
  +   scale = 127579.8, tolerance = 0.000001
  + ))
  > 
  > proc.time()
     user  system elapsed 
     1.92    0.26    2.18 
Running the tests in 'tests/multistage-rcpp.R' failed.
Complete output:
  > ## temporarily disable while it's being fixed
  > ## q("no")
  > 
  > ## 
  > ## Check that multistage() and multistage_rcpp() give same results
  > ## for different options.
  > ##
  > 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
  
  > 
  > # Check for a simple random sample ----
  > 
  >   data('api', package = 'survey')
  >   api_srs_design <- svydesign(
  +     data = apisrs,
  +     ids = ~ 1,
  +     weights = ~ 1
  +   )
  >   
  >   x <- as.matrix(api_srs_design$variables[,c('api00','api99')] / api_srs_design$prob)
  >   
  >   base_r_result <- survey:::multistage(x = x,
  +                                        clusters = api_srs_design$cluster,
  +                                        stratas = api_srs_design$strata,
  +                                        nPSUs = api_srs_design$fpc$sampsize, fpcs = api_srs_design$fpc$popsize,
  +                                        lonely.psu=getOption("survey.lonely.psu"),
  +                                        one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   rcpp_result <- survey:::multistage_rcpp(x = x,
  +                                           clusters = api_srs_design$cluster,
  +                                           stratas = api_srs_design$strata,
  +                                           nPSUs = api_srs_design$fpc$sampsize, fpcs = api_srs_design$fpc$popsize,
  +                                           lonely.psu=getOption("survey.lonely.psu"),
  +                                           one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   if (!isTRUE(all.equal(base_r_result, rcpp_result))) {
  +     stop("Differences between `multistage()` and `multistage_rcpp()` for SRS")
  +   }
  >   
  > # Check for a stratified simple random sample ----
  >   
  >   apistrat_design <- svydesign(
  +     data = apistrat,
  +     id =~ 1,
  +     strata =~ stype,
  +     weights =~ pw,
  +     fpc =~ fpc
  +   )
  >   
  >   x <- as.matrix(apistrat_design$variables[,c('api00','api99')] / apistrat_design$prob)
  >   
  >   base_r_result <- survey:::multistage(x = x,
  +                                        clusters = apistrat_design$cluster,
  +                                        stratas = apistrat_design$strata,
  +                                        nPSUs = apistrat_design$fpc$sampsize, fpcs = apistrat_design$fpc$popsize,
  +                                        lonely.psu=getOption("survey.lonely.psu"),
  +                                        one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   rcpp_result <- survey:::multistage_rcpp(x = x,
  +                                           clusters = apistrat_design$cluster,
  +                                           stratas = apistrat_design$strata,
  +                                           nPSUs = apistrat_design$fpc$sampsize, fpcs = apistrat_design$fpc$popsize,
  +                                           lonely.psu=getOption("survey.lonely.psu"),
  +                                           one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   if (!isTRUE(all.equal(base_r_result, rcpp_result))) {
  +     stop("Differences between `multistage()` and `multistage_rcpp()` for stratified sample")
  +   }
  >   
  >   ##_ Check whether expected true zeroes are actually computed as zeroes ----
  >   
  >   x <- as.matrix(model.matrix(~ -1 + stype, data = apistrat_design$variables) / (sum(apistrat_design$prob)/apistrat_design$prob))
  >   
  >   base_r_result <- survey:::multistage(x = x,
  +                                        clusters = apistrat_design$cluster,
  +                                        stratas = apistrat_design$strata,
  +                                        nPSUs = apistrat_design$fpc$sampsize, fpcs = apistrat_design$fpc$popsize,
  +                                        lonely.psu=getOption("survey.lonely.psu"),
  +                                        one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   rcpp_result <- survey:::multistage_rcpp(x = x,
  +                                           clusters = apistrat_design$cluster,
  +                                           stratas = apistrat_design$strata,
  +                                           nPSUs = apistrat_design$fpc$sampsize, fpcs = apistrat_design$fpc$popsize,
  +                                           lonely.psu=getOption("survey.lonely.psu"),
  +                                           one.stage=TRUE, stage = 1, cal = NULL)
  >   
  >   if (!all(rcpp_result == 0) | !all(diag(rcpp_result) >= 0)) {
  +     stop("Computations which should equal zero do not.")
  +   }
  > 
  > # Create a stratified, multistage design ----
  >   
  >   data(mu284)
  >   
  >   ## Create three strata, the third of which has only one PSU
  >   mu284_stratified <- rbind(
  +     transform(mu284, stratum = 1,
  +               y2 = y1 + c(0, 2, 1, 0, -1, 2, 0, 3, 0, 0, -1, -1, 1, -1, 0)),
  +     transform(mu284, stratum = 2,
  +               y2 = y1 + c(-1, 0, -1, 0, -1, 0, 0, 1, 2, -1, 1, 1, 0, 1, 0)),
  +     transform(mu284[1,], stratum = 3,
  +               y2 = y1 + 2)
  +   )
  >   
  >   ## Create domain variables which yield a lonely PSU or a lonely SSU
  >   mu284_stratified[['DOMAIN_W_LONELY_PSU']] <- mu284_stratified$stratum == 1 | (mu284_stratified$stratum == 2 & mu284_stratified$id1 == 19)
  >   mu284_stratified[['DOMAIN_W_LONELY_SSU']] <- mu284_stratified$stratum == 1 | (mu284_stratified$stratum == 2 & mu284_stratified$id1 == 19 & mu284_stratified$id2 == 1)
  >   
  >   ## Create a survey design object with no lonely PSUs
  >   dmu284_strat <- svydesign(data = subset(mu284_stratified, stratum != 3),
  +                             strata = ~ stratum, nest = TRUE,
  +                             id = ~ id1 + id2, fpc = ~ n1 + n2)
  >   ## Create a survey design with one lonely PSU
  >   dmu284_strat_w_lonely <- svydesign(data = mu284_stratified,
  +                                      strata = ~ stratum, nest = TRUE,
  +                                      id = ~ id1 + id2, fpc = ~ n1 + n2)
  >   ## Create subsetted design objects, subsetted to domain with a lonely PSU or SSU
  >   dmu284_strat_domain_lonely_psu <- subset(dmu284_strat, DOMAIN_W_LONELY_PSU)
  >   dmu284_strat_domain_lonely_ssu <- subset(dmu284_strat, DOMAIN_W_LONELY_SSU)
  >   
  >   
  > # Check same results for different values of 'lonely PSU' options ----
  >   
  >   lonely_psu_options <- c(certainty = 'certainty',
  +                           remove    = 'remove',
  +                           average   = 'average',
  +                           adjust    = 'adjust')
  >   one_stage_options <- c(TRUE, FALSE)
  >   
  >   design_lonely_psu_comparisons <- expand.grid('survey.lonely.psu' = lonely_psu_options,
  +                                                'one.stage' = one_stage_options,
  +                                                stringsAsFactors = FALSE)
  >   design_lonely_psu_comparisons[['results_match']] <- NA
  >   
  >   ##_ Check for lonely PSUs caused by design ----
  >   for (i in seq_len(nrow(design_lonely_psu_comparisons))) {
  +     survey.lonely.psu <- design_lonely_psu_comparisons[['survey.lonely.psu']][i]
  +     one.stage <- design_lonely_psu_comparisons[['one.stage']][i]
  +     options('survey.lonely.psu' = survey.lonely.psu)
  +     
  +     x <- as.matrix(dmu284_strat_w_lonely$variables[,c('y1','y2')] / dmu284_strat_w_lonely$prob)
  +     
  +     base_r_result <- survey:::multistage(x = x,
  +                                          clusters = dmu284_strat_w_lonely$cluster,
  +                                          stratas = dmu284_strat_w_lonely$strata,
  +                                          nPSUs = dmu284_strat_w_lonely$fpc$sampsize, fpcs = dmu284_strat_w_lonely$fpc$popsize,
  +                                          lonely.psu=getOption("survey.lonely.psu"),
  +                                          one.stage=one.stage, stage = 1, cal = NULL)
  +     
  +     rcpp_result <- survey:::multistage_rcpp(x = x,
  +                                             clusters = dmu284_strat_w_lonely$cluster,
  +                                             stratas = dmu284_strat_w_lonely$strata,
  +                                             nPSUs = dmu284_strat_w_lonely$fpc$sampsize, fpcs = dmu284_strat_w_lonely$fpc$popsize,
  +                                             lonely.psu=getOption("survey.lonely.psu"),
  +                                             one.stage=one.stage, stage = 1, cal = NULL)
  +     
  +       design_lonely_psu_comparisons[['results_match']][i] <- isTRUE(all.equal(base_r_result, rcpp_result))
  +   }
  >   
  >   print(design_lonely_psu_comparisons)
    survey.lonely.psu one.stage results_match
  1         certainty      TRUE          TRUE
  2            remove      TRUE          TRUE
  3           average      TRUE          TRUE
  4            adjust      TRUE          TRUE
  5         certainty     FALSE          TRUE
  6            remove     FALSE          TRUE
  7           average     FALSE          TRUE
  8            adjust     FALSE          TRUE
  >   if (!all(design_lonely_psu_comparisons$results_match)) {
  +     stop("Results for design lonely PSUs differ between base R and Rcpp implementations.")
  +   }
  > 
  >   ##_ Check for domain lonely PSUs caused by subsetting ----
  >   
  >   domain_lonely_psu_comparisons <- expand.grid('survey.lonely.psu' = lonely_psu_options,
  +                                                'survey.adjust.domain.lonely' = c(FALSE, TRUE),
  +                                                'one.stage' = one_stage_options,
  +                                                stringsAsFactors = FALSE)
  >   domain_lonely_psu_comparisons[['results_match']] <- NA
  >   
  >   for (i in seq_len(nrow(domain_lonely_psu_comparisons))) {
  +       
  +     options('survey.lonely.psu' = domain_lonely_psu_comparisons[['survey.lonely.psu']][i])
  +     options('survey.adjust.domain.lonely' = domain_lonely_psu_comparisons[['survey.adjust.domain.lonely']][i])
  +     one.stage <- domain_lonely_psu_comparisons[['one.stage']][i]
  +     
  +     
  +     x <- as.matrix(dmu284_strat_domain_lonely_psu$variables[,c('y1','y2')] / dmu284_strat_domain_lonely_psu$prob)
  +     suppressWarnings({
  +       base_r_result <- survey:::multistage(x = x,
  +                                            clusters = dmu284_strat_domain_lonely_psu$cluster,
  +                                            stratas = dmu284_strat_domain_lonely_psu$strata,
  +                                            nPSUs = dmu284_strat_domain_lonely_psu$fpc$sampsize,
  +                                            fpcs = dmu284_strat_domain_lonely_psu$fpc$popsize,
  +                                            lonely.psu=getOption("survey.lonely.psu"),
  +                                            one.stage=one.stage, stage = 1, cal = NULL)
  +     })
  +     
  +     rcpp_result <- survey:::multistage_rcpp(x = x,
  +                                             clusters = dmu284_strat_domain_lonely_psu$cluster,
  +                                             stratas = dmu284_strat_domain_lonely_psu$strata,
  +                                             nPSUs = dmu284_strat_domain_lonely_psu$fpc$sampsize,
  +                                             fpcs = dmu284_strat_domain_lonely_psu$fpc$popsize,
  +                                             lonely.psu=getOption("survey.lonely.psu"),
  +                                             one.stage=one.stage, stage = 1, cal = NULL)
  +     
  +     domain_lonely_psu_comparisons[['results_match']][i] <- isTRUE(all.equal(base_r_result, rcpp_result))
  +   }
  Warning messages:
  1: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  2: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  3: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  4: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  5: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  6: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  7: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  8: In arma_multistage(Y = as.matrix(x), samp_unit_ids = clusters, strata_ids = stratas,  :
    At least one stratum contains only one PSU at stage 1
  >   
  >   print(domain_lonely_psu_comparisons)
     survey.lonely.psu survey.adjust.domain.lonely one.stage results_match
  1          certainty                       FALSE      TRUE          TRUE
  2             remove                       FALSE      TRUE          TRUE
  3            average                       FALSE      TRUE          TRUE
  4             adjust                       FALSE      TRUE          TRUE
  5          certainty                        TRUE      TRUE          TRUE
  6             remove                        TRUE      TRUE          TRUE
  7            average                        TRUE      TRUE          TRUE
  8             adjust                        TRUE      TRUE          TRUE
  9          certainty                       FALSE     FALSE          TRUE
  10            remove                       FALSE     FALSE          TRUE
  11           average                       FALSE     FALSE          TRUE
  12            adjust                       FALSE     FALSE          TRUE
  13         certainty                        TRUE     FALSE          TRUE
  14            remove                        TRUE     FALSE          TRUE
  15           average                        TRUE     FALSE          TRUE
  16            adjust                        TRUE     FALSE          TRUE
  >   if (!all(domain_lonely_psu_comparisons$results_match)) {
  +     stop("Results for domain lonely PSUs differ between base R and Rcpp implementations.")
  +   }
  >   
  > # Checks for edge cases ----
  >   
  >   ## Zero records ----
  >   
  >   zero_row_design <- dmu284_strat |> subset(id1 == "Nonexistent")
  >   
  >   zero_row_result <- survey:::multistage_rcpp(
  +     x = as.matrix(zero_row_design$variables[,c('y1', 'y2')]),
  +     clusters = zero_row_design$cluster,
  +     stratas = zero_row_design$strata,
  +     nPSUs = zero_row_design$fpc$sampsize, 
  +     fpcs = zero_row_design$fpc$popsize,
  +     lonely.psu=getOption("survey.lonely.psu"),
  +     one.stage=one.stage, stage = 1, cal = NULL
  +   )
  >   
  >   if (any(is.na(zero_row_result)) || any(zero_row_result != 0)) {
  +     stop("`multistage_rcpp()` should return zeros for zero-row inputs.")
  +   }
  >   
  > 
  > proc.time()
     user  system elapsed 
     1.89    0.14    2.01 
Running the tests in 'tests/na_action.R' failed.
Complete output:
  > ## from Terry Therneau
  > 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
  
  > load("naa.rda")
  > 
  > fit1e <- svyglm( pseudo ~ age34 + ccr5 + factor(times),  design= adata.s,na.action=na.exclude)
  > fit1o <- svyglm( pseudo ~ age34 + ccr5 + factor(times),  design= adata.s)
  > all.equal(coef(fit1e),coef(fit1o))
  [1] TRUE
  > all.equal(vcov(fit1e),vcov(fit1o))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.46    0.20    1.65 
Running the tests in 'tests/nwts.R' failed.
Complete output:
  > 
  > ## examples from Breslow & Chatterjee: Applied Statistics 1999 No. 4, p458
  > ## data from Norman Breslow's web page.
  > 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
  
  > load("nwts.rda")
  > nwtsnb<-nwts
  > nwtsnb$case<-nwts$case-nwtsb$case
  > nwtsnb$control<-nwts$control-nwtsb$control
  > 
  > a<-rbind(nwtsb,nwtsnb)
  > a$in.ccs<-rep(c(TRUE,FALSE),each=16)
  > 
  > b<-rbind(a,a)
  > b$rel<-rep(c(1,0),each=32)
  > b$n<-ifelse(b$rel,b$case,b$control)
  > 
  > index<-rep(1:64,b$n)
  > 
  > nwt.exp<-b[index,c(1:3,6,7)]
  > nwt.exp$id<-1:4088
  > 
  > dccs2<-twophase(id=list(~id,~id),subset=~in.ccs,
  +                 strata=list(NULL,~interaction(instit,rel)),data=nwt.exp)
  > 
  > dccs8<-twophase(id=list(~id,~id),subset=~in.ccs,
  +                 strata=list(NULL,~interaction(instit,stage,rel)),data=nwt.exp)
  > 
  > gccs8<-calibrate(dccs2,phase=2,formula=~interaction(instit,stage,rel))
  > 
  > summary(svyglm(rel~factor(stage)*factor(histol),family=quasibinomial,design=dccs2))
  
  Call:
  svyglm(formula = rel ~ factor(stage) * factor(histol), design = dccs2, 
      family = quasibinomial)
  
  Survey design:
  twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
      subset = subset, data = data, pps = pps)
  
  Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)                     -2.5701     0.1288 -19.955  < 2e-16 ***
  factor(stage)2                   0.5482     0.1979   2.769 0.005708 ** 
  factor(stage)3                   0.4791     0.2032   2.359 0.018515 *  
  factor(stage)4                   1.0037     0.2592   3.872 0.000114 ***
  factor(histol)2                  1.3505     0.3107   4.346 1.51e-05 ***
  factor(stage)2:factor(histol)2   0.1152     0.4410   0.261 0.793876    
  factor(stage)3:factor(histol)2   0.5066     0.4241   1.194 0.232548    
  factor(stage)4:factor(histol)2   0.9785     0.6214   1.575 0.115615    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.000876)
  
  Number of Fisher Scoring iterations: 4
  
  > summary(svyglm(rel~factor(stage)*factor(histol),family=quasibinomial,design=dccs8))
  
  Call:
  svyglm(formula = rel ~ factor(stage) * factor(histol), design = dccs8, 
      family = quasibinomial)
  
  Survey design:
  twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
      subset = subset, data = data, pps = pps)
  
  Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)                    -2.71604    0.10827 -25.085  < 2e-16 ***
  factor(stage)2                  0.78141    0.14726   5.306 1.35e-07 ***
  factor(stage)3                  0.80093    0.15250   5.252 1.80e-07 ***
  factor(stage)4                  1.07293    0.17817   6.022 2.33e-09 ***
  factor(histol)2                 1.45836    0.31780   4.589 4.96e-06 ***
  factor(stage)2:factor(histol)2 -0.04743    0.43495  -0.109    0.913    
  factor(stage)3:factor(histol)2  0.28064    0.41298   0.680    0.497    
  factor(stage)4:factor(histol)2  0.90983    0.63774   1.427    0.154    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.000876)
  
  Number of Fisher Scoring iterations: 4
  
  > summary(svyglm(rel~factor(stage)*factor(histol),family=quasibinomial,design=gccs8))
  
  Call:
  svyglm(formula = rel ~ factor(stage) * factor(histol), design = gccs8, 
      family = quasibinomial)
  
  Survey design:
  calibrate(dccs2, phase = 2, formula = ~interaction(instit, stage, 
      rel))
  
  Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)                    -2.71604    0.10878 -24.968  < 2e-16 ***
  factor(stage)2                  0.78141    0.14729   5.305 1.35e-07 ***
  factor(stage)3                  0.80093    0.15212   5.265 1.68e-07 ***
  factor(stage)4                  1.07293    0.17905   5.993 2.77e-09 ***
  factor(histol)2                 1.45836    0.31757   4.592 4.88e-06 ***
  factor(stage)2:factor(histol)2 -0.04743    0.43432  -0.109    0.913    
  factor(stage)3:factor(histol)2  0.28064    0.41231   0.681    0.496    
  factor(stage)4:factor(histol)2  0.90983    0.63187   1.440    0.150    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.000876)
  
  Number of Fisher Scoring iterations: 4
  
  > 
  > ## check subsets of calibrated designs.
  > summary(svyglm(rel~factor(stage),
  +                family=quasibinomial,design=subset(dccs8,histol==1)))
  
  Call:
  svyglm(formula = rel ~ factor(stage), design = subset(dccs8, 
      histol == 1), family = quasibinomial)
  
  Survey design:
  subset(dccs8, histol == 1)
  
  Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -2.7160     0.1083 -25.085  < 2e-16 ***
  factor(stage)2   0.7814     0.1473   5.306 1.48e-07 ***
  factor(stage)3   0.8009     0.1525   5.252 1.97e-07 ***
  factor(stage)4   1.0729     0.1782   6.022 2.73e-09 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.001333)
  
  Number of Fisher Scoring iterations: 4
  
  Warning messages:
  1: In `[.twophase2`(x, r, ) : 1 strata have only one PSU in this subset.
  2: In summary.glm(g) :
    observations with zero weight not used for calculating dispersion
  3: In summary.glm(glm.object) :
    observations with zero weight not used for calculating dispersion
  4: In `[.twophase2`(design, nas == 0, ) :
    1 strata have only one PSU in this subset.
  5: In `[.twophase2`(design, nas == 0, ) :
    1 strata have only one PSU in this subset.
  > summary(svyglm(rel~factor(stage),
  +                family=quasibinomial,design=subset(gccs8,histol==1)))
  
  Call:
  svyglm(formula = rel ~ factor(stage), design = subset(gccs8, 
      histol == 1), family = quasibinomial)
  
  Survey design:
  subset(gccs8, histol == 1)
  
  Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -2.7160     0.1082 -25.105  < 2e-16 ***
  factor(stage)2   0.7814     0.1457   5.363 1.10e-07 ***
  factor(stage)3   0.8009     0.1504   5.324 1.34e-07 ***
  factor(stage)4   1.0729     0.1759   6.101 1.70e-09 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for quasibinomial family taken to be 1.001333)
  
  Number of Fisher Scoring iterations: 4
  
  Warning messages:
  1: In summary.glm(g) :
    observations with zero weight not used for calculating dispersion
  2: In summary.glm(glm.object) :
    observations with zero weight not used for calculating dispersion
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.79    0.54    3.32 
Running the tests in 'tests/quantile.R' failed.
Complete output:
  > 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
  
  > set.seed(42)
  > 
  > df<-data.frame(x=exp(rnorm(1000)))
  > df$y<-round(df$x,1)
  > ddf<-svydesign(id=~1,data=df)
  Warning message:
  In svydesign.default(id = ~1, data = df) :
    No weights or probabilities supplied, assuming equal probability
  > rdf<-as.svrepdesign(ddf)
  > 
  > SE(oldsvyquantile(~x,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE))
        0.01        0.1        0.5        0.9       0.99 
  0.01543792 0.01264453 0.03386225 0.16131948 2.10605799 
  > 
  > SE(oldsvyquantile(~x,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE))
                 x
  q0.01 0.01554946
  q0.1  0.01425828
  q0.5  0.03428012
  q0.9  0.16457488
  q0.99 1.88672412
  > 
  > SE(oldsvyquantile(~x,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,df=Inf))
        0.01        0.1        0.5        0.9       0.99 
  0.01545209 0.01265608 0.03388011 0.16145776 2.10061576 
  > 
  > SE(oldsvyquantile(~x,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,df=Inf))
                 x
  q0.01 0.01556240
  q0.1  0.01420086
  q0.5  0.03430175
  q0.9  0.16461019
  q0.99 1.88818131
  > 
  > 
  > oldsvyquantile(~y,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,ties="rounded",interval.type="betaWald")
  $quantiles
          0.01       0.1       0.5  0.9 0.99
  y 0.02352941 0.2230769 0.9340909 3.55  9.4
  
  $CIs
  , , y
  
               0.01       0.1       0.5      0.9      0.99
  (lower 0.01594200 0.2020115 0.8578143 3.207543  7.568142
  upper) 0.03287947 0.2467045 1.0059815 3.862241 14.978632
  
  
  > 
  > oldsvyquantile(~y,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE)
  Statistic:
          y
  q0.01 0.1
  q0.1  0.3
  q0.5  1.0
  q0.9  3.6
  q0.99 9.4
  SE:
                 y
  q0.01 0.02547977
  q0.1  0.02547977
  q0.5  0.02547977
  q0.9  0.15287859
  q0.99 1.88302540
  > 
  > oldsvyquantile(~y,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,ties="rounded",interval.type="betaWald",df=Inf)
  $quantiles
          0.01       0.1       0.5  0.9 0.99
  y 0.02352941 0.2230769 0.9340909 3.55  9.4
  
  $CIs
  , , y
  
               0.01       0.1       0.5      0.9      0.99
  (lower 0.01594200 0.2020115 0.8578143 3.207543  7.568142
  upper) 0.03287947 0.2467045 1.0059815 3.862241 14.978632
  
  
  > 
  > oldsvyquantile(~y,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE, df=Inf)
  Statistic:
          y
  q0.01 0.1
  q0.1  0.3
  q0.5  1.0
  q0.9  3.6
  q0.99 9.4
  SE:
                 y
  q0.01 0.02551067
  q0.1  0.02551067
  q0.5  0.02551067
  q0.9  0.15306404
  q0.99 1.88450896
  > 
  > 
  > 
  > df<-data.frame(x=exp(rnorm(20)))
  > df$y<-round(df$x,1)
  > ddf<-svydesign(id=~1,data=df)
  Warning message:
  In svydesign.default(id = ~1, data = df) :
    No weights or probabilities supplied, assuming equal probability
  > rdf<-as.svrepdesign(ddf)
  > SE(oldsvyquantile(~x,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE))
        0.01        0.1        0.5        0.9       0.99 
  0.03576771 0.07329496 0.30149917 2.03440784 1.88468597 
  > 
  > SE(oldsvyquantile(~x,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE))
                 x
  q0.01 0.00000000
  q0.1  0.04598541
  q0.5  0.25943731
  q0.9  2.12343073
  q0.99 1.91998924
  > 
  > SE(oldsvyquantile(~x,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,df=Inf))
        0.01        0.1        0.5        0.9       0.99 
  0.03637398 0.07470310 0.29938115 2.16381289 2.00708994 
  > 
  > SE(oldsvyquantile(~x,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,df=Inf))
                 x
  q0.01 0.00000000
  q0.1  0.04732818
  q0.5  0.24957966
  q0.9  2.25211815
  q0.99 2.04043154
  > 
  > 
  > oldsvyquantile(~y,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,ties="rounded",interval.type="betaWald")
  $quantiles
    0.01  0.1  0.5 0.9 0.99
  y  0.1 0.15 0.75 2.4 8.68
  
  $CIs
  , , y
  
             0.01       0.1       0.5      0.9       0.99
  (lower 0.100000 0.1000000 0.3674054 1.352279  0.5674674
  upper) 1.132533 0.3464541 1.6303784 8.473337 10.2000000
  
  
  > 
  > oldsvyquantile(~y,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE)
  Statistic:
           y
  q0.01 0.10
  q0.1  0.20
  q0.5  0.80
  q0.9  2.40
  q0.99 8.68
  SE:
                 y
  q0.01 0.00000000
  q0.1  0.02611112
  q0.5  0.26542051
  q0.9  2.10908603
  q0.99 1.91777710
  > 
  > oldsvyquantile(~y,ddf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE,ties="rounded",interval.type="betaWald",df=Inf)
  $quantiles
    0.01  0.1  0.5 0.9 0.99
  y  0.1 0.15 0.75 2.4 8.68
  
  $CIs
  , , y
  
             0.01       0.1       0.5      0.9       0.99
  (lower 0.100000 0.1000000 0.3674054 1.352279  0.5674674
  upper) 1.132533 0.3464541 1.6303784 8.473337 10.2000000
  
  
  > 
  > oldsvyquantile(~y,rdf, c(0.01,0.1,0.5,0.9,0.99),ci=TRUE, df=Inf)
  Statistic:
           y
  q0.01 0.10
  q0.1  0.20
  q0.5  0.80
  q0.9  2.40
  q0.99 8.68
  SE:
                 y
  q0.01 0.00000000
  q0.1  0.02551067
  q0.5  0.25244616
  q0.9  2.23558589
  q0.99 2.03881114
  > 
  > proc.time()
     user  system elapsed 
     2.39    0.31    2.70 
Running the tests in 'tests/raowuboot.R' failed.
Complete output:
  > ## regression test for bug reported by Richard Valliant
  > 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
  
  > s<-subbootweights(c(1,1),1:2, 50)
  > stopifnot(all(s$repweights$weights %in% c(0,2)))
  > 
  > proc.time()
     user  system elapsed 
     1.40    0.17    1.56 
Running the tests in 'tests/regpredict.R' failed.
Complete output:
  > 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.31    0.20    1.50 
Running the tests in 'tests/svyby_bug.R' failed.
Complete output:
  > 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
  
  > options(warn=2)
  > 
  > ## Caused warnings and unhelpful results in 4.1_1 (Guilherme Jacob)
  > data(api)
  > dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > svyby(~api99, ~stype, dclus1, svymean )
    stype    api99       se
  E     E 607.7917 22.81660
  H     H 595.7143 41.76400
  M     M 608.6000 32.56064
  > 
  > set.seed(123)
  > apiclus1$api99[ sample.int( nrow(apiclus1) , 5 ) ] <- NA
  > dclus1.na <-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
  > 
  > # subsetting w/ na.rm = FALSE...
  > svymean( ~api99 , subset( dclus1.na , stype == "E" ) , na.rm = FALSE )
        mean SE
  api99   NA NA
  > svymean( ~api99 , subset( dclus1.na , stype == "H" ) , na.rm = FALSE )
          mean     SE
  api99 595.71 41.764
  > svymean( ~api99 , subset( dclus1.na , stype == "M" ) , na.rm = FALSE )
        mean SE
  api99   NA NA
  > 
  > # ... looks like this:
  > svyby(~api99, ~stype, dclus1.na , svymean )
    stype    api99     se
  E     E       NA     NA
  H     H 595.7143 41.764
  M     M       NA     NA
  > 
  > # subsetting w/ na.rm = TRUE...
  > svymean( ~api99 , subset( dclus1.na , stype == "E" ) , na.rm = TRUE )
          mean     SE
  api99 606.53 22.901
  > svymean( ~api99 , subset( dclus1.na , stype == "H" ) , na.rm = TRUE )
          mean     SE
  api99 595.71 41.764
  > svymean( ~api99 , subset( dclus1.na , stype == "M" ) , na.rm = TRUE )
          mean     SE
  api99 613.39 33.224
  > 
  > # ... looks like this
  > svyby(~api99, ~stype, dclus1.na , svymean , na.rm = TRUE )
    stype    api99       se
  E     E 606.5319 22.90084
  H     H 595.7143 41.76400
  M     M 613.3913 33.22358
  > 
  > # Without missing values, this works:
  > svyby(~api99, ~stype, dclus1 , svymean , na.rm = TRUE , covmat = TRUE )
    stype    api99       se
  E     E 607.7917 22.81660
  H     H 595.7143 41.76400
  M     M 608.6000 32.56064
  > 
  > # ... but this breaks!
  > svyby(~api99, ~stype, dclus1.na , svymean , na.rm = TRUE , covmat = TRUE )
    stype    api99       se
  E     E 606.5319 22.90084
  H     H 595.7143 41.76400
  M     M 613.3913 33.22358
  > 
  > # ... and i don't think this is the expected behavior
  > svyby( ~api99, ~stype, dclus1.na , svymean , na.rm.all = TRUE , covmat = TRUE )
    stype    api99     se
  E     E       NA     NA
  H     H 595.7143 41.764
  M     M       NA     NA
  > svyby( ~api99, ~stype, dclus1.na , svymean , na.rm.all = TRUE , na.rm = TRUE , covmat = TRUE )
    stype    api99       se
  E     E 606.5319 22.90084
  H     H 595.7143 41.76400
  M     M 613.3913 33.22358
  > 
  > 
  > ## <TL> Now some more as tests
  > svyby(~api99, ~stype, dclus1.na , svytotal , na.rm = TRUE , covmat = TRUE )
    stype     api99       se
  E     E 2894629.0 792809.8
  H     H  282283.9 101492.9
  M     M  477513.4 107555.5
  > svyby(~api99, ~stype, dclus1.na , svyratio , na.rm = TRUE , denominator=~api00, covmat = TRUE )
    stype api99/api00 se.api99/api00
  E     E   0.9367750    0.006365702
  H     H   0.9630485    0.010520785
  M     M   0.9682246    0.009390982
  > 
  > ff<-function(f,d,...,na.rm=TRUE) svyglm(f,d,...)
  > svyby(api99~1, ~stype, dclus1.na , ff , na.rm = TRUE , covmat = TRUE )
    stype (Intercept) se.(Intercept)
  E     E    606.5319       22.90084
  H     H    595.7143       41.76400
  M     M    613.3913       33.22358
  > 
  > proc.time()
     user  system elapsed 
     1.70    0.17    1.85 
Running the tests in 'tests/svyby_se.R' failed.
Complete output:
  > 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)
  > 
  > a<-svyby(~enroll,~stype,design=dstrat,svytotal,vartype=c("ci","se"))
  > b<-svyby(~enroll,~stype,design=dstrat,svytotal,vartype=c("se","ci"))
  > 
  > 
  > stopifnot(all.equal(SE(a),SE(b)))
  > 
  > proc.time()
     user  system elapsed 
     1.60    0.21    1.82 
Running the tests in 'tests/svyivreg.R' failed.
Complete output:
  > 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
  
  > library(AER)
  Loading required package: car
  Loading required package: carData
  Loading required package: lmtest
  Loading required package: zoo
  
  Attaching package: 'zoo'
  
  The following objects are masked from 'package:base':
  
      as.Date, as.Date.numeric
  
  Loading required package: sandwich
  > 
  > load("cigsw.rda")
  > 
  > des<-svydesign(id=~1, weights=~wt, data=cigsw)
  > m<-svyivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), design=des)
  > 
  > all.equal(as.vector(coef(m)), c(10.42009 , -1.588135, 0.6140887),tolerance=1e-6)
  [1] TRUE
  > all.equal(as.vector(SE(m)), c( 1.047699, .3394232,  .3614382 ),tolerance=1e-6)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.56    0.28    1.82 
Running the tests in 'tests/toy_example_for_postStratify.R' failed.
Complete output:
  > 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
  
  > 
  > # Dummy data for testing nonresponse adjustments
  > toy <- data.frame(id = 1:10,
  +                   dummy = 1,
  +                   class = c(rep("A", 5), rep("B", 5)),
  +                   responded = c(rep(TRUE, 8), FALSE, FALSE),
  +                   weight = rep(100, 10))
  > 
  > # With jackknife replicate weights
  > toy_repweights <- matrix(rep(1000/9, 100), nrow = 10)
  > diag(toy_repweights) <- rep(0, 10)
  > 
  > # Scramble up which person is in which jackknife group
  > toy_repweights <- toy_repweights[sample(1:10, size = 10), ]
  > 
  > toy_design <- svrepdesign(variables = toy[, 1:4],
  +                           weights = toy$weight,
  +                           repweights = toy_repweights,
  +                           type = "JK1",
  +                           scale = 0.9)
  > 
  > # Get the sum of the weights for the full sample
  > poptotals <- as.data.frame(svyby(formula = ~dummy,
  +                                  by = ~class,
  +                                  design = toy_design,
  +                                  FUN = survey::svytotal))
  > 
  > poptotals$Freq <- poptotals$dummy
  > poptotals$dummy <- NULL
  > poptotals$se <- NULL
  > 
  > # Adjust the weights of the responding sample to match to the full sample
  > adjusted <-  postStratify(design = subset(toy_design, responded),
  +                           strata = ~class,
  +                           population = poptotals)
  > 
  > # This works for the weights...
  > svyby(formula = ~dummy,
  +       by = ~class,
  +       design = adjusted,
  +       FUN = survey::svytotal)
    class dummy se
  A     A   500  0
  B     B   500  0
  > 
  > # ...and some, but not all, the replicate weights
  > stopifnot(all.equal(colSums(toy_design$repweights), colSums(adjusted$repweights)))
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.43    0.15    1.59 
Running the tests in 'tests/twophase.R' failed.
Complete output:
  > 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
  
  > 
  > ## two-phase simple random sampling.
  > data(pbc, package="survival")
  > pbc$id<-1:nrow(pbc)
  > pbc$randomized<-with(pbc, !is.na(trt) & trt>-9)
  > (d2pbc<-twophase(id=list(~id,~id), data=pbc, subset=~I(!randomized)))
  Two-phase sparse-matrix design:
   twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
      subset = subset, data = data, pps = pps)
  Phase 1:
  Independent Sampling design (with replacement)
  svydesign(ids = ~id)
  Phase 2:
  Independent Sampling design
  svydesign(ids = ~id, fpc = `*phase1*`)
  > m<-svymean(~bili, d2pbc)
  > all.equal(as.vector(coef(m)),with(pbc, mean(bili[!randomized])))
  [1] TRUE
  > all.equal(as.vector(SE(m)),
  +           with(pbc, sd(bili[!randomized])/sqrt(sum(!randomized))),
  +           tolerance=0.002)
  [1] TRUE
  > 
  > ## two-stage sampling as two-phase
  > data(mu284)
  > ii<-with(mu284, c(1:15, rep(1:5,n2[1:5]-3)))
  > mu284.1<-mu284[ii,]
  > mu284.1$id<-1:nrow(mu284.1)
  > mu284.1$sub<-rep(c(TRUE,FALSE),c(15,34-15))
  > dmu284<-svydesign(id=~id1+id2,fpc=~n1+n2, data=mu284)
  > ## first phase cluster sample, second phase stratified within cluster
  > (d2mu284<-twophase(id=list(~id1,~id),strata=list(NULL,~id1),
  +                    fpc=list(~n1,NULL),data=mu284.1,subset=~sub,method="approx"))
  Two-phase design: twophase(id = list(~id1, ~id), strata = list(NULL, ~id1), fpc = list(~n1, 
      NULL), data = mu284.1, subset = ~sub, method = "approx")
  Phase 1:
  1 - level Cluster Sampling design
  With (5) clusters.
  svydesign(ids = ~id1, fpc = ~n1)
  Phase 2:
  Stratified Independent Sampling design
  svydesign(ids = ~id, strata = ~id1, fpc = `*phase1*`)
  > (d22mu284<-twophase(id=list(~id1,~id),strata=list(NULL,~id1),
  +                    fpc=list(~n1,NULL),data=mu284.1,subset=~sub,method="full"))
  Two-phase sparse-matrix design:
   twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
      subset = subset, data = data, pps = pps)
  Phase 1:
  1 - level Cluster Sampling design
  With (5) clusters.
  svydesign(ids = ~id1, fpc = ~n1)
  Phase 2:
  Stratified Independent Sampling design
  svydesign(ids = ~id, strata = ~id1, fpc = `*phase1*`)
  > summary(d2mu284)
  Two-phase design: twophase(id = list(~id1, ~id), strata = list(NULL, ~id1), fpc = list(~n1, 
      NULL), data = mu284.1, subset = ~sub, method = "approx")
  Phase 1:
  1 - level Cluster Sampling design
  With (5) clusters.
  svydesign(ids = ~id1, fpc = ~n1)
  Probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0.1     0.1     0.1     0.1     0.1     0.1 
  Population size (PSUs): 50 
  Phase 2:
  Stratified Independent Sampling design
  svydesign(ids = ~id, strata = ~id1, fpc = `*phase1*`)
  Probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.3333  0.3750  0.4286  0.4674  0.6000  0.6000 
  Stratum Sizes: 
             19 31 45 47 50
  obs         3  3  3  3  3
  design.PSU  3  3  3  3  3
  actual.PSU  3  3  3  3  3
  Population stratum sizes (PSUs): 
  19 31 45 47 50 
   5  7  8  5  9 
  Data variables:
  [1] "id1" "n1"  "id2" "y1"  "n2"  "id"  "sub"
  > t1<-svytotal(~y1, dmu284)
  > t2<-svytotal(~y1, d2mu284)
  > t22<-svytotal(~y1,d22mu284)
  > m1<-svymean(~y1, dmu284)
  > m2<-svymean(~y1, d2mu284)
  > m22<-svymean(~y1, d22mu284)
  > all.equal(coef(t1),coef(t2))
  [1] TRUE
  > all.equal(coef(t1),coef(t22))
  [1] TRUE
  > all.equal(coef(m1),coef(m2))
  [1] TRUE
  > all.equal(coef(m1),coef(m22))
  [1] TRUE
  > all.equal(as.vector(SE(m1)),as.vector(SE(m2)))
  [1] TRUE
  > all.equal(as.vector(SE(m1)),as.vector(SE(m22)))
  [1] TRUE
  > all.equal(as.vector(SE(t1)),as.vector(SE(t2)))
  [1] TRUE
  > all.equal(as.vector(SE(t1)),as.vector(SE(t22)))
  [1] TRUE
  > 
  > ## case-cohort design
  > ##this example requires R 2.3.1 or later for cch and data.
  > library("survival")
  > data(nwtco, package="survival")
  > ## unstratified, equivalent to Lin & Ying (1993)
  > print(dcchs<-twophase(id=list(~seqno,~seqno), strata=list(NULL,~rel),
  +                  subset=~I(in.subcohort | rel), data=nwtco))
  Two-phase sparse-matrix design:
   twophase2(id = id, strata = strata, probs = probs, fpc = fpc, 
      subset = subset, data = data, pps = pps)
  Phase 1:
  Independent Sampling design (with replacement)
  svydesign(ids = ~seqno)
  Phase 2:
  Stratified Independent Sampling design
  svydesign(ids = ~seqno, strata = ~rel, fpc = `*phase1*`)
  > cch1<-svycoxph(Surv(edrel,rel)~factor(stage)+factor(histol)+I(age/12),
  +                design=dcchs)
  > dcchs2<-twophase(id=list(~seqno,~seqno), strata=list(NULL,~rel),
  +                  subset=~I(in.subcohort | rel), data=nwtco,method="approx")
  > cch1.2<-svycoxph(Surv(edrel,rel)~factor(stage)+factor(histol)+I(age/12),
  +                design=dcchs)
  > all.equal(coef(cch1),coef(cch1.2))
  [1] TRUE
  > all.equal(SE(cch1),SE(cch1.2))
  [1] TRUE
  > ## Using survival::cch 
  > subcoh <- nwtco$in.subcohort
  > selccoh <- with(nwtco, rel==1|subcoh==1)
  > ccoh.data <- nwtco[selccoh,]
  > ccoh.data$subcohort <- subcoh[selccoh]
  > cch2<-cch(Surv(edrel, rel) ~ factor(stage) + factor(histol) + I(age/12),
  +           data =ccoh.data, subcoh = ~subcohort, id=~seqno,
  +           cohort.size=4028, method="LinYing", robust=TRUE)
  > 
  > print(all.equal(as.vector(coef(cch1)),as.vector(coef(cch2))))
  [1] TRUE
  > ## cch has smaller variances by a factor of 1.0005 because
  > ## there is a (n/(n-1)) in the survey phase1 varianace
  > print(all.equal(as.vector(SE(cch1)), as.vector(SE(cch2)),tolerance=0.0006))
  [1] TRUE
  > 
  > 
  > ## bug report from Takahiro Tsuchiya for version 3.4
  > ## We used to not match Sarndal exactly, because our old phase-one
  > ## estimator had a small bias for finite populations
  > rei<-read.table(tmp<-textConnection(
  + "  id   N n.a h n.ah n.h   sub  y
  + 1   1 300  20 1   12   5  TRUE  1
  + 2   2 300  20 1   12   5  TRUE  2
  + 3   3 300  20 1   12   5  TRUE  3
  + 4   4 300  20 1   12   5  TRUE  4
  + 5   5 300  20 1   12   5  TRUE  5
  + 6   6 300  20 1   12   5 FALSE NA
  + 7   7 300  20 1   12   5 FALSE NA
  + 8   8 300  20 1   12   5 FALSE NA
  + 9   9 300  20 1   12   5 FALSE NA
  + 10 10 300  20 1   12   5 FALSE NA
  + 11 11 300  20 1   12   5 FALSE NA
  + 12 12 300  20 1   12   5 FALSE NA
  + 13 13 300  20 2    8   3  TRUE  6
  + 14 14 300  20 2    8   3  TRUE  7
  + 15 15 300  20 2    8   3  TRUE  8
  + 16 16 300  20 2    8   3 FALSE NA
  + 17 17 300  20 2    8   3 FALSE NA
  + 18 18 300  20 2    8   3 FALSE NA
  + 19 19 300  20 2    8   3 FALSE NA
  + 20 20 300  20 2    8   3 FALSE NA
  + "), header=TRUE)
  > close(tmp)
  > 
  > des.rei <- twophase(id=list(~id,~id), strata=list(NULL,~h),
  +                     fpc=list(~N,NULL), subset=~sub, data=rei, method="approx")
  > tot<- svytotal(~y, des.rei)
  > des.rei2 <- twophase(id=list(~id,~id), strata=list(NULL,~h),
  +                     fpc=list(~N,NULL), subset=~sub, data=rei)
  > tot2<- svytotal(~y, des.rei2)
  > 
  > ## based on Sarndal et al (9.4.14)
  > rei$w.ah <- rei$n.ah / rei$n.a
  > a.rei <- aggregate(rei, by=list(rei$h), mean, na.rm=TRUE)
  > a.rei$S.ysh <- tapply(rei$y, rei$h, var, na.rm=TRUE)
  > a.rei$y.u <- sum(a.rei$w.ah * a.rei$y)
  > V <- with(a.rei, sum(N * (N-1) * ((n.ah-1)/(n.a-1) - (n.h-1)/(N-1)) * w.ah * S.ysh / n.h))
  > V <- V + with(a.rei, sum(N * (N-n.a) * w.ah * (y - y.u)^2 / (n.a-1)))
  > 
  > a.rei$f.h<-with(a.rei, n.h/n.ah)
  > Vphase2<-with(a.rei, sum(N*N*w.ah^2* ((1-f.h)/n.h) *S.ysh))
  > 
  > a.rei$f<-with(a.rei, n.a/N)
  > a.rei$delta.h<-with(a.rei, (1/n.h)*(n.a-n.ah)/(n.a-1))
  > Vphase1<-with(a.rei, sum(N*N*((1-f)/n.a)*( w.ah*(1-delta.h)*S.ysh+ ((n.a)/(n.a-1))*w.ah*(y-y.u)^2)))
  > 
  > V
  [1] 36522.63
  > Vphase1
  [1] 24072.63
  > Vphase2
  [1] 12450
  > vcov(tot)
           y
  y 35911.05
  attr(,"phases")
  attr(,"phases")$phase1
           [,1]
  [1,] 23461.05
  
  attr(,"phases")$phase2
        y
  y 12450
  
  > vcov(tot2)
           [,1]
  [1,] 36522.63
  attr(,"phases")
  attr(,"phases")$phase1
           [,1]
  [1,] 24072.63
  
  attr(,"phases")$phase2
        [,1]
  [1,] 12450
  
  > ## phase 2 identical
  > all.equal(Vphase2,drop(attr(vcov(tot),"phases")$phase2))
  [1] TRUE
  > all.equal(Vphase2,drop(attr(vcov(tot2),"phases")$phase2))
  [1] TRUE
  > ## phase 1 differs by 2.6% for old twophase estimator
  > Vphase1/attr(vcov(tot),"phases")$phase1
           [,1]
  [1,] 1.026068
  > all.equal(Vphase1,as.vector(attr(vcov(tot2),"phases")$phase1))
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.29    0.35    2.59