##########################################
### Example: Friench skiers
### Test of independence
### Likelihood Ratio Statistics 
### Fisher's Exact Test
##########################################

### Here is one way to read the data vector of values with labels for the table
ski<-matrix(c(31, 17, 109, 122), ncol=2, dimnames=list(Treatment=c("Placebo", "VitaminC"), Cold=c("Cold", "NoCold")))
ski

### Pearson's Chi-squared test  with Yates' continuity correction
result<-chisq.test(ski)
result

###Let's look at the obseved, expected values and the residuals
result$observed
result$expected
result$residuals

###Pearson's Chi-squared test  withOUT Yates' continuity correction
result<-chisq.test(ski, correct=FALSE)
result
result$observed
result$expected
result$residuals

### Let us look at the Percentage, Row Percentage and Column Percentage 
### of the total observations contained in each cell.

Contingency_Table=list(Frequency=ski,Expected=result$expected,Percentage=prop.table(ski),RowPercentage=prop.table(ski,1),ColPercentage=prop.table(ski,2))
Contingency_Table

Percentage=100*ski/sum(ski)
RowSums=rowSums(ski)
RowPercentage=100*rbind(ski[1,]/RowSums[1],ski[2,]/RowSums[2])
ColSums=colSums(ski)
ColPercentage=100*cbind(ski[,1]/ColSums[1],ski[,2]/ColSums[2])
Percentage
RowPercentage
ColPercentage


### Pearson's Chi-squared test  with Yates' continuity correction
result=chisq.test(ski)
result
result$observed
result$expected
result$residuals

### Likelihood Ratio Chi-Squared Statistic
G2=2*sum(ski*log(ski/result$expected))
G2
pvalue=1-pchisq(2*sum(ski*log(ski/result$expected)),df=1)
pvalue

### OR USE OUR function LRstats()
### You first must compile (run) this function
LRstats(ski)

### Fisher's Exact Test
Fisher_Exact_TwoSided=fisher.test(ski,alternative = "two.sided")
Fisher_Exact_Less=fisher.test(ski,alternative = "less")
Fisher_Exact_Greater=fisher.test(ski,alternative = "greater")
rbind(Fisher_Exact_TwoSided,Fisher_Exact_Less,Fisher_Exact_Greater)

### Column 1 Risk Estmates

risk1_col1=ski[1,1]/RowSums[1]
risk2_col1=ski[2,1]/RowSums[2]
rho1=risk1_col1/risk2_col1
total1=ColSums[1]/sum(RowSums)
diff1=risk2_col1-risk1_col1
rbind(risk1_col1,risk2_col1,total1,diff1)

### The confidence interval for the difference in proportions for column 1

SE_diff1=sqrt(risk1_col1*(1-risk1_col1)/RowSums[1]+risk2_col1*(1-risk2_col1)/RowSums[2])
CI_diff1=cbind(diff1-qnorm(0.975)*SE_diff1,diff1+qnorm(0.975)*SE_diff1)
SE_diff1
CI_diff1

### Column 2 Risk Estmates

risk1_col2=ski[1,2]/RowSums[1]
risk2_col2=ski[2,2]/RowSums[2]
total2=ColSums[2]/sum(RowSums)
diff2=risk2_col2-risk1_col2
rbind(risk1_col2,risk2_col2,total2,diff2)

### The confidence interval for the difference in proportions for column 2

SE_diff2=sqrt(risk1_col2*(1-risk1_col2)/RowSums[1]+risk2_col2*(1-risk2_col2)/RowSums[2])
CI_diff2=cbind(diff2-qnorm(0.975)*SE_diff2,diff2+qnorm(0.975)*SE_diff2)
SE_diff2
CI_diff2

### Estimate of the Odds of the two rows

odds1=(ski[2,1]/RowSums[2])/(ski[1,1]/RowSums[1])
odds2=(ski[2,2]/RowSums[2])/(ski[1,2]/RowSums[1])

### Odds Ratio

oddsratio=odds1/odds2
odds1
odds2
oddsratio

### Confidence Interval of the odds ratio
log_CI=cbind(log(oddsratio)-qnorm(0.975)*sqrt(sum(1/ski)),log(oddsratio)+qnorm(0.975)*sqrt(sum(1/ski)))
CI_oddsratio=exp(log_CI)
CI_oddsratio

#################################
###using the 'vcd' package
install.packages("vcd")
library(vcd)
## To get the deviance statistics, pearson X^2, and a few others
assocstats(ski)

oddsratio(ski, log=FALSE)
lor=oddsratio(ski) ## OR on the log scale
lor
confint(lor) ## CI on the log scale
exp(confint(lor)) ## CI on the basic scale