R Under development (unstable) (2026-02-23 r89457 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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

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

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

> ##
> ## tests calibration using a working model 
> ##
> 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(nwtco)
>  nwtco$incc2<-as.logical(with(nwtco, ifelse(rel | instit==2,1,rbinom(nrow(nwtco),1,.1))))
> 
> dccs2<-twophase(id=list(~seqno,~seqno),strata=list(NULL,~interaction(rel,instit)),
+     data=nwtco, subset=~incc2)
>     
> impmodel<-svyglm(I(histol==2)~instit*stage+age, design=dccs2, family=quasibinomial())
> nwtco$histhat<-1+predict(impmodel,newdata=nwtco)
> 
> phase1model<-glm(rel~histhat*factor(stage),data=nwtco)
> 
> newmethod<- calibrate(dccs2, formula=phase1model,phase=2,calfun=cal.raking)
> 
> infs<-survey:::estfuns(phase1model)
> colnames(infs)<-paste0("h",1:ncol(infs))
> nwtco<-cbind(nwtco,infs)
> 
> dccs2h<- twophase(id=list(~seqno,~seqno),strata=list(NULL,~interaction(rel,instit)),
+     data=nwtco, subset=~incc2)
> oldmethod<-calibrate(dccs2h, ~h1+h2+h3+h4+h5+h6+h7+h8,calfun="raking",phase=2)
> 
> a<-svyglm(I(histol==2)~instit*stage+age, design=newmethod, family=quasibinomial())
> b<-svyglm(I(histol==2)~instit*stage+age, design=oldmethod, family=quasibinomial())
> 
> all.equal(coef(a),coef(b))
[1] TRUE
> all.equal(SE(a), SE(b),tol=1e-3)
[1] "Mean relative difference: 0.003138839"
> 
> proc.time()
   user  system elapsed 
   1.79    0.23    2.01