有双向弗里德曼检验吗?

机器算法验证 r 排名 弗里德曼测试
2022-04-12 16:23:18

我对来自同一患者的 3 个不同位置的 5 种不同蛋白质的浓度(从 0 到 5)进行了排名。总共有 15 次测量每位患者和 61 名患者,因此 915 次观察。

我想知道的是:

  • 3 个位置的 5 种蛋白质浓度相同。
  • 任何蛋白质都以较高浓度存在于一个特定位置。

我想我需要一个双向弗里德曼方差分析,因为我有 2 个分类变量(蛋白质和位置)和 1 个序数变量(排名浓度)。我的问题是:

  1. 如何在 R 中运行测试?
  2. 引导是我唯一的解决方案吗?
  3. 互动呢?
  4. 我读过关于比例赔率的文章。能帮上忙吗?

为了更清楚,我这里是部分数据:

  Protein   Location Concentration
    Prot1   Loc1    0
    Prot1   Loc1    0
    Prot1   Loc1    1
    Prot1   Loc1    0
    Prot1   Loc1    0
    Prot1   Loc1    2
    Prot1   Loc1    1
    Prot1   Loc1    1
    Prot1   Loc1    1
    Prot1   Loc2    1
    Prot1   Loc2    0
    Prot1   Loc2    0
    Prot1   Loc2    1
    Prot1   Loc2    1
    Prot1   Loc2    3
    Prot1   Loc2    0
    Prot1   Loc2    0
    Prot1   Loc2    2
    Prot1   Loc2    0
    Prot1   Loc2    0
    Prot1   Loc3    0
    Prot1   Loc3    0
    Prot1   Loc3    0
    Prot1   Loc3    1
    Prot1   Loc3    1
    Prot1   Loc3    2
    Prot1   Loc3    0
    Prot1   Loc3    1
    Prot1   Loc3    1
    Prot1   Loc3    0
    Prot1   Loc3    0
    Prot2   Loc1    1
    Prot2   Loc1    1
    Prot2   Loc1    2
    Prot2   Loc1    0
    Prot2   Loc1    0
    Prot2   Loc1    1
    Prot2   Loc1    0
    Prot2   Loc1    0
    Prot2   Loc1    0
    Prot2   Loc1    2
    Prot2   Loc1    0
    Prot2   Loc2    2
    Prot2   Loc2    2
    Prot2   Loc2    1
    Prot2   Loc2    1
    Prot2   Loc2    0
    Prot2   Loc2    1
    Prot2   Loc2    3
    Prot2   Loc2    0
    Prot2   Loc2    0
    Prot2   Loc2    3
    Prot2   Loc2    0
    Prot2   Loc3    3
    Prot2   Loc3    1
    Prot2   Loc3    2
    Prot2   Loc3    1
    Prot2   Loc3    0
    Prot2   Loc3    1
    Prot2   Loc3    0
    Prot2   Loc3    0
    Prot2   Loc3    0
    Prot2   Loc3    1
    Prot2   Loc3    0
    Prot3   Loc1    1
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc1    1
    Prot3   Loc1    2
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc1    0
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc2    1
    Prot3   Loc2    5
    Prot3   Loc2    1
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc2    0
    Prot3   Loc3    1
    Prot3   Loc3    0
    Prot3   Loc3    0
    Prot3   Loc3    0
    Prot3   Loc3    0
    Prot3   Loc3    5
    Prot3   Loc3    3
    Prot3   Loc3    2
    Prot3   Loc3    0
    Prot3   Loc3    0
    Prot3   Loc3    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc1    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc2    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot4   Loc3    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc1    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc2    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
    Prot5   Loc3    0
3个回答

您的数据是有序评级,因此您需要某种形式的有序逻辑回归。但我也收集到您的数据不是独立的(“……每位患者 15 次测量……”),因此也需要考虑到这一点。因此,这里的适当方法是混合效应序数逻辑回归。在 R 中,混合效应 OLR 模型可以与序数包相匹配。

这是您的数据的简短演示:

