# SECTION 4.1


library(BiasedUrn)



# Table 2.8 data:

income<-c("<15000","15000-25000","25000-40000",">40000") 
jobsat<-c("VD","LD","MS","VS") 
data<-c(1,2,1,0,3,3,6,1,10,10,14,9,6,7,12,11) 
(table.2.8<-cbind(expand.grid(income=income,jobsat=jobsat),count=data) )
(temp<-tapply(table.2.8$count,table.2.8[,1:2], sum) )


# A function written by Laura Thompson to compute GAMMA  

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) 
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 * pr, "%", sep = ""), CI = paste(c("[", max(CIa[1], -1), 
", ", min(CIa[2], 1), "]"), collapse = "")) 
} 

Gamma2.f(temp)

# Change table to numeric

inc		<-c(7.5,20,32.5,60)
jbsat		<-1:4
tbl.2.8	<-cbind(expand.grid(income=inc,jobsat=jbsat),count = table.2.8$count )
rep.tbl	<-tbl.2.8[,1:2][rep(1:nrow(table.2.8),table.2.8$count),] 
(a	<- cor.test(rep.tbl[,1],rep.tbl[,2]))

(M.sq	<- (a$est^2)*(nrow(res)-1) )
1-pchisq(3.807461,1) 			# M^2 ~~ chisq, 1 df 

cor.test(rep.tbl[,1],rep.tbl[,2])
cor.test(rep.tbl[,1],rep.tbl[,2],method = "spearman")
cor.test(rep.tbl[,1],rep.tbl[,2],method = "kendall")


# Different numeric values

inc		<-c(7.5,20,32.5,60)
jbsat		<-c(1,2,5,10)
tbl.2.8	<-cbind(expand.grid(income=inc,jobsat=jbsat),count = table.2.8$count )
rep.tbl	<-tbl.2.8[,1:2][rep(1:nrow(table.2.8),table.2.8$count),] 
(a	<- cor.test(rep.tbl[,1],rep.tbl[,2]))

(M.sq	<- (a$est^2)*(nrow(res)-1) )
1-pchisq(M.sq,1) 			# M^2 ~~ chisq, 1 df 

cor.test(rep.tbl[,1],rep.tbl[,2],method = "spearman")
cor.test(rep.tbl[,1],rep.tbl[,2],method = "kendall")




# Table 3.7 data:

Alcohol.consumption<-c("0","<1","1-2","3-5","6<=") 
Malformation<-c("Absent","Present") 
data<-c(17066,14464,788,126,37,48,38,5,1,1) 
(table.3.7<-cbind(expand.grid(Alcohol.consumption=Alcohol.consumption,Malformation=Malformation),count=data) )
(temp<-tapply(table.3.7$count,table.3.7[,2:1], sum) )

chisq.test(temp) 


# score configuration  1:

levels(table.3.7$Alcohol.consumption)<-c(0,0.5,1.5,4,7) 
levels(table.3.7$Malformation)<-0:1 

res<-table.3.7[,1:2][rep(1:nrow(table.3.7),table.3.7$count),] 
(m.sq.1 <- (cor(res)[2,1]^2)*(nrow(res)-1) )
1-pchisq(m.sq.1,1)

x <- as.numeric(as.vector(res[,1]))
y <- as.numeric(as.vector(res[,2]))
cor.test(x,y)
summary(lm(y~x))

plot(x,y,xlab="Alcohol consumption score",ylab="Malformation") 
abline(lm(y~x),col=2)

(best.reg <- sapply(split(y,x),mean))
points(unique(x),best.reg,col=3)

plot(x,y,xlab="Alcohol consumption score",ylab="Malformation",ylim=c(0,0.027)) 
abline(lm(y~x),col=2)
points(unique(x),best.reg,col=3)


# score configuration  2:

levels(table.3.7$Alcohol.consumption)<-1:5 
levels(table.3.7$Malformation)<-0:1 

res<-table.3.7[,1:2][rep(1:nrow(table.3.7),table.3.7$count),] 
(m.sq.2 <- (cor(res)[2,1]^2)*(nrow(res)-1) )
1-pchisq(m.sq.2,1)

x <- as.numeric(as.vector(res[,1]))
plot(x,y,xlab="Alcohol consumption score",ylab="Malformation",ylim=c(0,0.027)) 
abline(lm(y~x),col=2)
points(unique(x),best.reg,col=3)



# score configuration  3:

( levels(table.3.7$Alcohol.consumption)<-unique(rank(res[,1])))
 levels(table.3.7$Malformation)<-0:1 

res<-table.3.7[,1:2][rep(1:nrow(table.3.7),table.3.7$count),] 
(m.sq.3 <- (cor(res)[2,1]^2)*(nrow(res)-1) )
1-pchisq(m.sq.3,1)

