R语言cochran_armitage检验

Cochran Armitage检验是一种线性趋势检验，常用于自变量是有序分类变量，而因变量是二分类变量的资料，可以用来检验自变量和因变量存不存在线性趋势。

注意和Cochran-Mantel-Haenszel检验区分，CMH检验是研究两个分类变量之间关联性的一种检验方法。但有时数据除了我们研究的变量外，还混杂或隐含了其它的变量，如果将这些变量纳入分析中，则有可能得出完全不同的结论，著名的Simpson悖论就是这个问题的典型案例。

换句话说，在2 x 2表格数据的基础上，引入了第三个分类变量，称之为混杂变量。混杂变量的引入使得CMH检验可以用于分析分层样本，作为生物统计学领域的一种常用技术，该检验常用于疾病对照研究。

现在我们想要了解某种药物剂量和疗效之间的关系，药物剂量有50mg，100mg，200mg，300mg，500mg，5个水平，疗效分为有效/无效两个水平。这种情况可以使用Cochran Armitage检验。

df <- matrix(c(13, 136, 17, 125, 16, 104, 32, 149, 9, 45), 
             nrow = 5, byrow = T,
             dimnames = list("Dose" = c("50", "100", "200", "300", "500"),
                             "effect" = c("Yes", "No"))
             )

df
##      effect
## Dose  Yes  No
##   50   13 136
##   100  17 125
##   200  16 104
##   300  32 149
##   500   9  45

首先可以计算一下不同药物剂量下的有效率是多少：

df[,1]/rowSums(df)
##         50        100        200        300        500 
## 0.08724832 0.11971831 0.13333333 0.17679558 0.16666667

可以看到随着药物剂量增加，有效率整体也是增加的，下面使用CAM检验验证一下。

使用DescTools包中的CochranArmitageTest()函数进行检验：

DescTools::CochranArmitageTest(df)
## 
## 	Cochran-Armitage test for trend
## 
## data:  df
## Z = 2.2116, dim = 5, p-value = 0.02699
## alternative hypothesis: two.sided

结果显示P值为p-value = 0.02699，小于0.05，可以认为疗效会随着药物剂量增加而增加。

现在的df是一个频数统计表类型的数据，我们可以把它变成每行一个患者的数据，然后进行logistic回归看看结果。

df1 <- rstatix::counts_to_cases(df)
psych::headTail(df1)
##     Dose effect
## 1     50    Yes
## 2     50    Yes
## 3     50    Yes
## 4     50    Yes
## ... <NA>   <NA>
## 643  500     No
## 644  500     No
## 645  500     No
## 646  500     No

把Dose变成数值型：

df1$Dose <-  as.numeric(factor(df1$Dose))

summary(glm(effect~Dose, data = df1,family = binomial()))
## 
## Call:
## glm(formula = effect ~ Dose, family = binomial(), data = df1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1762   0.4435   0.4912   0.6001   0.6616  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.48493    0.29598   8.396   <2e-16 ***
## Dose        -0.21544    0.08985  -2.398   0.0165 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 510.57  on 645  degrees of freedom
## Residual deviance: 504.71  on 644  degrees of freedom
## AIC: 508.71
## 
## Number of Fisher Scoring iterations: 4

logistic回归的结果也显示，剂量的P值是小于0.05的。

下面是CMH检验的一个补充。

默认的CMH检验只能进行3个变量的检验，vcdExtra中的CMHtest()可以进行两个变量的CMH检验。

vcdExtra::CMHtest(df, types = "cor")
## Cochran-Mantel-Haenszel Statistics for Dose by effect 
## 
##           AltHypothesis  Chisq Df     Prob
## cor Nonzero correlation 5.8217  1 0.015829