prtdf = read.table(text="Protein   Location Concentration
                         Prot1     Loc1     0
                         ...
                         Prot5     Loc3     0", header=T)

这些数据有几个问题。首先,它们不是很平衡(这实际上没什么大不了的):

with(prtdf, table(Protein, Location))
#        Location
# Protein Loc1 Loc2 Loc3
#   Prot1    9   11   11
#   Prot2   11   11   11
#   Prot3   11   11   11
#   Prot4   11   11   11
#   Prot5   11   11   11

至关重要的是,他们缺少患者 ID 指示器。我会编一个,假设每个类别中的顺序是一致的,并且根据患者 ID(这在现实中很可能是完全错误的,所以要预先警告):

prtdf$ID = c(1:9, rep(1:11, times=14))
prtdf    = prtdf[,c(4,1:3)]
head(prtdf, 10)
#    ID Protein Location Concentration
# 1   1   Prot1     Loc1             0
# 2   2   Prot1     Loc1             0
# 3   3   Prot1     Loc1             1
# 4   4   Prot1     Loc1             0
# 5   5   Prot1     Loc1             0
# 6   6   Prot1     Loc1             2
# 7   7   Prot1     Loc1             1
# 8   8   Prot1     Loc1             1
# 9   9   Prot1     Loc1             1
# 10  1   Prot1     Loc2             1
tail(prtdf, 12)
#     ID Protein Location Concentration
# 152 11   Prot5     Loc2             0
# 153  1   Prot5     Loc3             0
# 154  2   Prot5     Loc3             0
# 155  3   Prot5     Loc3             0
# 156  4   Prot5     Loc3             0
# 157  5   Prot5     Loc3             0
# 158  6   Prot5     Loc3             0
# 159  7   Prot5     Loc3             0
# 160  8   Prot5     Loc3             0
# 161  9   Prot5     Loc3             0
# 162 10   Prot5     Loc3             0
# 163 11   Prot5     Loc3             0

接下来,我们需要确保ID并被Concentration适当地归类为因素。(另请注意,您缺少任何4's in Concentration。)

with(prtdf, table(Concentration))
# Concentration
#   0   1   2   3   4   5 
# 120  26  10   5   0   2 
prtdf$Concentration = factor(prtdf$Concentration, levels=0:5, ordered=T)
prtdf$ID            = factor(prtdf$ID,            levels=1:11)

现在我们可以尝试拟合一个模型:

library(ordinal)
mod = clmm(Concentration~Protein*Location+(1|ID), data=prtdf, Hess=T, nAGQ=17)
# Warning message:
#   (1) Hessian is numerically singular: parameters are not uniquely determined 
# In addition: Absolute convergence criterion was met, but relative criterion was not met

那崩溃了。问题似乎是所有的Concentrationsin"Prot4""Prot5"都是 0:

aggregate(Concentration~Protein*Location, data=prtdf, function(x){ mean(as.numeric(x)) })
#    Protein Location Concentration
# 1    Prot1     Loc1      1.666667
# 2    Prot2     Loc1      1.636364
# 3    Prot3     Loc1      1.363636
# 4    Prot4     Loc1      1.000000
# 5    Prot5     Loc1      1.000000
# 6    Prot1     Loc2      1.727273
# 7    Prot2     Loc2      2.181818
# 8    Prot3     Loc2      1.636364
# 9    Prot4     Loc2      1.000000
# 10   Prot5     Loc2      1.000000
# 11   Prot1     Loc3      1.545455
# 12   Prot2     Loc3      1.818182
# 13   Prot3     Loc3      2.000000
# 14   Prot4     Loc3      1.000000
# 15   Prot5     Loc3      1.000000
table(prtdf[prtdf$Protein%in%c("Prot4","Prot5"), "Concentration"])
#  0  1  2  3  4  5 
# 66  0  0  0  0  0 

我们将从分析中简单地排除这些级别:

mod2 = clmm(Concentration~Protein*Location+(1|ID), data=prtdf, Hess=T, nAGQ=7,
            subset=!prtdf$Protein%in%c("Prot4","Prot5"))
summary(mod2)
# Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
# quadrature approximation with 7 quadrature points 
# 
# formula: Concentration ~ Protein * Location + (1 | ID)
# data:    prtdf
# subset:  !prtdf$Protein %in% c("Prot4", "Prot5")
# 
#  link  threshold nobs logLik  AIC    niter     max.grad cond.H
#  logit flexible  97   -106.48 238.97 760(2283) 1.06e-04 2.0e+02
# 
# Random effects:
# Groups Name        Variance Std.Dev.
# ID     (Intercept) 0.5756   0.7587  
# Number of groups:  ID 11 
# 
# Coefficients:
#                           Estimate Std. Error z value Pr(>|z|)
# ProteinProt2              -0.14342    0.85466  -0.168    0.867
# ProteinProt3              -0.96646    0.92526  -1.045    0.296
# LocationLoc2               0.03622    0.86002   0.042    0.966
# LocationLoc3              -0.17634    0.84245  -0.209    0.834
# ProteinProt2:LocationLoc2  0.96548    1.19393   0.809    0.419
# ProteinProt3:LocationLoc2  0.02491    1.31402   0.019    0.985
# ProteinProt2:LocationLoc3  0.57413    1.18357   0.485    0.628
# ProteinProt3:LocationLoc3  1.13724    1.27533   0.892    0.373
# 
# Threshold coefficients:
#     Estimate Std. Error z value
# 0|1   0.1591     0.6595   0.241
# 1|2   1.6771     0.6904   2.429
# 2|3   2.7789     0.7615   3.649
# 3|5   4.1442     0.9793   4.232

现在这确实返回了一个结果,但是由于您的变量是因子(或多级分类变量),因此各个级别的 p 值并不重要。您想了解整个因素的重要性。特别是,我认为您可能有兴趣了解交互是否重要。我们可以通过拟合加法模型(即没有交互项)并执行嵌套模型测试来测试:

mod2a = clmm(Concentration~Protein+Location+(1|ID), data=prtdf, Hess=T, nAGQ=7,
             subset=!prtdf$Protein%in%c("Prot4","Prot5"))
anova(mod2a, mod2)
# Likelihood ratio tests of cumulative link models:
#   
#       formula:                                      link: threshold:
# mod2a Concentration ~ Protein + Location + (1 | ID) logit flexible
# mod2  Concentration ~ Protein * Location + (1 | ID) logit flexible
# 
#       no.par    AIC  logLik LR.stat df Pr(>Chisq)
# mod2a      9 233.02 -107.51                      
# mod2      13 238.97 -106.48   2.053  4      0.726

对于这些数据,交互作用似乎并不显着。如果您还想测试加法模型中的变量,可以这样方便地完成:

drop1(mod2a, test="Chisq")
# Single term deletions
# 
# Model:
# Concentration ~ Protein + Location + (1 | ID)
#          Df    AIC    LRT Pr(>Chi)
# <none>      233.02                
# Protein   2 232.44 3.4238   0.1805
# Location  2 229.74 0.7212   0.6973

它们在这些数据中也不显着。

为您的问题提供明确的答案:尽管弗里德曼检验只是单向检验,但序数逻辑回归是 Kruskal-Wallis 检验的推广,混合效应 OLR 是 OLR 和弗里德曼检验的推广。自举在这里不太可能对您有所帮助。序数逻辑回归通常称为比例优势模型。

您可以使用具有相互作用的简单线性回归来确定蛋白质、位置(及其相互作用)与浓度的关系。以下是您提供的数据的输出:

> summary(lm(Concentration~Protein*Location, data=prtdf))

Call:
lm(formula = Concentration ~ Protein * Location, data = prtdf)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1818 -0.5455  0.0000  0.0000  4.3636 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)  
(Intercept)                0.66667    0.27918   2.388   0.0182 *
ProteinProt2              -0.03030    0.37645  -0.080   0.9360  
ProteinProt3              -0.30303    0.37645  -0.805   0.4221  
ProteinProt4              -0.66667    0.37645  -1.771   0.0786 .
ProteinProt5              -0.66667    0.37645  -1.771   0.0786 .
LocationLoc2               0.06061    0.37645   0.161   0.8723  
LocationLoc3              -0.12121    0.37645  -0.322   0.7479  
ProteinProt2:LocationLoc2  0.48485    0.51890   0.934   0.3516  
ProteinProt3:LocationLoc2  0.21212    0.51890   0.409   0.6833  
ProteinProt4:LocationLoc2 -0.06061    0.51890  -0.117   0.9072  
ProteinProt5:LocationLoc2 -0.06061    0.51890  -0.117   0.9072  
ProteinProt2:LocationLoc3  0.30303    0.51890   0.584   0.5601  
ProteinProt3:LocationLoc3  0.75758    0.51890   1.460   0.1464  
ProteinProt4:LocationLoc3  0.12121    0.51890   0.234   0.8156  
ProteinProt5:LocationLoc3  0.12121    0.51890   0.234   0.8156  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8375 on 148 degrees of freedom
Multiple R-squared:  0.2019,    Adjusted R-squared:  0.1263 
F-statistic: 2.673 on 14 and 148 DF,  p-value: 0.001637

它表明蛋白质 4 和 5 的浓度略高,但它们在蛋白质浓度的任何位置之间没有差异。此外,位置和蛋白质之间没有相互作用,因此没有一种蛋白质优先集中在任何一个位置。

在阅读了更多之后,我想我得到了自己问题的答案。

由于响应变量是有序的,而其他因素是分类的,我不应该使用简单的线性回归,所以我使用了逻辑回归模型。

R中有两个包可以处理序数变量:lrm {rms}polr {MASS}.

我之所以选择是lrm因为它显示了每个估计的 p 值,但两个结果都是相同的。

遵循@rnso的符号(谢谢!):

> library (rms)

> lrm (Concentration ~ Protein * Location, data = prtdf)

Logistic Regression Model

lrm(formula = Concentration ~ Protein * Location, data = prtdf)

Frequencies of Responses

  0   1   2   3   5 
120  26  10   5   2 

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs           163    LR chi2      60.55    R2       0.380    C       0.805    
max |deriv| 0.002    d.f.            14    g        5.173    Dxy     0.610    
                     Pr(> chi2) <0.0001    gr     176.462    gamma   0.644    
                                           gp       0.111    tau-a   0.263    
                                           Brier    0.083                     

                              Coef     S.E.    Wald Z Pr(>|Z|)
y>=1                           -0.0288  0.5961 -0.05  0.9615  
y>=2                           -1.4135  0.6168 -2.29  0.0219  
y>=3                           -2.4530  0.6863 -3.57  0.0004  
y>=5                           -3.7741  0.9132 -4.13  <0.0001 
Protein=Prot2                  -0.1834  0.8232 -0.22  0.8237  
Protein=Prot3                  -0.9596  0.8933 -1.07  0.2827  
Protein=Prot4                 -10.4962 58.1844 -0.18  0.8568  
Protein=Prot5                 -10.4962 58.1844 -0.18  0.8568  
Location=Loc2                  -0.1326  0.8291 -0.16  0.8729  
Location=Loc3                  -0.3028  0.8180 -0.37  0.7113  
Protein=Prot2 * Location=Loc2   1.0297  1.1526  0.89  0.3717  
Protein=Prot3 * Location=Loc2   0.1827  1.2607  0.14  0.8848  
Protein=Prot4 * Location=Loc2   0.1326 82.2851  0.00  0.9987  
Protein=Prot5 * Location=Loc2   0.1326 82.2851  0.00  0.9987  
Protein=Prot2 * Location=Loc3   0.6113  1.1413  0.54  0.5923  
Protein=Prot3 * Location=Loc3   1.0436  1.2382  0.84  0.3993  
Protein=Prot4 * Location=Loc3   0.3028 82.2850  0.00  0.9971  
Protein=Prot5 * Location=Loc3   0.3028 82.2850  0.00  0.9971  

因此,对于我给出的截断数据集,在蛋白质浓度、位置或它们的组合之间没有发现显着差异。

我使用了加州大学洛杉矶分校数字研究和教育研究所的“序数逻辑回归”指南。

还有什么建议吗?