sink(file="WNARout.txt",split=TRUE,append=TRUE)
load(system.file("doc","nwts.rda",package="survey")
nwtsc <- read.table("NWTSc_dat.txt",header=T)
nwtsb <- read.table("NWTSb_dat.txt",header=T)
nwtexb<-read.csv(file="nwtexb.txt")
nwtexb$in.ccs<-!is.na(nwtexb$histol)
describe(nwtexb)
xtabs(~instit+cc, data=nwtexb)
xtabs(~instit+cc, subset=!is.na(histol),data=nwtexb)
nwtex<-read.table("nwtex.txt",header=T)

# Check contents of nwtexb
nwtexb$in.ccs<-!is.na(nwtexb$histol)
nwtexb$id<-1:4088
with(nwtexb,table(instit,cc,in.ccs))

# Fit models using survey package
dccs2<-twophase(id=list(~id,~id),subset=~in.ccs,strata=list(NULL,~interaction(instit,cc)),data=nwtexb)
sum.dccs2<-summary(dccs2)
fit2<-svyglm(cc~stage*histol,family=binomial,design=dccs2)
sum2<-summary(fit2)
round(sum2$coefficients,4)
dccs8<-twophase(id=list(~id,~id),subset=~in.ccs,strata=list(NULL,~interaction(stage,instit,cc)),data=nwtexb)
summary(dccs8)
fit8<-svyglm(cc~stage*histol,family=binomial,design=dccs8)
sum8<-summary(fit8)
round(sum8$coefficients,4)
pop.types<-xtabs(~stage+instit+cc,data=nwtexb)
sample.types<-xtabs(~stage+instit+cc,subset=in.ccs,data=nwtexb)
dccsa8<-calibrate(dccs2,formula=~interaction(stage,instit,cc))
fita8<-svyglm(cc~stage*histol,family=binomial,design=dccsa8)
suma8<-summary(fita8)
round(suma8$coefficients,4)
dccse8<-estWeights(dccs2,formula=~interaction(stage,instit,cc))
fite8<-svyglm(cc~factor(stage)*factor(histol),family=binomial,design=dccse8)
summary(fite8)
dccsb8<-calibrate(dccs2,formula=~factor(stage))
fitb8<-svyglm(cc~factor(stage)*factor(histol),family=binomial,design=dccsb8)
summary(fitb8)
dccsf8<-estWeights(dccs2,formula=~factor(stage))
fitf8<-svyglm(cc~factor(stage)*factor(histol),family=binomial,design=dccsf8)
summary(fitf8)

# Construct expanded dataset for NWTS
nwtsn <- nwts
b <- rbind(nwts,nwtsn)
b$cc <- rep(c(1,0), each=16)
b$n <- ifelse(b$cc,b$case,b$control)
index <- rep(1:32, b$n)
nwtex <- b[index, c(1:3,6,7)]
nwtex$id <- 1:4088
nwtex<-nwtex[,c(1:4,6)]

# and then for NWTSc
nwtsnc <- nwtsc
nwtsnc$case <- nwts$case - nwtsc$case
nwtsnc$control <- nwts$control - nwtsc$control
a <- rbind(nwtsc,nwtsnc)
a$in.ccs <- rep(c(TRUE,FALSE), each=16)
b <- rbind(a,a)
b$cc <- rep(c(1,0), each=32)
b$n <- ifelse(b$cc,b$case,b$control)
index <- rep(1:64, b$n)
nwtexc <- b[index, c(1:3,6,7)]
nwtexc$id <- 1:4088
nwtexc$histol[nwtexc$in.ccs==FALSE]<-NA

# stage, instit and histol are all factors in all datasets
nwtex$stage<-factor(nwtex$stage,labels=c("I","II","III","IV"))
nwtexb$stage<-factor(nwtexb$stage,labels=c("I","II","III","IV"))
nwtexc$stage<-factor(nwtexc$stage,labels=c("I","II","III","IV"))
nwtsb$stage<-factor(nwtsb$stage,labels=c("I","II","III","IV"))

nwtex$instit<-factor(nwtex$instit,labels=c("FH","UH"))
nwtexb$instit<-factor(nwtexb$instit,labels=c("FH","UH"))
nwtexc$instit<-factor(nwtexc$instit,labels=c("FH","UH"))
nwtsb$instit<-factor(nwtsb$instit,labels=c("FH","UH"))

nwtex$histol<-factor(nwtex$histol,labels=c("FH","UH"))
nwtexb$histol<-factor(nwtexb$histol,labels=c("FH","UH"))
nwtexc$histol<-factor(nwtexc$histol,labels=c("FH","UH"))
nwtsb$histol<-factor(nwtsb$histol,labels=c("FH","UH"))

# fit interaction model to main cohort
fit.main<-glm(cc~stage*histol,family=binomial,data=nwtex)
sum.main<-summary(fit.main)
round(sum.main$coefficients,4)

# load tps
source("tps.R")
nn<-xtabs(~instit+cc,data=nwtex)
nn0<-nn[,1]
nn1<-nn[,2]
nwtex$group<-interaction(nwtex$instit,nwtex$stage)
nns<-xtabs(~group+cc,data=nwtex)
nn0s<-nns[,1]
nn1s<-nns[,2]
nwtexb$group<-interaction(nwtexb$instit,nwtexb$stage)
nwtexba<-nwtexb[nwtexb$in.ccs,]

# WL analyses using tps
fitwl2<-tps(cc~stage*histol,data=nwtexba,nn0=nn0,nn1=nn1, group=nwtexba$instit, method="WL")
sumwl2<-summary(fitwl2)
round(sumwl2$coefficients,4)

fitwl8<-tps(cc~stage*histol,data=nwtexba,nn0=nn0s,nn1=nn1s, group=nwtexba$group, method="WL")
sumwl8<-summary(fitwl8)
round(sumwl8$coefficients,4)


# ML analyses using tps
# Must unfactor the factors for this
nwtexba$stage<-as.numeric(nwtexba$stage)
nwtexba$histol<-as.numeric(nwtexba$histol)
nwtexba$instit<-as.numeric(nwtexba$instit)
nwtexba$group<-as.numeric(nwtexba$group)


fitml2<-tps(cc~factor(stage)*factor(histol),data=nwtexba,nn0=nn0,nn1=nn1, group=nwtexba$instit, method="ML")
summl2<-summary(fitml2)
round(summl2$coefficients,4)

fitml8<-tps(cc~factor(stage)*factor(histol),data=nwtexba,nn0=nn0s,nn1=nn1s, group=nwtexba$group, method="ML")
summl8<-summary(fitml8)
round(summl8$coefficients,4)

fitml2<-
tps(cbind(case,control)~factor(stage)*factor(histol),
data=nwtsb,nn0=nn0,nn1=nn1, group=nwtsb$instit, method="ML")
summl2<-summary(fitml2)
round(summl2$coefficients,4)