在相同数据上获得不同的 cor() 输出？

Question

我在计算相关性时遇到了问题。 我有两个数据集： d3：一个有 1259 个观测值和 264 个变量。 d4：一个有 1196 个观察值和 185 个变量 - 这些只是血液测试。 两个数据对参与者具有相同的唯一 ID，因此可以合并它们。

合并时（称为：d）我有 1150 个观察值（因为它们没有验血，并且验血数据有很多控制样本等等。当删除没有信息 x 的参与者时（因为这是纳入标准），我们最终会有 1144 个观测值。

在 d4 中，我有两个样本在所有变量中都存在缺失值。他们包括在 1144 名参与者中。

所以在汇总数据管理之后，我想计算信息 x（我们称之为 Edu）和所有血液样本 (184) 之间的相关性。

我通过 d 的子集创建一个新数据集：

dbio <- d %>%
  select(Edu, BMP_6:CCL16)

我使用这些代码： 前两个函数：

cor.prob <- function (X, dfr = nrow(X) - 2) {
  R <- cor(X, use="complete.obs", method = c("pearson"))
  above <- row(R) < col(R)
  r2 <- R[above]^2
  Fstat <- r2 * dfr/(1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr)
  R[row(R) == col(R)] <- NA
  R
}

flattenSquareMatrix <- function(m) {
  if( (class(m) != "matrix") | (nrow(m) != ncol(m))) stop("Must be a square matrix.") 
  if(!identical(rownames(m), colnames(m))) stop("Row and column names must be equal.")
  ut <- upper.tri(m)
  data.frame(i = rownames(m)[row(m)[ut]],
             j = rownames(m)[col(m)[ut]],
             cor=t(m)[ut],
             p=m[ut])
}

然后我开始进行关联：

kor_masterlist_dbio <- flattenSquareMatrix (cor.prob(dbio))

kor_masterlist_dbio_ordnet <- kor_masterlist_dbio[order(-abs(kor_masterlist_dbio$cor)),]

kor_masterlist_dbio$threshold <- as.factor(kor_masterlist_dbio$p < 0.05)

kor_masterlist_dbio_Edu <- subset(kor_masterlist_dbio_ordnet, i == "Edu" & j != "ID")

kor_masterlist_dbio_Edu$threshold <- as.factor(kor_masterlist_dbio_Edu$p < 0.05)


kor_masterlist_dbio_Edu[kor_masterlist_dbio_Edu$threshold==T,]

kor_masterlist_dbio_Edu

现在的问题： 我在 cor() 中使用“complete.obs”。如果我通过以下代码将他们从数据中取出，则所有血液测试都缺失（NA）的两名参与者：

d <- d[rowSums(is.na(d[,3:6]))!=4,]
And end up with d: n = 1142.

然后我从相关性中得到不同的结果。我最终减少了两次血液测试，因为 p 值高于 0.05。我不明白这两个参与者缺少值时的区别。

当使用“pariwise”ind cor() 时，我得到了与 1142 或 1144 名参与者相同的结果。但是最后 4 次血液测试与我使用“complete.obs”得到的不同

其余的相关系数和 p 值略有不同。但仍然是相同的排名。

我希望有人可以帮助我。无论是否包括参与者，我都应该得到相同的结果。

我尝试制作以下可重现的小示例来显示我的问题。如果你有/没有运行它，你会得到不同的结果：

d <- d[rowSums(is.na(d[,3:6]))!=4,]

ID <- c(1,2,3,4,5,6,7,8,9,10)
Edu <- c(4,7,9,13,15,18,11,10,12,8)
CLL <- c(NA,0.33,0.45,NA,0.76,0.34,0.49,0.86,0.43,0.26)
BLL <- c(NA,0.93,0.15,NA,0.86,0.14,0.29,0.36,0.93,0.76)
ALL <- c(NA,0.53,0.65,NA,0.26,0.54,0.99,0.76,0.63,NA)
DLL <- c(NA,0.13,0.95,NA,0.36,0.74,0.19,NA,0.83,0.56)

d <- data.frame(ID, Edu, CLL, BLL, ALL, DLL)
d <- d[rowSums(is.na(d[,3:6]))!=4,]

cor.prob <- function (X, dfr = nrow(X) - 2) {
  R <- cor(X, use="complete.obs", method = c("pearson"))
  above <- row(R) < col(R)
  r2 <- R[above]^2
  Fstat <- r2 * dfr/(1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr)
  R[row(R) == col(R)] <- NA
  R
}

flattenSquareMatrix <- function(m) {
  if( (class(m) != "matrix") | (nrow(m) != ncol(m))) stop("Must be a square matrix.") 
  if(!identical(rownames(m), colnames(m))) stop("Row and column names must be equal.")
  ut <- upper.tri(m)
  data.frame(i = rownames(m)[row(m)[ut]],
             j = rownames(m)[col(m)[ut]],
             cor=t(m)[ut],
             p=m[ut])
}

