R notes for Thu Nov 11 10:54:31 2004
====================================

Stepwise fitting of models
--------------------------

table.4.3 <- crabs
names(table.4.3) <- c("C","S","W","Sa","Wt","Sa.bin")

table.4.3$C.fac<-factor(table.4.3$C, levels=c("5","4","3","2"))
table.4.3$S.fac<-factor(table.4.3$S, levels=c("3","2","1"))

crab.fit.logist.full<-glm(Sa.bin~C.fac+S.fac+W+I(Wt/1000), 
   family=binomial,data=table.4.3)

crab.fit.logist.stuffed<-glm(Sa.bin~C.fac*S.fac*W, family=binomial,
   data=table.4.3)

res<-step(crab.fit.logist.stuffed, list(lower = ~1, upper =  
   formula(crab.fit.logist.stuffed) ), scale=1, trace=T, direction="backward") 

res$anova

library(mass)
# or library(MASS)

res<-stepAIC(crab.fit.logist.stuffed, list(lower = ~1, 
   upper =  formula(crab.fit.logist.stuffed) ), scale=1, 
   trace=T, direction="backward") 

res$anova

# regression diagnostics here?


Common Odds Ratio and Cochrane-Mantel-Haenszel test
---------------------------------------------------


Data for multicenter clinical trial of cream vs placebo in curing an
infection, Table 6.9, p. 230, Agresti.

table.6.9 <- read.table("table6-9.txt",header=T)

table.6.9$Treatment <- as.factor(table.6.9$Treatment)
table.6.9$Response <- as.factor(table.6.9$Response)

levels(table.6.9$Treatment)<-c("Drug","Control")
levels(table.6.9$Response)<-c("Success","Failure")

# rewrite the data frame as a contingency table...

table.6.9.array<-xtabs(Freq~Treatment+Response+Center, data=table.6.9)

# Following function taken from Laura Thompson's notes...

oddsratio.L<-
function (x, stratum = NULL, Log = TRUE, conf.level = 0.95, correct=T) 
{
    l <- length(dim(x))
    if (l > 2 && is.null(stratum)) 
        stratum <- 3:l
    if (l - length(stratum) > 2) 
        stop("All but 2 dimensions must be specified as strata.")
    if (l == 2 && dim(x) != c(2, 2)) 
        stop("Not a 2 x 2 - table.")
    if (!is.null(stratum) && dim(x)[-stratum] != c(2, 2)) 
        stop("Need strata of 2 x 2 - tables.")
    lor <- function(y, correct, Log) {
        if(correct) y<-y + 0.5
        or <- y[1, 1] * y[2, 2]/y[1, 2]/y[2, 1]
        if (Log) 
            log(or)
        else or
    }
    ase <- function(y,correct) sqrt(sum(1/(ifelse(correct,y + 0.5,y))))
    if (is.null(stratum)) {
        LOR <- lor(x, correct)
        ASE <- ase(x)
    }
    else {
        LOR <- apply(x, stratum, lor, correct=correct, Log=Log)
        ASE <- apply(x, stratum, ase, correct=correct)
    }
    I <- ASE * qnorm((1 + conf.level)/2)
    Z <- LOR/ASE
    structure(LOR, ASE = if (Log) 
        ASE, lwr = if (Log) 
        LOR - I
    else exp(log(LOR) - I), upr = if (Log) 
        LOR + I
    else exp(log(LOR) + I), Z = if (Log) 
        Z, P = if (Log) 
        1 - pnorm(abs(Z)), log = Log, class = "oddsratio")
}

oddsratio.L(table.6.9.array, correct=F, Log=F)

length(unique(table.6.9$Center))

oddsratio.L(table.6.9.array, correct=F, Log=F)[1:8]

#################################
# CMH -- this is a function available in R or Splus

mantelhaen.test(table.6.9.array, correct=F)

mantelhaen.test(table.6.9.array, correct=F,exact=T)

####################################


# logit model test

# get counts for each tx within each center (16 counts)
n<-aggregate(table.6.9$Freq, 
     list(table.6.9$Treatment,table.6.9$Center), FUN=sum)$x

# attach them to the table, replicated within succ/fail columns
table.6.9$n<-rep(n,rep(2,16))

# now fit the model for the cond prob of success given tx and center,
# and compare with the null model of no common OR

res<-glm(Freq/n~Center+Treatment, family=binomial, 
     data=table.6.9,  weights=n, subset=Response=="Success")

summary(res)

anova(res)

anova(res,test="Chisq")

# heterogeneity of conditional ORs

res2<-update(res, .~.-Treatment)
res3<-update(res, .~.+Center:Treatment)

anova(res2,res3,test="Chisq")

# res has only a common tx x ctrl interaction
# res2 has no tx effect (a null model)
# res3 has center x tx x ctrl interaction

anova(res,res3,test="Chisq")

=======================================================================