[R] Pearson's Chi-squared Test

[Peter Dalgaard BSA]

Jorge Magalhães <jmagalhaes at oninetspeed.pt> writes:

>  How i can perform a Pearson's Chi-squared Test  in this data set:

[that was garbled. Here's the original]

                      |         Outcome
   ---------+---------+--------------------------+
   Treatment|  Sex    |None    |Some    |Marked  |  Total
   ---------+---------+--------+--------+--------+
   Active   |  Female |      6 |      5 |     16 |     27
            |  Male   |      7 |      2 |      5 |     14
   ---------+---------+--------+--------+--------+
   Placebo  |  Female |     19 |      7 |      6 |     32
            |  Male   |     10 |      0 |      1 |     11
   ---------+---------+--------+--------+--------+
   Total                    42       14       28       84

> 
>  if i do:
>  y<- matrix(c(5,6,16, 19,7,6,7,2,5,10,0,1), nc=3)
> 
>  and
> 
>  chisq.test(y) i found that X-squared=20.3176 but the true value is 15.075
>  (http://www.math.yorku.ca/SCS/Courses/grcat/grc2.html#H2_43:Sample)

No you don't, and I don't see the value 15.075 anywhere on that page.

Please check your input, rather than relying on r-helpers to do so:
y<- matrix(c(6,5,16, 19,7,6,7,2,5,10,0,1), nc=3,byrow=T)
             !!!                                !!!!!!!

will give you the 20.3176, but that test (for independence in the 4x3
table) isn't on the page you cite.
 
With correct entry of your "x" table (see below) you might have
calculated the 2x3 marginal over sex and gotten

> chisq.test(margin.table(x,c(1,3)))

        Pearson's Chi-squared test

data:  margin.table(x, c(1, 3))
X-squared = 13.055, df = 2, p-value = 0.001463

which is indeed in the SAS output.

>  Now i tried:
> 
>  the Cochran-Mantel-Haenszel Chi-Squared Test for Count Data:
> 
>  x <- array(c(6, 19, 7, 10,
>          	5, 7, 2, 0,
>  		16, 6, 5, 1),
>        dim = c(2, 2, 3),
>        dimnames = list(
>            Active = c("Female", "Male"),
>            Placebo = c("Female", "Male"),
>            Outcome.Level = c("None", "Some", "Marked")))
> 
>  mantelhaen.test(x)
> 
>  and now X-squared=2.0863
> 
>  What is wrong?

Printing that table ought to have given you a clue. Here's the right
way:

 x <- array(c(6, 19, 7, 10,
              5,  7, 2, 0,
             16,  6, 5, 1),
        dim = c(2, 2, 3),
        dimnames = list(
           treat=c("Active","Placebo"),
           sex=c("Female", "Male"),
           Outcome.Level = c("None", "Some", "Marked")))

However, for the Mantel-Haenszel test you'd want to stratify by sex,
not outcome, as happens when that is the last dimension of the table,
so use aperm to sort the indices differently:

mantelhaen.test(aperm(x,c(3,1,2))

gives 14.6323, as in SAS.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907