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