x <- as.numeric(as.vector(res[,1]))
plot(x,y,xlab="Alcohol consumption score",ylab="Malformation",ylim=c(0,0.027)) 
abline(lm(y~x),col=2)
points(unique(x),best.reg,col=3)






#=====================================================

# SECTION 4.2:  exact tests

dhyper(0:4,4,4,4)

dhyper(3,4,4,4) + dhyper(4,4,4,4)

0.5*dhyper(3,4,4,4) + dhyper(4,4,4,4)

(fisher.test(matrix(c(3,1,1,3),byrow=T,ncol=2), alternative="greater")) 


#	Haashara:   mid-p simulation

#	1. Let's draw cdf of hyper-geometric p-value 

(x.vec		<- 4:0)
(p.value		<- cumsum(dhyper(x.vec,4,4,4)))	#	the values that the p-value statistic assumes
(p.value.pr		<- dhyper(x.vec,4,4,4))			# 	probability statistic assumes each value

#	Q: what is CDF_p (1/2) ?   A: Prob(p-value <= 0.5) = ?  

p.seq		<- seq(0,1,length = 10000)
p.value.CDF	<- rep(NA,10000)
for(i in 1:10000)
		p.value.CDF[i]	<- sum(p.value.pr[p.value <= p.seq[i]])

plot(p.seq,p.value.CDF,type="l",col=3,lwd=2,xlab = "p-value value",ylab = "p-value CDF")
points(p.value,p.value,pch = 3,cex=2)
abline(0,1,col=2)


# Excercise: how do we compute expectation of p-value  

#	2. Now let's add mid-p CDF

(mid.p.value		<- cumsum(dhyper(x.vec,4,4,4)) - dhyper(x.vec,4,4,4)/2) #	the values that the p-value statistic assumes
(mid.p.value.pr			<- dhyper(x.vec,4,4,4))			# 	probability statistic assumes each value

mid.p.value.CDF	<- rep(NA,10000)
for(i in 1:10000)
		mid.p.value.CDF[i]	<- sum(p.value.pr[mid.p.value <= p.seq[i]])

plot(p.seq,p.value.CDF,type="l",col=3,lwd=2,xlab = "p-value value",ylab = "p-value CDF")
abline(0,1,col=2)
lines(p.seq,mid.p.value.CDF,type="l",col=4,lwd=2)



# What is expectation of mid-p-value  

t.test(sample(p.value,size = 10^6,replace = T,prob = p.value.pr))
t.test(sample(mid.p.value,size = 10^6,replace = T,prob = p.value.pr))



#	Now we return to multi-hypergeometrics.......


(n.ij  <- cbind(c(1,6,1),c(2,3,5),c(2,2,2)))
(n     <- sum(n.ij))
(n.pj  <- apply(n.ij,2,sum))
(n.ip  <- apply(n.ij,1,sum))

prod(factorial(n.pj)) *  prod(factorial(n.ip)) / ( factorial(n) * prod(factorial(n.ij)) ) 

dhyper(1,5,19,8) *  dhyper(6,11,8,7)  * dhyper(2,4,12,10) * dhyper(3,5,7,8)
dhyper(2,5,19,6) *  dhyper(2,3,15,10) * dhyper(6,11,8,7)  * dhyper(3,5,7,8)





#	EXACT CI for OR in 2 by 2 tables


(tbl <- rbind(c(25,25),c(10,40)))

fisher.test(tbl) 

(25*40) / (25*10)

#	P-value:

plot(0:35,dhyper(0:35,50,50,35))
2*sum(dhyper(25:35,50,50,35))

#	For CI we need to consider non-central distribution


? dFNCHypergeo
dFNCHypergeo(24,50,50,35,1)
dhyper(24,50,50,35)

theta.ser	<- seq(0,20,length = 10000)
LE.pvec		<- rep(NA,10000)
GE.pvec		<- rep(NA,10000)
for(i in 1:10000)
{
	LE.pvec[i]	<- sum(dFNCHypergeo(25:35,50,50,35,theta.ser[i]))
	GE.pvec[i]	<- sum(dFNCHypergeo(0:25,50,50,35,theta.ser[i]))
}

plot(theta.ser,LE.pvec,type="l",col=2,lwd=2)
lines(theta.ser,GE.pvec,type="l",col=3,lwd=2)

abline(h = 0.05,lty=2)
abline(h = 0.025,lty=2)

max(theta.ser[LE.pvec <= 0.025])
min(theta.ser[GE.pvec <= 0.025])


fisher.test(tbl,or=1.52,alternative="greater") 
fisher.test(tbl,or=1.53,alternative="greater") 
fisher.test(tbl,or=10.87,alternative="less") 
fisher.test(tbl,or=10.86,alternative="less")