###################        
#### Chapter 2 ####
###################


# =============================
## Comparing two proportions ##
# =============================

# Example (heart attack v. no attack)
x<-c(104,189) # aspirin, placebo
n<-c(11037,11034)

# test H0:p1=p2 (equal probabilities of heart attack)
prop.test(x,n,p=c(.5,.5))
prop.test(x,n)$p.value

# Test H0:p1=p2 v. H1:p1<p2
prop.test(x,n,alt="less")$p.value

# sample difference of proportions:
temp<-prop.test(x,n)  
names(temp$estimate)<-NULL  # optional
temp$estimate[[1]]-temp$estimate[[2]]

# relative risk
temp$estimate[[2]]/temp$estimate[[1]]

# odds ratio
x[2]*(n[1]-x[1])/(x[1]*(n[2]-x[2]))


# ==================
## partial tables ##
# ==================

# Create a cross-classified table out of a Table
vic.race<-c("white","black")
def.race<-vic.race
death.penalty<-c("yes", "no")
datalabel<-list(defendant=def.race,death=death.penalty,victim=vic.race)
table.2.6<- expand.grid(defendant=def.race,death=death.penalty,victim=vic.race) 
Data<-c(53, 11, 414, 37, 0, 4, 16, 139)
table.2.6<-cbind(table.2.6,count=Data)
temp <- xtabs(count~defendant+death+victim ,data=table.2.6)
temp

# conditional odds ratios
apply(temp,3,function(x) x[1,1]*x[2,2]/(x[2,1]*x[1,2]))
summary(oddsratio(temp, log=F, stratum=3))
summary.table(temp)



# ========================
# association in tables
# ========================

income<-c("<15000","15000-25000","25000-40000",">40000")
jobsat<-c("VD","LD","MS","VS")
table.2.8<-expand.grid(income=income,jobsat=jobsat) 
Data<-c(1,2,1,0,3,3,6,1,10,10,14,9,6,7,12,11)
table.2.8<-cbind(table.2.8,count=Data)
(temp<-xtabs(count~income+jobsat,table.2.8))
library(vcd)
assocstats(temp)


# ========================
## Confidence intervals ##
# =========================


library(methods)
Drug<-c("Placebo","Aspirin")
Infarction<-c("yes","no")
table.3.1<-expand.grid(drug=Drug,infarction=Infarction) 
Data<-c(28,18,656,658)
table.3.1<-cbind(table.3.1,count=Data)
tapply(table.3.1$count,table.3.1[,1:2], sum)

# wald confidence interval
Wald.ci<-function(Table, aff.response, alpha=.05){
    # Gives two-sided Wald CI's for odds ratio, difference in proportions and relative risk.
    # Table is a 2x2 table of counts with rows giving the treatment populations
    # aff.response is a string like "c(1,1)" giving the cell of the beneficial response and the treatment category 
    # alpha is significance level
    
    pow<-function(x, a=-1) x^a 
    z.alpha<-qnorm(1-alpha/2)
    
    if(is.character(aff.response)) where<-eval(parse(text=aff.response)) else where<-aff.response
    Next<-as.numeric(where==1) + 1
    
    # OR
    odds.ratio<-Table[where[1],where[2]]*Table[Next[1],Next[2]]/(Table[where[1],Next[2]]*Table[Next[1],where[2]])
    se.OR<-sqrt(sum(pow(Table)))
    ci.OR<-exp(log(odds.ratio) + c(-1,1)*z.alpha*se.OR)
    
    # difference of proportions 
    p1<-Table[where[1],where[2]]/(n1<-Table[where[1],Next[2]] + Table[where[1],where[2]])
    p2<-Table[Next[1],where[2]]/(n2<-Table[Next[1],where[2]]+Table[Next[1],Next[2]])
    
    se.diff<-sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
    ci.diff<-(p1-p2) + c(-1,1)*z.alpha*se.diff
    
    # relative risk
    RR<-p1/p2
    se.RR<-sqrt((1-p1)/(p1*n1) + (1-p2)/(p2*n2))
    ci.RR<-exp(log(RR) + c(-1,1)*z.alpha*se.RR)
    
    list(OR=list(odds.ratio=odds.ratio, CI=ci.OR), proportion.difference=list(diff=p1-p2, CI=ci.diff), relative.risk=list(relative.risk=RR,CI=ci.RR))
    }
Wald.ci(temp, "c(1, 1)")        # or use Wald.ci(temp, c(1, 1))



# =========================
## Tests of Independence ##
# =========================


# Table 3.2 (p.80)

# set up data and expected counts
religion.counts<-c(178,138,108,570,648,442,138,252,252)
table.3.2<-cbind(expand.grid(list(Religious.Beliefs=c("Fund", "Mod", "Lib"), Highest.Degree=c("<HS","HS or JH", "Bachelor or Grad"))),count=religion.counts)
table.3.2.array<-tapply(table.3.2$count,table.3.2[,1:2], sum)


# Pearson Chi-squared test
(res<-chisq.test(table.3.2.array))
res$expected
#                  Highest.Degree
# Religious.Beliefs      <HS HS or JH Bachelor or Grad
#              Fund 137.8078 539.5304         208.6618
#              Mod  161.4497 632.0910         244.4593
#              Lib  124.7425 488.3786         188.8789
(res<-chisq.test(table.3.2, sim=T)) # simulate p-value



# LRT
2*sum(table.3.2.array*log(table.3.2.array/res$expected))
# or
table.3.2.df<-cbind(expand.grid(list(Religious.Beliefs=c("Fund", "Mod", "Lib"), Highest.Degree=c("<HS","HS or JH", "Bachelor or Grad"))),count=religion.counts)
res<-glm(count~Religious.Beliefs+Highest.Degree,data=table.3.2.df,family=poisson)
predict(res,type="response")


