r-在相同的数据上获得不同的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

现在的问题是：我使用"；complete.obs"；在cor.((中。两名参与者的所有血液测试都缺失(NA(，如果我通过以下代码将他们从数据中删除：

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

然后，我从相关性中得到不同的结果。我最后少做了两次血液测试，这是一个显著的结果，p值超过0.05。当这两位参与者的价值观缺失时，我不明白其中的区别。

当使用"；"对等"；ind-cor((我对1142或1144名参与者得到了相同的结果。但最后4次血液测试与我用"；complete.obs"；

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

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

我试着做下面这个可重复的小例子来说明我的问题。如果你使用/不使用运行它，你会得到不同的结果

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包(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包(v2.0.0(于2021-07-09创建}

Answer 1

引入了差异

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

由于CCD_ 1(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"而不是dfr0中受益？

事实上，如果我们使用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])
…