> ##########################################
> ### 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
          Cold
Treatment  Cold NoCold
  Placebo    31    109
  VitaminC   17    122
> 
> ### Pearson's Chi-squared test  with Yates' continuity correction
> result<-chisq.test(ski)
> result

	Pearson's Chi-squared test with Yates' continuity correction

data:  ski 
X-squared = 4.1407, df = 1, p-value = 0.04186

> 
> ###Let's look at the obseved, expected values and the residuals
> result$observed
          Cold
Treatment  Cold NoCold
  Placebo    31    109
  VitaminC   17    122
> result$expected
          Cold
Treatment      Cold  NoCold
  Placebo  24.08602 115.914
  VitaminC 23.91398 115.086
> result$residuals
          Cold
Treatment       Cold     NoCold
  Placebo   1.408787 -0.6421849
  VitaminC -1.413846  0.6444908
> 
> ###Pearson's Chi-squared test  withOUT Yates' continuity correction
> result<-chisq.test(ski, correct=FALSE)
> result

	Pearson's Chi-squared test

data:  ski 
X-squared = 4.8114, df = 1, p-value = 0.02827

> result$observed
          Cold
Treatment  Cold NoCold
  Placebo    31    109
  VitaminC   17    122
> result$expected
          Cold
Treatment      Cold  NoCold
  Placebo  24.08602 115.914
  VitaminC 23.91398 115.086
> result$residuals
          Cold
Treatment       Cold     NoCold
  Placebo   1.408787 -0.6421849
  VitaminC -1.413846  0.6444908
> 
> ### 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
$Frequency
          Cold
Treatment  Cold NoCold
  Placebo    31    109
  VitaminC   17    122

$Expected
          Cold
Treatment      Cold  NoCold
  Placebo  24.08602 115.914
  VitaminC 23.91398 115.086

$Percentage
          Cold
Treatment       Cold   NoCold
  Placebo  0.1111111 0.390681
  VitaminC 0.0609319 0.437276

$RowPercentage
          Cold
Treatment       Cold    NoCold
  Placebo  0.2214286 0.7785714
  VitaminC 0.1223022 0.8776978

$ColPercentage
          Cold
Treatment       Cold    NoCold
  Placebo  0.6458333 0.4718615
  VitaminC 0.3541667 0.5281385

> 
> 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
          Cold
Treatment      Cold  NoCold
  Placebo  11.11111 39.0681
  VitaminC  6.09319 43.7276
> RowPercentage
         Cold   NoCold
[1,] 22.14286 77.85714
[2,] 12.23022 87.76978
> ColPercentage
             [,1]     [,2]
Placebo  64.58333 47.18615
VitaminC 35.41667 52.81385
> 
> 
> ### Pearson's Chi-squared test  with Yates' continuity correction
> result=chisq.test(ski)
> result

	Pearson's Chi-squared test with Yates' continuity correction

data:  ski 
X-squared = 4.1407, df = 1, p-value = 0.04186

> result$observed
          Cold
Treatment  Cold NoCold
  Placebo    31    109
  VitaminC   17    122
> result$expected
          Cold
Treatment      Cold  NoCold
  Placebo  24.08602 115.914
  VitaminC 23.91398 115.086
> result$residuals
          Cold
Treatment       Cold     NoCold
  Placebo   1.408787 -0.6421849
  VitaminC -1.413846  0.6444908
> 
> ### Likelihood Ratio Chi-Squared Statistic
> G2=2*sum(ski*log(ski/result$expected))
> G2
[1] 4.871697
> pvalue=1-pchisq(2*sum(ski*log(ski/result$expected)),df=1)
> pvalue
[1] 0.02730064
> 
> ### OR USE OUR function LRstats()
> ### You first must compile (run) this function
> LRstats(ski)
[1] 4.87169682 0.02730064
> 
> ### 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)
                      p.value    conf.int  estimate null.value alternative
Fisher_Exact_TwoSided 0.03849249 Numeric,2 2.035861 1          "two.sided"
Fisher_Exact_Less     0.9910067  Numeric,2 2.035861 1          "less"     
Fisher_Exact_Greater  0.02052272 Numeric,2 2.035861 1          "greater"  
                      method                               data.name
Fisher_Exact_TwoSided "Fisher's Exact Test for Count Data" "ski"    
Fisher_Exact_Less     "Fisher's Exact Test for Count Data" "ski"    
Fisher_Exact_Greater  "Fisher's Exact Test for Count Data" "ski"    
> 
> ### 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)
               Placebo
risk1_col1  0.22142857
risk2_col1  0.12230216
total1      0.17204301
diff1      -0.09912641
> 
> ### 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
   Placebo 
0.04476243 
> CI_diff1
               [,1]        [,2]
VitaminC -0.1868592 -0.01139366
> 
> ### 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)
              Placebo
risk1_col2 0.77857143
risk2_col2 0.87769784
total2     0.82795699
diff2      0.09912641
> 
> ### 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
   Placebo 
0.04476243 
> CI_diff2
               [,1]      [,2]
VitaminC 0.01139366 0.1868592
> 
> ### 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
 VitaminC 
0.5523323 
> odds2
VitaminC 
1.127318 
> oddsratio
 VitaminC 
0.4899524 
> 
> ### 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
              [,1]      [,2]
VitaminC 0.2569419 0.9342709
> 
> #################################
> ###using the 'vcd' package
> install.packages("vcd")
trying URL 'http://lib.stat.cmu.edu/R/CRAN/bin/macosx/leopard/contrib/2.14/vcd_1.2-13.tgz'
Content type 'application/x-gzip' length 1349845 bytes (1.3 Mb)
opened URL
==================================================
downloaded 1.3 Mb


The downloaded packages are in
	/var/folders/rl/g9pzw2wj2ql6mwznzb_z6h680000gn/T//RtmpJkvj0g/downloaded_packages
> library(vcd)
> ## To get the deviance statistics, pearson X^2, and a few others
> assocstats(ski)
                    X^2 df P(> X^2)
Likelihood Ratio 4.8717  1 0.027301
Pearson          4.8114  1 0.028272

Phi-Coefficient   : 0.131 
Contingency Coeff.: 0.13 
Cramer's V        : 0.131 
> 
> oddsratio(ski, log=FALSE)
[1] 2.041015
> lor=oddsratio(ski) ## OR on the log scale
> lor
[1] 0.713447
> confint(lor) ## CI on the log scale
            lwr     upr
[1,] 0.07477366 1.35212
> exp(confint(lor)) ## CI on the basic scale
         lwr      upr
[1,] 1.07764 3.865613