#Heart Condition/ Arterial Complication - Independence Tests

AC<-rep(c("YES","NO"),2)
HC<-rep(c("YES","NO"),each=2)
N<-c(50,110,2,35)
ACHC<-data.frame(factor(AC),factor(HC),N)
names(ACHC)<-c("AC","HC","N")

#Chisquare/LR Test; Manual


nmat<-matrix(N,nrow=2,byrow=T)
row.totals<-c(sum(nmat[1,]),sum(nmat[2,]))
col.totals<-c(sum(nmat[,1]),sum(nmat[,2]))

fit.n11<-row.totals[1]*col.totals[1]/sum(N)
fit.n12<-row.totals[1]*col.totals[2]/sum(N)
fit.n21<-row.totals[2]*col.totals[1]/sum(N)
fit.n22<-row.totals[2]*col.totals[2]/sum(N)

fitted.nmat<-matrix(c(fit.n11,fit.n12,fit.n21,fit.n22),nrow=2,byrow=T)

chi.stat<-sum((nmat-fitted.nmat)^2/fitted.nmat)
lr.stat<-2*sum(nmat*log(nmat/fitted.nmat))

chisq.p<-1-pchisq(chi.stat,1)
lr.p<-1-pchisq(lr.stat,1)

print(c(chisq.p,lr.p))

#-----------------------------------------------------------

#Chisquare test using the special SPLUS function
chisq.test(matrix(N,nrow=2,byrow=T),correct=F)
chisq.test(matrix(N,nrow=2,byrow=T),correct=T)

#-----------------------------------------------------------


#Fisher's Exact test using the special SPLUS function
fisher.test(matrix(N,nrow=2,byrow=T))

#-----------------------------------------------------------

#Chisquare /LR using Monte Carlo exact testing, assuming row and column totals fixed.
#Use the rhyper function to simulate new tables
#under independence null hypothesis

nsim<-5000

m<-160
n<-37
k<-52

Nvec<-rhyper(nsim,m,n,k)

chi.statvec<-rep(0,nsim)
lr.statvec<-rep(0,nsim)

for(isim in 1:nsim){

	newN<-Nvec[isim]
	nmat<-matrix(c(newN,m-newN,k-newN,n-(k-newN)),nrow=2,byrow=T)
	row.totals<-c(sum(nmat[1,]),sum(nmat[2,]))
	col.totals<-c(sum(nmat[,1]),sum(nmat[,2]))

	fit.n11<-row.totals[1]*col.totals[1]/sum(N)
	fit.n12<-row.totals[1]*col.totals[2]/sum(N)
	fit.n21<-row.totals[2]*col.totals[1]/sum(N)
	fit.n22<-row.totals[2]*col.totals[2]/sum(N)

	fitted.nmat<-matrix(c(fit.n11,fit.n12,fit.n21,fit.n22),nrow=2,byrow=T)

	chi.statvec[isim]<-sum((nmat-fitted.nmat)^2/fitted.nmat)
	lr.statvec[isim]<-2*sum(nmat*log(nmat/fitted.nmat))

}

#Possible Cell entries n11
table(Nvec)

#Compute Chisq/ EXACT p by Monte Carlo approximation
chisq.exactp<-length(chi.statvec[chi.statvec >= chi.stat])/nsim

lr.exactp<-length(lr.statvec[chi.statvec >= lr.stat])/nsim

print(c(chisq.exactp,lr.exactp))