线性模型残差相关性测试的自动化

Question

是否有方法/代码可用于自动化我正在使用的相关性测试过程？自动化代码应该能够比较数据集中的所有列，并生成相关矩阵作为输出。目前，我正在手动更改代码中的列名称并将结果写入 Excel 工作表。然而，由于有超过 60 个观察柱，这很费力，所以就有了这个问题。

可以在此处找到这种相关性分析方法的详细原因。

我使用 IRIS 数据集准备了一个代码示例。 IRIS 数据集最初不包含重复，但是每个观测和物种可提供 50 个重复。该代码与我用于分析数据集的代码相同。

mm1 <- lm(Sepal.Length ~ factor(Replication) + Species, data=iris)
res1 <- mm1$residuals
mm2 <- lm(Sepal.Width ~ factor(Replication) + Species, data=iris)
res2 <- mm2$residuals
cor.test(res1,res2)

Answer 1

我忽略你的

Replication

，因为你没有提供包含它的数据。不过，以下应该很容易调整。

DF <- iris

library(data.table)
setDT(DF)
DFlong <- melt(DF, id.vars = "Species")

library(nlme)
mm <- lmList(value ~ Species | variable, data = DFlong, na.action = na.fail)

res <- matrix(residuals(mm), nrow = nrow(DF), dimnames = list(NULL, rownames(coef(mm))))

library(psych)
cors <- corr.test(res)
cors$r
#             Sepal.Length Sepal.Width Petal.Length Petal.Width
#Sepal.Length    1.0000000   0.5302358    0.7561642   0.3645064
#Sepal.Width     0.5302358   1.0000000    0.3779162   0.4705346
#Petal.Length    0.7561642   0.3779162    1.0000000   0.4844589
#Petal.Width     0.3645064   0.4705346    0.4844589   1.0000000

cors <- corr.test(res)
print(cors, digits = 2, short = FALSE)
# Call:corr.test(x = res)
# Correlation matrix 
#              Sepal.Length Sepal.Width Petal.Length Petal.Width
# Sepal.Length         1.00        0.53         0.76        0.36
# Sepal.Width          0.53        1.00         0.38        0.47
# Petal.Length         0.76        0.38         1.00        0.48
# Petal.Width          0.36        0.47         0.48        1.00
# Sample Size 
# [1] 150
# Probability values (Entries above the diagonal are adjusted for multiple tests.) 
#              Sepal.Length Sepal.Width Petal.Length Petal.Width
# Sepal.Length            0           0            0           0
# Sepal.Width             0           0            0           0
# Petal.Length            0           0            0           0
# Petal.Width             0           0            0           0
# 
#  Confidence intervals based upon normal theory.  To get bootstrapped values, try cor.ci
#             raw.lower raw.r raw.upper raw.p lower.adj upper.adj
# Spl.L-Spl.W      0.40  0.53      0.64     0      0.36      0.67
# Spl.L-Ptl.L      0.68  0.76      0.82     0      0.65      0.84
# Spl.L-Ptl.W      0.22  0.36      0.50     0      0.22      0.50
# Spl.W-Ptl.L      0.23  0.38      0.51     0      0.21      0.52
# Spl.W-Ptl.W      0.34  0.47      0.59     0      0.30      0.61
# Ptl.L-Ptl.W      0.35  0.48      0.60     0      0.31      0.63

cors$p
#             Sepal.Length  Sepal.Width Petal.Length  Petal.Width
#Sepal.Length 0.000000e+00 1.495184e-11 2.860911e-28 4.524678e-06
#Sepal.Width  2.990368e-12 0.000000e+00 3.724687e-06 3.699152e-09
#Petal.Length 4.768185e-29 1.862343e-06 0.000000e+00 1.340992e-09
#Petal.Width  4.524678e-06 1.233051e-09 3.352479e-10 0.000000e+00
cors$p.adj
#[1] 1.495184e-11 2.860911e-28 3.724687e-06 4.524678e-06 3.699152e-09 1.340992e-09

如果进行多次比较，则应使用调整后的 p 值。