pdf("g03.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class553")
setwd("~/Taught1/960-553/Data")
<pre>
#<a name=block1> Block 1</a>
#*****************************/
#* Proportion Difference CI. */
#*****************************/
library(DescTools)#For BinomDiffCI
#Give the traditional Wald interval.
#*************************************************************/
#*  Prostate data                                            */
#* From Andrews and Herzberg (1985), Data: A Collection of   */
#* Problems from Many Fields for the Student and Research    */
#* Worker, Table 46.  Observations represent subjects in a   */
#* prostate cohort study, randomized to one of four dose     */
#* levels of diethylstilbestrol.  Rx records dose in four    */
#* ordered categories, with 1 being placebo.  Disease stage  */
#* is 3 or 4.  monfol is months of followup.  surv is 0 if   */
#* alive after 50 mo, and codes cause of death otherwise.    */
#* http://lib.stat.cmu.edu/datasets/Andrews/T46.1            */
#* The on-line version of the data set adds 3 fields before  */
#* the first field in the book.  Variables of interest are   */
#* stage, rx, monfol, and surv in fields 5, 6, 10, 11 of the */
#* online version, resp.  Causes of death are given by var-  */
#* ious positive integers in surv; I recode these to 1.  The */
#* data file has more columns than we are interested in.  Put*/
#* place-holding variables in for variables not of interest  */
#* between variables of interest.  Data were previously pub- */
#* lished by Byar and Corle (1977 Chronic Disease) and Byar  */
#* and Green (1980 Bull. Cancer).  Lower value of dichoto-   */
#* mized dose begins with blank to make it alphabetized      */
#* before high.                                              */
#*************************************************************/
# If we omit the ones past where the data of interest stops
# in the "what" list, and use flush=T, R will ignore them. 
prostate<-read.table("T46.1")[,c(5,6,10,11)]
dimnames(prostate)[[2]]<-c("stage","rx","monfol","surv")
prostate$alive<-(prostate$surv==0)+0                      
prostate$dose<-c(" low","high")[1+(prostate$rx>1)]         
stage4<-prostate[prostate$stage==4,]
mat<-table(stage4$rx,stage4$alive)
BinomDiffCI( mat[1,1],mat[1,1]+mat[1,2],
    mat[2,1],mat[2,1]+mat[2,2],method="wald")
#Be very careful how you use the above command.  It
#gives you an interval for the first proportion minus
#the second, which I would have found the reverse to
#be more natural.
#Give Newcombe's interval with continuity correction.
BinomDiffCI( mat[1,1],mat[1,1]+mat[1,2],
    mat[2,1],mat[2,1]+mat[2,2],method="scorecc")
library(TheorStat)#For bivariatebincover
bivariatebincover()
#<a name=block2> Block 2</a>
#************************/
#* Odds ratio CI.       */
#************************/
# fisher.test does exact test and gives exact CI.
fisher.test(stage4$dose,stage4$alive)$conf.int
mat<-table(stage4$dose,stage4$alive)
# I don't think that there's a built-in that gives the 
# asymptotic confidence interval. Get interval on log scale.
(ci<-log((mat[1,1]*mat[2,2])/(mat[2,1]*mat[1,2]))+
   c(-1,1)*1.96*sqrt(sum(1/mat)))
# Exponentiate to get on original scale
exp(ci)
</pre>