################# ## Section 4.1 ## ################# ## library("osDesign") data("Ohio") Ohio ################# ## Section 4.2 ## ################# ## XM <- cbind(Int=1, Ohio[,1:3]) XI <- cbind(XM, SbyR=XM[,3]*XM[,4]) XM ## cbind(Int=1, Age1=as.numeric(Ohio[,1] == 1), Age2=as.numeric(Ohio[,1] == 2), Ohio[,2:3]) ## fitM <- glm(cbind(Death, N-Death) ~ factor(Age) + Sex + Race, data=Ohio, family=binomial) summary(fitM) ## fitI <- glm(cbind(Death, N-Death) ~ factor(Age) + Sex * Race, data=Ohio, family=binomial) summary(fitI) ################# ## Section 5.1 ## ################# ## ocAge <- tpsSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, nII0=c(50,50,50), nII1=c(50,50,50), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) ocAge ################# ## Section 5.2 ## ################# ## ## Percent Median Bias ## ocInterCC <- ccSim(B=50000, betaTruth=fitI$coef, X=XI, N=Ohio$N, nCC=200, r=1, betaNames=c("Int", "Age1", "Age2", "Sex", "Race", "SxR"), monitor=1000) ## ocInterS <- tpsSim(B=50000, betaTruth=fitI$coef, X=XI, N=Ohio$N, strata=3, nII0=c(50,50), nII1=c(50,50), betaNames=c("Int", "Age1", "Age2", "Sex", "Race", "SxR"), monitor=1000) ## ocInterR <- tpsSim(B=50000, betaTruth=fitI$coef, X=XI, N=Ohio$N, strata=4, nII0=c(50,50), nII1=c(50,50), betaNames=c("Int", "Age1", "Age2", "Sex", "Race", "SxR"), monitor=1000) ## ocInterSR <- tpsSim(B=50000, betaTruth=fitI$coef, X=XI, N=Ohio$N, strata=c(3,4), nII0=c(25,25,25,25), nII1=c(25,25,25,25), betaNames=c("Int", "Age1", "Age2", "Sex", "Race", "SxR"), monitor=1000) ## resPB <- c(NA, NA, ocInterCC$results$betaMedianPB[1,4], NA, NA, ocInterCC$results$betaMedianPB[1,6]) resPB <- rbind(resPB, c(NA, NA, ocInterCC$results$betaMedianPB[2,4], NA, NA, ocInterCC$results$betaMedianPB[2,6])) resPB <- rbind(resPB, cbind(rbind(ocInterS$results$betaMedianPB[3:5,4], ocInterR$results$betaMedianPB[3:5,4], ocInterSR$results$betaMedianPB[3:5,4]), rbind(ocInterS$results$betaMedianPB[3:5,6], ocInterR$results$betaMedianPB[3:5,6], ocInterSR$results$betaMedianPB[3:5,6]))) round(resPB, 1) ################# ## Section 5.3 ## ################# ## ocAll <- tpsSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=0, nII=c(250, 250), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) ocAll$digits <- 1 ocAll ################# ## Section 5.4 ## ################# ## ocAgeCC <- tpsSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, cohort=FALSE, NI=c(5000,5000), nII0=c(50,50,50), nII1=c(50,50,50), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) ocAgeCC$digits <- 1 ocAgeCC ################# ## Section 5.5 ## ################# ## resultsCC <- ccSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, nCC=500, r=c(0.67, 1.5), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) resultsCC$digits <- 1 resultsCC ################# ## Section 6.1 ## ################# ## phaseI(betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4) phaseI(betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=c(2,4)) phaseI(betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, nII0=c(100, 100, 100), nII1=c(100, 100, 100)) phaseI(betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4, cohort=FALSE, NI=c(5000,5000)) ################# ## Section 6.2 ## ################# ## powerRaceTPS <- tpsPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4, nII=seq(from=100, to=1000, by=100), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) ## pdf("FigurePower_TPS_Race.pdf", width=10, height=10) par(mfrow=c(2,2)) plotPower(powerRaceTPS, coefNum=2, xAxis=seq(from=100, to=1000, by=100), main=expression("Age effect (65-74 vs. 55-64 years), " * beta[A1]), legendXY=c(800, 60)) plotPower(powerRaceTPS, coefNum=3, xAxis=seq(from=100, to=1000, by=100), main=expression("Age effect (75-84 vs. 55-64 years), " * beta[A2]), legendXY=c(800, 60)) plotPower(powerRaceTPS, coefNum=4, xAxis=seq(from=100, to=1000, by=100), main=expression("Sex effect, " * beta[S]), legendXY=c(800, 60)) plotPower(powerRaceTPS, coefNum=5, xAxis=seq(from=100, to=1000, by=100), main=expression("Race effect, " * beta[R]), legendXY=c(800, 60)) dev.off() ################# ## Section 6.3 ## ################# ## powerSexTPS <- tpsPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=3, nII=seq(from=100, to=1000, by=100), monitor=1000) nLvls <- length(powerSexTPS$nII) ## pdf("FigurePower_TPS_Compare.pdf", width=12, height=12) par(mfrow=c(2,2)) plotPower(powerRaceTPS, coefNum=2, include="ML", xAxis=seq(from=100, to=1000, by=100), main=expression("Age effect (65-74 vs. 55-64 years), " * beta[A1])) lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (3*nLvls) + c(1:nLvls)), 2], lwd=2, lty=2, pch=2, type="o") lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (0*nLvls) + c(1:nLvls)), 2], lwd=2, lty=3, pch=3, type="o") legend(650, 65, c("TPS: Race", "TPS: Sex", "Case-control"), lwd=c(2,2,2), lty=c(1,2,3), pch=c(1,2,3)) plotPower(powerRaceTPS, coefNum=3, include="ML", xAxis=seq(from=100, to=1000, by=100), main=expression("Age effect (75-84 vs. 55-64 years), " * beta[A2])) lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (3*nLvls) + c(1:nLvls)), 3], lwd=2, lty=2, pch=2, type="o") lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (0*nLvls) + c(1:nLvls)), 3], lwd=2, lty=3, pch=3, type="o") legend(650, 65, c("TPS: Race", "TPS: Sex", "Case-control"), lwd=c(2,2,2), lty=c(1,2,3), pch=c(1,2,3)) plotPower(powerRaceTPS, coefNum=4, include="ML", xAxis=seq(from=100, to=1000, by=100), main=expression("Sex effect, " * beta[S])) lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (3*nLvls) + c(1:nLvls)), 4], lwd=2, lty=2, pch=2, type="o") lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (0*nLvls) + c(1:nLvls)), 4], lwd=2, lty=3, pch=3, type="o") legend(650, 65, c("TPS: Race", "TPS: Sex", "Case-control"), lwd=c(2,2,2), lty=c(1,2,3), pch=c(1,2,3)) plotPower(powerRaceTPS, coefNum=5, include="ML", xAxis=seq(from=100, to=1000, by=100), main=expression("Race effect, " * beta[R])) lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (3*nLvls) + c(1:nLvls)), 5], lwd=2, lty=2, pch=2, type="o") lines(powerSexTPS$nII, powerSexTPS$betaPower[(1 + (0*nLvls) + c(1:nLvls)), 5], lwd=2, lty=3, pch=3, type="o") legend(650, 65, c("TPS: Race", "TPS: Sex", "Case-control"), lwd=c(2,2,2), lty=c(1,2,3), pch=c(1,2,3)) dev.off() ################# ## Section 6.4 ## ################# ## powerCC <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1, nCC=seq(from=100, to=500, by=50), betaNames=c("Int", "Age1", "Age2", "Sex", "Race"), monitor=1000) ## pdf("FigurePower_CC.pdf", width=10, height=10) par(mfrow=c(2,2)) plotPower(powerCC, coefNum=2, xAxis=seq(from=100, to=500, by=100), main=expression("Age effect (65-74 vs. 55-64 years), " * beta[A1])) plotPower(powerCC, coefNum=3, xAxis=seq(from=100, to=500, by=100), main=expression("Age effect (75-84 vs. 55-64 years), " * beta[A2])) plotPower(powerCC, coefNum=4, xAxis=seq(from=100, to=500, by=100), main=expression("Sex effect, " * beta[S])) plotPower(powerCC, coefNum=5, xAxis=seq(from=100, to=500, by=100), main=expression("Race effect, " * beta[R])) dev.off() ## powerCC100 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1, nCC=200, monitor=1000) powerCC150 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1.5, nCC=250, monitor=1000) powerCC200 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=2, nCC=300, monitor=1000) powerCC250 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=2.5, nCC=350, monitor=1000) powerCC300 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=3, nCC=400, monitor=1000) powerCC350 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=3.5, nCC=450, monitor=1000) powerCC400 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=4, nCC=500, monitor=1000) powerCC450 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=4.5, nCC=550, monitor=1000) powerCC500 <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=5, nCC=600, monitor=1000) ## powerCCr <- rbind(powerCC100$betaPower[2,-1], powerCC150$betaPower[2,-1], powerCC200$betaPower[2,-1], powerCC250$betaPower[2,-1], powerCC300$betaPower[2,-1], powerCC350$betaPower[2,-1], powerCC400$betaPower[2,-1], powerCC450$betaPower[2,-1], powerCC500$betaPower[2,-1]) ## pdf("FigurePower_CCr.pdf", width=11, height=11) par(cex=1.5) plot(c(1,5), c(0,100), xlab="Control:Case ratio, r", ylab="Power", axes=F, type="n") axis(1, at=seq(from=1, to=5, by=1)) axis(2, at=seq(from=0, to=100, by=20)) lines(seq(from=1, to=5, by=0.5), powerCCr[,1], lwd=2, type="o", pch=1) lines(seq(from=1, to=5, by=0.5), powerCCr[,2], lwd=2, type="o", pch=2) lines(seq(from=1, to=5, by=0.5), powerCCr[,3], lwd=2, type="o", pch=3) lines(seq(from=1, to=5, by=0.5), powerCCr[,4], lwd=2, type="o", pch=4) legend(2, 90, c(expression("Age effect (65-74 vs. 55-64 years), " * beta[A1]), expression("Age effect (75-84 vs. 55-64 years), " * beta[A2]), expression("Sex effect, " * beta[S]), expression("Race effect, " * beta[R])), lwd=rep(2,4), pch=1:4) dev.off() ################# ## Section 6.5 ## ################# ## newBetaM <- fitM$coef newBetaM[2:3] <- newBetaM[2:3] * 0.5 phaseI(betaTruth=newBetaM, X=XM, N=Ohio$N, strata=c(2,4)) ## newBetaM[1] newBetaM[1] <- beta0(betaX=newBetaM[-1], X=XM, N=Ohio$N, rhoY=sum(Ohio$Death)/sum(Ohio$N)) newBetaM[1] ## Timing ## ## calcTimeDiff <- function(start, end, units="s") { startD <- as.numeric(substring(start, 9, 10)) startH <- as.numeric(substring(start, 12, 13)) startM <- as.numeric(substring(start, 15, 16)) startS <- as.numeric(substring(start, 18, 19)) endD <- as.numeric(substring(end, 9, 10)) endH <- as.numeric(substring(end, 12, 13)) endM <- as.numeric(substring(end, 15, 16)) endS <- as.numeric(substring(end, 18, 19)) value <- (endD*86400 + endH*3600 + endM*60 + endS) - (startD*86400 + startH*3600 + startM*60 + startS) if(units == "m") value <- value / 60 if(units == "h") value <- value / 60^2 return(value) } ## ## start <- date() ocAge <- tpsSim(B=500, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, nII0=c(50,50,50), nII1=c(50,50,50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() ocAge <- tpsSim(B=1000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, nII0=c(50,50,50), nII1=c(50,50,50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() ocAge <- tpsSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=2, nII0=c(50,50,50), nII1=c(50,50,50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## ## start <- date() ocAll <- tpsSim(B=500, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=0, nII=c(250, 250), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() ocAll <- tpsSim(B=1000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=0, nII=c(250, 250), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() ocAll <- tpsSim(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=0, nII=c(250, 250), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## ## start <- date() powerRaceTPS <- tpsPower(B=500, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4, nII=seq(from=100, to=1000, by=100), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() powerRaceTPS <- tpsPower(B=1000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4, nII=seq(from=100, to=1000, by=100), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() powerRaceTPS <- tpsPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, strata=4, nII=seq(from=100, to=1000, by=100), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## ## start <- date() powerCC <- ccPower(B=500, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1, nCC=seq(from=100, to=500, by=50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() powerCC <- ccPower(B=1000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1, nCC=seq(from=100, to=500, by=50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m") ## start <- date() powerCC <- ccPower(B=10000, betaTruth=fitM$coef, X=XM, N=Ohio$N, r=1, nCC=seq(from=100, to=500, by=50), monitor=1000) end <- date() calcTimeDiff(start, date(), units="m")