# using vcd
library(vcd)
summary(assoc.stats(table.3.2.array))
expected(table.3.2.array)
mar.table(table.3.2.array)


## Follow -up tests ##
# pearson residuals (same as splus)
# partitioning chi-square (same as splus)

## Two-way Table with Ordered Classifications ##
# linear trend (same as splus)
# monotone trend
levels(table.2.8$income)<-c(7.5,20,32.5,60)
levels(table.2.8$jobsat)<-1:4

Gamma2.f<-function(x, pr=0.95)
{
    # x is a matrix of counts.  You can use output of crosstabs or xtabs in R.
    # A matrix of counts can be formed from a data frame by using design.table.
    
    # Confidence interval calculation and output from Greg Rodd
    
    # Check for using S-PLUS and output is from crosstabs (needs >= S-PLUS 6.0)
    if(is.null(version$language) && inherits(x, "crosstabs")) { oldClass(x)<-NULL; attr(x, "marginals")<-NULL}
        
    n <- nrow(x)
    m <- ncol(x)
    pi.c<-pi.d<-matrix(0,nr=n,nc=m)
        
    row.x<-row(x)
    col.x<-col(x)
    
    for(i in 1:(n)){
        for(j in 1:(m)){
            pi.c[i, j]<-sum(x[row.x<i & col.x<j]) + sum(x[row.x>i & col.x>j])
            pi.d[i, j]<-sum(x[row.x<i & col.x>j]) + sum(x[row.x>i & col.x<j])
        }
    }

    C <- sum(pi.c*x)/2
    D <- sum(pi.d*x)/2
    
    psi<-2*(D*pi.c-C*pi.d)/(C+D)^2
    sigma2<-sum(x*psi^2)-sum(x*psi)^2
    
    gamma <- (C - D)/(C + D)
    pr2 <- 1 - (1 - pr)/2
    CIa <- qnorm(pr2) * sqrt(sigma2) * c(-1, 1) + gamma

    list(gamma = gamma, C = C, D = D, sigma = sqrt(sigma2), Level = paste(
        100 * pr2, "%", sep = ""), CI = paste(c("[", max(CIa[1], -1), 
        ", ", min(CIa[2], 1), "]"), collapse = ""))     
}

temp<-xtabs(formula = count ~ income + jobsat, data = table.2.8)
res<-Gamma2.f(temp)


# =====================================
## Small samples tests of independence
# =====================================


# Fisher's Exact test

# Table 3.8 (p.92)
# fisher.test() gives a two-sided test
(fisher.test(matrix(c(3,1,1,3),byrow=T,ncol=2), alternative="greater"))

#         Fisher's Exact Test for Count Data
# 
# data:  matrix(c(3, 1, 1, 3), byrow = T, ncol = 2) 
# p-value = 0.2429
# alternative hypothesis: true odds ratio is greater than 1 
# 95 percent confidence interval:
#  0.3135693       Inf 
# sample estimates:
# odds ratio 
#   6.408309 

# ====================================
# Fisher's Exact test for I x J tables
# ====================================

table.3.9<-matrix(c(25,25,12,0,1,3),byrow=T,ncol=3)

Gamma2.f(table.3.9)
# $gamma:
# [1] 0.8716578

# $C:
# [1] 175

# $D:
# [1] 12


library(combinat)
num<-prod(fact(rep(1,2)%*%table.3.9))*prod(fact(rep(1,3)%*%t(table.3.9)))
den<-fact(sum(table.3.9))*prod(fact(table.3.9))
term1<-num/den
temp<-matrix(c(25,26,11,0,0,4),byrow=T,ncol=3)
num<-prod(fact(rep(1,2)%*%temp))*prod(fact(rep(1,3)%*%t(temp)))
den<-fact(sum(temp))*prod(fact(temp))
term2<-num/den

term1+term2
# [1] 0.01830808

fisher.test(table.3.9)$p.value/2
# [1] 0.01704545



# ===================================
## Conditional associations -CMH Test
# ===================================

table.6.9<-data.frame(scan(file="c:/program files/insightful/splus61/users/CDA/clinical trials table 69.ssc", what=list(Center="",Treatment="",Response="",Freq=0)))
 a 1 1 11 
 a 1 2 25  
 a 2 1 10  
 a 2 2 27
 b 1 1 16 
 b 1 2 4   
 b 2 1 22  
 b 2 2 10
 c 1 1 14 
 c 1 2 5   
 c 2 1 7   
 c 2 2 12
 d 1 1 2  
 d 1 2 14  
 d 2 1 1   
 d 2 2 16
 e 1 1 6  
 e 1 2 11  
 e 2 1 0   
 e 2 2 12
 f 1 1 1  
 f 1 2 10  
 f 2 1 0   
 f 2 2 10
 g 1 1 1  
 g 1 2 4   
 g 2 1 1   
 g 2 2 8
 h 1 1 4  
 h 1 2 2   
 h 2 1 6   
 h 2 2 1


levels(table.6.9$Treatment)<-c("Drug","Control")
levels(table.6.9$Response)<-c("Success","Failure")
table.6.9.array<-xtabs(Freq~Treatment+Response+Center, data=table.6.9)

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)

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



# ============================================
# McNemar Test - Compare Dependent Proportions
# ============================================

## Comparing dependent proportions
table.10.1<-matrix(c(794,150,86,570),byrow=T,ncol=2)
mcnemar.test(table.10.1,correct=F)











#