kor_masterlist_d <- flattenSquareMatrix (cor.prob(d))
kor_masterlist_d_order <- kor_masterlist_d[order(-abs(kor_masterlist_d$cor)),]
kor_masterlist_d_Edu <- subset(kor_masterlist_d_order, i == "Edu" & j != "ID")
kor_masterlist_d_Edu
#>      i   j        cor         p
#> 8  Edu ALL -0.3319602 0.4217894
#> 3  Edu CLL  0.2646661 0.5264383
#> 12 Edu DLL  0.2609108 0.5325405
#> 5  Edu BLL -0.2492912 0.5515835

^{由reprex package (v2.0.0) 于 2021-07-09 创建}

ID <- c(1,2,3,4,5,6,7,8,9,10)
Edu <- c(4,7,9,13,15,18,11,10,12,8)
CLL <- c(NA,0.33,0.45,NA,0.76,0.34,0.49,0.86,0.43,0.26)
BLL <- c(NA,0.93,0.15,NA,0.86,0.14,0.29,0.36,0.93,0.76)
ALL <- c(NA,0.53,0.65,NA,0.26,0.54,0.99,0.76,0.63,NA)
DLL <- c(NA,0.13,0.95,NA,0.36,0.74,0.19,NA,0.83,0.56)

d <- data.frame(ID, Edu, CLL, BLL, ALL, DLL)

cor.prob <- function (X, dfr = nrow(X) - 2) {
  R <- cor(X, use="complete.obs", method = c("pearson"))
  above <- row(R) < col(R)
  r2 <- R[above]^2
  Fstat <- r2 * dfr/(1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr)
  R[row(R) == col(R)] <- NA
  R
}

flattenSquareMatrix <- function(m) {
  if( (class(m) != "matrix") | (nrow(m) != ncol(m))) stop("Must be a square matrix.") 
  if(!identical(rownames(m), colnames(m))) stop("Row and column names must be equal.")
  ut <- upper.tri(m)
  data.frame(i = rownames(m)[row(m)[ut]],
             j = rownames(m)[col(m)[ut]],
             cor=t(m)[ut],
             p=m[ut])
}

kor_masterlist_d <- flattenSquareMatrix (cor.prob(d))
kor_masterlist_d_order <- kor_masterlist_d[order(-abs(kor_masterlist_d$cor)),]
kor_masterlist_d_Edu <- subset(kor_masterlist_d_order, i == "Edu" & j != "ID")
kor_masterlist_d_Edu
#>      i   j        cor         p
#> 8  Edu ALL -0.3319602 0.3487063
#> 3  Edu CLL  0.2646661 0.4599218
#> 12 Edu DLL  0.2609108 0.4665494
#> 5  Edu BLL -0.2492912 0.4873203

^{由reprex package (v2.0.0) 于 2021-07-09 创建}

Answer 1

区别用

介绍

  Fstat <- r2 * dfr/(1 - r2)

由于dfr 的值不同（8 和 6）。

> r2
 [1] 0.234375000 0.011456074 0.070048129 0.002401998 0.062146106 0.062751663
 [7] 0.096834764 0.110197604 0.165400138 0.216547595 0.057255529 0.068074453
[13] 0.030140179 0.136955494 0.005027238
> r2*8/(1-r2)
 [1] 2.44897959 0.09271069 0.60259574 0.01926226 0.53011332 0.53562464
 [7] 0.85773686 0.99076023 1.58543173 2.21121379 0.48586255 0.58437675
[13] 0.24861473 1.26951037 0.04042111
> r2*6/(1-r2)
 [1] 1.83673469 0.06953302 0.45194680 0.01444669 0.39758499 0.40171848
 [7] 0.64330264 0.74307017 1.18907380 1.65841034 0.36439691 0.43828256
[13] 0.18646105 0.95213278 0.03031583

8 或 6 都不正确，因为它基于 nrow(X)，而对于 use="complete.obs"，它必须基于完整观察的数量。这可以通过将函数定义更改为cor.prob <- function (X, dfr = sum(complete.cases(X)) - 2) { … 来实现。因此，无论有无先验，都会产生相同的结果
d <- d[rowSums(is.na(d[,3:6]))!=4,]

但是如果我选择使用pairwise.complete.obs 而不是complete.obs，那么我会保持原样吗？而且，由于我的数据中有 NA 值不同的位置，那么我将受益于 "pairwise" 而不是 "complete.obs"？

确实，如果我们使用pairwise.complete.obs，我们会获得只有部分列是NA 的观察结果。但是，由于我们对各个列有不同数量的观察，因此单个 dfr 值是不合适的；相反，我们可以使用dfr 矩阵：

library(psych)
cor.prob <- function (X, dfr = pairwiseCount(X) - 2)
{
  …
  Fstat <- r2 * dfr[above]/(1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr[above])
  …