Which statistical test should I use for 1 ordinal and two categorical variables?

Question

It would be great if someone could help me out

I would like to know if there is an association between 3 variables, mental health, having a health plan and income level. Below is a contingency table.

> ftable(te)
                incomex   lt10   lt15   lt20   lt25   lt35   lt50   lt75   mt75
goodmh hlthpln1                                                                
FALSE  Yes                5458   5164   4992   5234   5110   5953   5692   7832
       No                 2066   1461   1694   1633   1200    836    441    304
TRUE   Yes               12638  15392  21125  27425  35982  49004  55541 104420
       No                 4307   3961   6197   6536   5729   4911   2944   2437

I was thinking of doing a chi square test, however, I am not sure how this works with an ordinal variable, and more than two variables. I thought a bonferroni correction is needed as well, when you do multiple testing.

I also read about a Cochran–Armitage test for trend.

What test do you think is the best in this situation? And how do you suggest to write this in R?

Answer 1

Which statistical test? ... with your data many tests are possible, it depends on the question(s) you are asking from the data ... and you did not tell us! Nevertheless, with such complex data it is better to think in terms of models, and for that, the first you should do is to transform the data from a table into long form. I did that:

library(tidyverse)
structure(tab)
# A tibble: 32 x 4
   hlthpln1 incomex count goodmh
   <fct>      <dbl> <int> <fct> 
 1 Yes           10  5458 FALSE 
 2 Yes           15  5164 FALSE 
 3 Yes           20  4992 FALSE 
 4 Yes           25  5234 FALSE 
 5 Yes           35  5110 FALSE 
 6 Yes           50  5953 FALSE 
 7 Yes           75  5692 FALSE 
 8 Yes           85  7832 FALSE 
 9 No            10  2066 FALSE 
10 No            15  1461 FALSE 
# … with 22 more rows

Some examples: If you are interested in how the probability of having a health plan depends on the other variables, maybe a logistic regression

  mod.log <- mgcv::gam( hlthpln1  ~ s(incomex, by=goodmh, k=5), weight=count, data=tab, family=binomial)

where the ordinal variable incomex is represented via a spline. Or, if you are interested in modeling incomex as a function of the other factors, maybe an ordinal regression:

mod.polr <- MASS::polr( as.ordered(incomex)  ~  hlthpln1*goodmh, data=tab, Hess=TRUE, weight=count)

There are certainly other possibilities ...

To get the data frame tab, use dget on

structure(list(hlthpln1 = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("No", "Yes"
), class = "factor"), incomex = c(10, 15, 20, 25, 35, 50, 75, 
85, 10, 15, 20, 25, 35, 50, 75, 85, 10, 15, 20, 25, 35, 50, 75, 
85, 10, 15, 20, 25, 35, 50, 75, 85), count = c(5458L, 5164L, 
4992L, 5234L, 5110L, 5953L, 5692L, 7832L, 2066L, 1461L, 1694L, 
1633L, 1200L, 836L, 441L, 304L, 12638L, 15392L, 21125L, 27425L, 
35982L, 49004L, 55541L, 104420L, 4307L, 3961L, 6197L, 6536L, 
5729L, 4911L, 2944L, 2437L), goodmh = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("FALSE", 
"TRUE"), class = "factor")), row.names = c(NA, -32L), class = c("tbl_df", 
"tbl", "data.frame"))