R：如何修复 Bonferroni 校正的错误编码

Question

我被困在如何在 RStudio 中为 Bonferroni 校正和 Pearson 相关矩阵的原始 P 值进行编码。我是一名学生，并且是 R 的新手。我也对如何获取数据的均值、SD 和 n 表感到迷茫。当我计算 Pearson 相关矩阵时，我只得到了 r 值，而不是原始概率值。我不确定如何编写代码以在 RStudio 中获得它。然后我尝试计算 Bonferroni 校正并收到一条错误消息，指出无法强制列表对象键入 double。如何修复我的代码，这样它就消失了？我还尝试为数据创建一个包含平均值、SD 和 n 的表，但我陷入了如何继续的困境。

我的数据如下：

Tree Height DBA Leaf Diameter
45.3    14.9    0.76
75.2    26.6    1.06
70.1    22.9    1.19
95      31.8    1.59
107.8   35.5    0.43
93      26.2    1.49
91.5    29      1.19
78.5    29.2    1.1
85.2    30.3    1.24
50      16.8    0.67
47.1    12.8    0.98
73.2    28.4    1.2

我已经安装了 dplyr、tidyr、multcomp、multcompview 的软件包 我从 excel CSV（逗号分隔）文件中读取数据，这将创建数据>dataHW8_1 12obs。 3 个变量

summary(dataHW8_1)

然后我创建了数据的散点图 plot(dataHW8_1$Tree_Height,dataHW8_1$DBA,main="散点图树高与胸高 (DBA) 的树干直径",xlab="树高 (cm)",ylab="DBA (cm)") plot(dataHW8_1$Tree_Height,dataHW8_1$Leaf_Diameter,main="散点图树高与叶径",xlab="树高(cm)",ylab="叶径(cm)") plot(dataHW8_1$DBA,dataHW8_1$Leaf_Diameter,main="散点图树干直径在胸高 (DBA) 与叶直径",xlab="DBA (cm)",ylab="叶直径 (cm)")

然后我注意到数据不是线性的，所以我使用 log() 函数对其进行了转换 dataHW8_1log = log(dataHW8_1)

然后我使用转换后的数据重新创建了我的散点图

plot(dataHW8_1log$Tree_Height,dataHW8_1log$DBA,main="Scatterplot of 
Transformed (log)Tree Height Vs Trunk Diameter at Breast Height 
(DBA)",xlab="Tree Height (cm)",ylab="DBA (cm)")
plot(dataHW8_1log$Tree_Height,dataHW8_1log$Leaf_Diameter,main="Scatterplot 
of Transformed (log)Tree Height Vs Leaf Diameter",xlab="Tree Height 
(cm)",ylab="Leaf Diameter (cm)")
plot(dataHW8_1log$DBA,dataHW8_1log$Leaf_Diameter,main="Scatterplot of 
Transformed (log) Trunk Diameter at Breast Height (DBA) Vs Leaf 
Diameter",xlab="DBA (cm)",ylab="Leaf Diameter (cm)")

然后我创建了散点图的矩阵图

pairs(dataHW8_1log)

然后我使用 Pearson 方法计算了相关系数 这没有给出 P 值的未校正矩阵------你是怎么做到的？

cor(dataHW8_1log,method="pearson")

我被困在如何获得数据的原始概率（未校正的 P 值）矩阵

然后我计算了 Bonferroni 校正-----你是怎么做的？

Data$Bonferroni = 
p.adjust(dataHW8_1log, 
method = "bonferroni")

这样做给了我以下错误：

 Error in p.adjust(dataHW8_1log, method = "bonferroni") : 
 (list) object cannot be coerced to type 'double'

我尝试使用 lapply 进行修复，但这并没有解决我的问题

然后我尝试制作一个均值、SD、n 表，但我只能创建以下代码，并且卡在从那里去哪里------你是怎么做到的？

(,data = dataHW8_1log,
      FUN = function(x) c(Mean = mean(x, na.rm = T),
                          n = length(x),
                          sd = sd(x, na.rm = T))

我在网上尝试了以下示例，但没有一个可以帮助我正确编码 Bonferroni 校正。如果有人可以帮助解释我做错了什么以及如何制作矩阵/表格，我将不胜感激。

Answer 1

这是一个使用 50 行 x 10 列示例数据框的示例。

# 50 rows x 10 columns sample dataframe
df <- as.data.frame(matrix(runif(500), ncol = 10));

我们可以显示成对散点图。

# Pairwise scatterplot
pairs(df);

我们现在可以使用cor.test 来获取单个比较的 p 值。我们使用便利函数cor.test.p 对所有成对比较执行此操作。为了给予应得的功劳，函数cor.test.p取自this SO post，并将dataframe作为参数，同时返回一个未校正的p值矩阵。

# cor.test on dataframes
# From: https://stackoverflow.com/questions/13112238/a-matrix-version-of-cor-test
cor.test.p <- function(x) {
    FUN <- function(x, y) cor.test(x, y)[["p.value"]];
    z <- outer(
      colnames(x),
      colnames(x),
      Vectorize(function(i,j) FUN(x[,i], x[,j])));
    dimnames(z) <- list(colnames(x), colnames(x));
    return(z);
}

# Uncorrected p-values from pairwise correlation tests
pval <- cor.test.p(df);

我们现在通过对每一行（或每一列，因为矩阵是对称的）应用 Bonferroni 校正来校正多重假设检验，我们就完成了。请注意，p.adjust 将 p 值向量作为参数。

# Multiple hypothesis-testing corrected p-values
# Note: pval is a symmetric matrix, so it doesn't matter if we correct 
# by column or by row
padj <- apply(pval, 2, p.adjust, method = "bonferroni");
padj;
#V1 V2        V3 V4        V5        V6 V7 V8        V9 V10
#V1   0  1 1.0000000  1 1.0000000 1.0000000  1  1 1.0000000   1
#V2   1  0 1.0000000  1 1.0000000 1.0000000  1  1 1.0000000   1
#V3   1  1 0.0000000  1 0.9569498 1.0000000  1  1 1.0000000   1
#V4   1  1 1.0000000  0 1.0000000 1.0000000  1  1 1.0000000   1
#V5   1  1 0.9569498  1 0.0000000 1.0000000  1  1 1.0000000   1
#V6   1  1 1.0000000  1 1.0000000 0.0000000  1  1 0.5461443   1
#V7   1  1 1.0000000  1 1.0000000 1.0000000  0  1 1.0000000   1
#V8   1  1 1.0000000  1 1.0000000 1.0000000  1  0 1.0000000   1
#V9   1  1 1.0000000  1 1.0000000 0.5461443  1  1 0.0000000   1
#V10  1  1 1.0000000  1 1.0000000 1.0000000  1  1 1.0000000   0