> ### Example: Heart Disease Example Lesson 3 ## > ### Simple line by line R code > ### Nice R code that corresponds to SAS code and output > ####################################################### > > ## enter data > heart <-c(12,8,31,41,307,246,439,245) > heart<-matrix(heart,4,2) > heart=t(heart) > > ## run the chi-squared test of independence & save it into a new object > result<-chisq.test(heart) > result Pearson's Chi-squared test data: heart X-squared = 35.0285, df = 3, p-value = 1.202e-07 > > ## Let's look at the obseved, expected values and the residuals > result\$observed [,1] [,2] [,3] [,4] [1,] 12 8 31 41 [2,] 307 246 439 245 > result\$expected [,1] [,2] [,3] [,4] [1,] 22.08277 17.58315 32.53574 19.79834 [2,] 296.91723 236.41685 437.46426 266.20166 > result\$residuals [,1] [,2] [,3] [,4] [1,] -2.1456212 -2.2853872 -0.26923882 4.764917 [2,] 0.5851431 0.6232594 0.07342547 -1.299464 > > ### Likelihood Ratio Test > LR=2*sum(heart*log(heart/result\$expected)) > LR [1] 31.92124 > LRchisq=1-pchisq(LR,df=(4-1)*(2-1)) > LRchisq [1] 5.43736e-07 > > > ##make sure you have function LRstats() > LRstats(heart) [1] 3.192124e+01 5.437360e-07 > > > ## Let's calculate the conditional probabilities > ## the following function gives the desired marginal, in this case, the counts for the serum groups > serum<-margin.table(heart,2) > serum [1] 319 254 470 286 > > ## let's look at the counts for the four groups with CHD > heart[1,] [1] 12 8 31 41 > > ## then counts for the four groups with NOCHD, which is the second column of data in the dataframe we created above > heart[2,] [1] 307 246 439 245 > > ### conditional probabilities are: > heart[2,]/serum [1] 0.9623824 0.9685039 0.9340426 0.8566434 > heart[1,]/serum [1] 0.03761755 0.03149606 0.06595745 0.14335664 > > ######################################## > ### Nice R code that corresponds to SAS code and output > ####################################################### > heart=matrix(c(12,307,8,246,31,439,41,245), ncol=4, dimnames=list(CHD=c("chd", "nochd"), serum=c("0-199", "200-199","220-259","260+"))) > heart serum CHD 0-199 200-199 220-259 260+ chd 12 8 31 41 nochd 307 246 439 245 > count=heart > > ### Chi-Square Independence Test > result=chisq.test(count) > result\$expected serum CHD 0-199 200-199 220-259 260+ chd 22.08277 17.58315 32.53574 19.79834 nochd 296.91723 236.41685 437.46426 266.20166 > > ### Let us look at the Percentage, Row Percentage and Column Percentage > ### of the total observations contained in each cell. > > Contingency_Table=list(Frequency=count,Expected=result\$expected,Deviation=count-result\$expected,Percentage=prop.table(count),RowPercentage=prop.table(count,1),ColPercentage=prop.table(count,2)) > Contingency_Table \$Frequency serum CHD 0-199 200-199 220-259 260+ chd 12 8 31 41 nochd 307 246 439 245 \$Expected serum CHD 0-199 200-199 220-259 260+ chd 22.08277 17.58315 32.53574 19.79834 nochd 296.91723 236.41685 437.46426 266.20166 \$Deviation serum CHD 0-199 200-199 220-259 260+ chd -10.08277 -9.583145 -1.535741 21.20166 nochd 10.08277 9.583145 1.535741 -21.20166 \$Percentage serum CHD 0-199 200-199 220-259 260+ chd 0.009029345 0.006019564 0.02332581 0.03085026 nochd 0.231000752 0.185101580 0.33032355 0.18434913 \$RowPercentage serum CHD 0-199 200-199 220-259 260+ chd 0.1304348 0.08695652 0.3369565 0.4456522 nochd 0.2481811 0.19886823 0.3548909 0.1980598 \$ColPercentage serum CHD 0-199 200-199 220-259 260+ chd 0.03761755 0.03149606 0.06595745 0.1433566 nochd 0.96238245 0.96850394 0.93404255 0.8566434 > > ###### Computing various measures of association > library(vcd) > assocstats(heart) X^2 df P(> X^2) Likelihood Ratio 31.921 3 5.4374e-07 Pearson 35.028 3 1.2015e-07 Phi-Coefficient : 0.162 Contingency Coeff.: 0.16 Cramer's V : 0.162 > > ### For the Pearson correlation coefficent > ### and Mantel-Haenszel, > ### for IxJ tables, you can also use > ### pears.cor() function. > ### Mak sure you run this function first! > ### c(1,2) and c(1,2,3,4), are the vectors of score values > pears.cor(heart, c(1,2),c(1,2,3,4)) [1] -0.1403189 26.1475279 > ### and this should give you, r=-0.14, M2=26.1475 > > ##Gamma > Gamma.f(heart) \$gamma [1] -0.4265663 \$C [1] 24227 \$D [1] 60271