## confidence intrervals ##


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))

# score confidence interval 

gfun<-function(p,y0){
    sum((c(-(28/p[1])+(656/(1-p[1]))-p[3], -(18/p[2])+(658/(1-p[2]))+p[3], p[1]-p[2])-y0)^2)
}
optim(fn=gfun,par=c(28/(28+656),18/(18+658),7), method="BFGS",y0=c(0,0,-.1),control=list(trace=2,REPORT=10))

score.ci<-function(Delta) {
    temp<-optim(fn=gfun,par=c(28/(28+656),18/(18+658),7), method="BFGS",y0=c(0,0,Delta), 
    control=list(ndeps=c(rep(.000000001,3))))
    p<-temp$par[1:2]
    z<-(p[1]-p[2]-Delta)/sqrt(p[1]*(1-p[1])/684 + p[2]*(1-p[2])/676)
}
Delta<-seq(-.2,.2,.0005)
Delta<-seq(-.05,.05,.0001)
z<-sapply(Delta,score.ci)
Delta[abs(z)< qnorm(.975)]

# profile likelihood ci

library(Bhat)

# neg. log-likelihood of "new" multinomial model
nlogf <- function (p) { 
    p1p<-p[1]
    pp1<-p[2]
    theta<-p[3]
    n11 <- table.3.1$count[1]
    n21<-table.3.1$count[2]
    n12<-table.3.1$count[3]
    n22<-table.3.1$count[4]
    p22<-(-1 + p1p + pp1 + 2*theta - (p1p + pp1)*theta - sqrt(-4*p1p*(-1 + pp1)*(-1 + theta) + 
        (1 + p1p + pp1*(-1 + theta) - p1p*theta)^2))/(2*(-1 + theta))
    p11 <- (1 + p1p*(-1 + theta) + pp1*(-1 + theta) - sqrt(-4*p1p*(-1 + pp1)*(-1 + theta) + 
        (1 + p1p + pp1*(-1 + theta) - p1p*theta)^2))/(2*(-1 + theta))
    p21 <- (-1 + p1p + pp1*(-1 + theta) - p1p*theta + sqrt(-4*p1p*(-1 + pp1)*(-1 + theta) + 
        (1 + p1p + pp1*(-1 + theta) - p1p*theta)^2))/(2*(-1 + theta))
    p12 <- (-1 + pp1 + p1p*(-1 + theta) - pp1*theta + sqrt(-4*p1p*(-1 + pp1)*(-1 + theta) + 
        (1 + p1p + pp1*(-1 + theta) - p1p*theta)^2))/(2*(-1 + theta))
                
    -(n11*log(p11) + n12*log(p12) + n21*log(p21) + n22*log(p22))
    }

# check answer with MLEs
nlogf(c((28+656)/(684+676),(28+18)/(684+676),28*658/(18*656)))

# neg LLH of multinomial
logLH<-function(p){
    n11 <- table.3.1$count[1]
    n21<-table.3.1$count[2]
    n12<-table.3.1$count[3]
    n22<-table.3.1$count[4]
    p11<-p[1]
    p12<-p[2]
    p21<-p[3]
-(n11*log(p11) + n12*log(p12) + n21*log(p21) + n22*log(1-p11-p12-p21))
}

# check answer with MLEs
logLH(c(28/(684+676), 656/(684+676),18/(684+676)))

# check MLEs are what we expect
x <- list(label=c("p11","p12","p21"),est=c(28/(684+676),656/(684+676),18/(684+676)),
low=c(0.001,0.001,0.001),upp=c(1,1,1))
dfp(x,f=logLH)

# Now get CI on theta
x <- list(label=c("p1p","pp1","theta"),est=c((28+656)/(684+676),(28+18)/(684+676),1.56),low=c(0,0,0),upp=c(1,1,100))
plkhci(x,nlogf,"theta")

## 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
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


# Exact Small-sample confidence intervals for 2x2 tables

library(ctest)
(fisher.test(matrix(c(3,1,1,3),byrow=T,ncol=2), alternative="two.sided", or=1))


        Fisher's Exact Test for Count Data

data:  matrix(c(3, 1, 1, 3), byrow = T, ncol = 2) 
p-value = 0.4857
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
   0.2117329 621.9337505 
sample estimates:
odds ratio 
  6.408309 

# Using non-central hypergeometric function

f<-function(x, alpha,t0){
    resl<-hypergeometric(4,4,8,x[1])
    resu<-hypergeometric(4,4,8,x[2])
    sum(c(1-resl$p(t0-1) - alpha/2, resu$p(t0)-alpha/2)^2)
}

optim(par=c(.22, 622), fn=f, method="BFGS", alpha=.05, t0=3, control=list(parscale=c(1,100)))

$par
[1]   0.2117342 626.2385687

$value
[1] 1.345465e-13

$counts
function gradient 
      91       80 

$convergence
[1] 0

$message
NULL

res<-hypergeometric(4,4,8,.212)     # lower endpoint
1-res$p(2)
[1] 0.02505818

res<-hypergeometric(4,4,8,626.24)        # Upper endpoint
res$p(3)    # p is the CDF returned by hypergeometric; we evaluate it at t=3
[1] 0.0250017



# mid-p-value


f<-function(x, alpha,t0){
    resl<-hypergeometric(4,4,8,x[1])
    resu<-hypergeometric(4,4,8,x[2])
    sum(c(1-resl$p(t0-1) -.5*resl$d(t0) - alpha/2, resu$p(t0-1) + .5*resu$d(t0) -alpha/2)^2)
}

optim(par=c(.22, 622), fn=f, method="BFGS", alpha=.05, t0=3, control=list(parscale=c(1,100)))

$par
[1]   0.3100547 308.5567363

res<-hypergeometric(4,4,8,308.5567363)
res$p(2) + .5*res$d(3)
[1] 0.02500000

res<-hypergeometric(4,4,8,0.3100547)
1-res$p(2) - .5*res$d(3)
[1] 0.02499998