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