R vs SAS CMH

Cochran-Mantel-Haenszel Test

The CMH procedure tests for conditional independence in partial contingency tables for a 2 x 2 x K design. However, it can be generalized to tables of X x Y x K dimensions.

Naming Convention

For the remainder of this document, we adopt the following naming convention when referring to variables of a contingency table:

X = exposure
Y = response
K = control

Scale

The scale of the exposure (X) and response (Y) variables dictate which test statistic is computed for the contingency table. Each test statistic is evaluated on different degrees of freedom (df):

General association statistic (X and Y both nominal) results in (X-1) * (Y-1) dfs
Row mean scores statistic (X is nominal and Y is ordinal) results in X-1 dfs
Nonzero correlation statistic (X and Y both ordinal) results in 1 df

Testing Strategy

Data

To begin investigating the differences in the SAS and R implementations of the CMH test, we decided to use the CDISC Pilot data set, which is publicly available on the PHUSE Test Data Factory repository. We applied very basic filtering conditions upfront (see below) and this data set served as the basis of the examples to follow.

# A tibble: 6 × 36
  STUDYID  SITEID SITEGR1 USUBJID TRTSDT     TRTEDT     TRTP  TRTPN   AGE AGEGR1
  <chr>    <chr>  <chr>   <chr>   <date>     <date>     <chr> <dbl> <dbl> <chr> 
1 CDISCPI… 701    701     01-701… 2014-01-02 2014-07-02 Plac…     0    63 <65   
2 CDISCPI… 701    701     01-701… 2012-08-05 2012-09-01 Plac…     0    64 <65   
3 CDISCPI… 701    701     01-701… 2013-07-19 2014-01-14 Xano…    81    71 65-80 
4 CDISCPI… 701    701     01-701… 2014-03-18 2014-03-31 Xano…    54    74 65-80 
5 CDISCPI… 701    701     01-701… 2014-07-01 2014-12-30 Xano…    81    77 65-80 
6 CDISCPI… 701    701     01-701… 2013-02-12 2013-03-09 Plac…     0    85 >80   
# ℹ 26 more variables: AGEGR1N <dbl>, RACE <chr>, RACEN <dbl>, SEX <chr>,
#   ITTFL <chr>, EFFFL <chr>, COMP24FL <chr>, AVISIT <chr>, AVISITN <dbl>,
#   VISIT <chr>, VISITNUM <dbl>, ADY <dbl>, ADT <date>, PARAMCD <chr>,
#   PARAM <chr>, PARAMN <dbl>, AVAL <dbl>, ANL01FL <chr>, DTYPE <chr>,
#   AWRANGE <chr>, AWTARGET <dbl>, AWTDIFF <dbl>, AWLO <dbl>, AWHI <dbl>,
#   AWU <chr>, QSSEQ <dbl>

Schemes

In order to follow a systematic approach to testing, and to cover variations in the CMH test, we considered the traditional 2 x 2 x K design as well as scenarios where the generalized CMH test is employed (e.g. 5 x 3 x 3).

We present 5 archetype test scenarios that illustrate diverging results, possibly related to sparse data and possibly considered edge cases.

Number	Schema	Variables	Relevant Test	Description
1	2x2x2	X = TRTP, Y = SEX, K = AGEGR1	General Association	TRTP and AGEGR1 were limited to two categories, overall the the groups were rather balanced
3	2x2x3	X = TRTP, Y = SEX, K = RACE	General Association	Gives back NaN in R because RACE is very imbalanced
6	2x5x2	X = TRTP, Y = AVAL, K = SEX	Row Means	Compare Row Means results for R and SAS because Y is ordinal
9	3x5x17	X = TRTP, Y = AVAL, K = SITEID	Row Means	SITEID has many strata and provokes sparse groups, AVAL is ordinal, therefore row means statistic applies here, R threw an error
10	5x3x3	X = AVAL, Y = AGEGR1, K = TRTP	Correlation	X and Y are ordinal variables and therefore the correlation statistics has to be taken here

Results

Here the results can be seen:

CMH Statistics

Loading required package: vcd

Loading required package: grid

Loading required package: gnm


Attaching package: 'vcdExtra'

The following object is masked from 'package:dplyr':

    summarise

scenarios this is a test

As it can be seen, there are two schemata where R does not provide any results:

	Test	Chi-Square		df		p-value
	Test	SAS	R	SAS	R	SAS	R
1	Correlation	0.01767040	0.0008689711	1	1	0.89424870	0.9764831
	Row Means	0.01767040	2.4820278527	1	2	0.89424870	0.2890910
	General Association	0.01767040	2.4820278527	1	2	0.89424870	0.2890910
3^*	Correlation	0.00278713	NaN	1	1	0.95789662	NaN
	Row Means	2.38606985	NaN	2	2	0.30329938	NaN
	General Association	2.38606985	NaN	2	2	0.30329938	NaN
6	Correlation	1.14720000	0.1115439738	1	1	0.28410000	0.7383931
	Row Means	1.14720000	2.6632420358	1	2	0.28410000	0.2640489
	General Association	2.56720000	6.5238474681	4	8	0.63260000	0.5887637
9^†	Correlation	0.08544312	NA	1	NA	0.77005225	NA
	Row Means	2.47631367	NA	2	NA	0.28991809	NA
	General Association	7.03387844	NA	8	NA	0.53298189	NA
10	Correlation	2.73816085	0.8715295423	1	1	0.09797747	0.3505322
	Row Means	4.40701092	3.0445270087	4	4	0.35371641	0.5504018
	General Association	5.73053819	5.7305381934	8	8	0.67738613	0.6773861
^* Reason for NaN in schema 3: Stratum k = AMERICAN INDIAN OR ALASKA NATIVE can not be compared because there are only values for one treatment and one gender.
^† Reason for Error 4: For large sparse table (many strata) CMHTest will occasionally throw an error in solve.default(AVA) because of singularity

Summary and Recommendation

Having explored the available R packages to calculate the CMH statistics, using R can only be recommended if the analysis design is equivalent to 2 x 2 x K. Then, the base mantelhaen.test() function as well as the vcdExtra package show reliable results which are equal to the output of the SAS function. The same is true for the common odds ratio even though there is a marked difference in decimals.

For the generalized version of the cmh test no R package can be recommended so far. SAS and R outputs differ substantially (possibly due to the underlying subroutines or functions) and the vcdExtra package seems to deliver inconsistent results outside the 2 x 2 x K design.

References

Accessible Summary: https://online.stat.psu.edu/stat504/lesson/4/4.4

An Introduction to Categorical Data Analysis 2nd Edition (Agresti): http://users.stat.ufl.edu/~aa/

SAS documentation (Specification): https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_freq_examples07.html

SAS documentation (Theoretical Basis + Formulas): https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/procstat/procstat_freq_details92.html

Original Paper 1: https://doi.org/10.2307%2F3001616

Original Paper 2: https://doi.org/10.1093/jnci/22.4.719