在相同的数据上获得不同的 cor() 输出？
Getting different cor() output on same data?
我在计算相关性时遇到了问题。
我有两个数据集：
d3：一个有 1259 个观察值和 264 个变量。
d4：一个有 1196 个观察值和 185 个变量——这些只是血液测试。
两个数据都有相同的参与者唯一 ID，因此可以合并它们。
合并时（称为：d）我有 1150 个观察值（因为没有验血的数据已经出来了，验血数据有很多对照样本等等。当删除没有信息 x 的参与者时（因为这是纳入标准）我们最终得到 1144 个观察值。
在 d4 中，我有两个样本在所有变量中都缺少值。他们被包括在1144名参与者中。
所以在 sum datamanagement 之后我想计算信息 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 值略有不同。但还是一样的排名。
希望有人能帮助我。无论是否包括参与者，我都应该得到相同的结果。
我尝试制作以下可重现的小示例来说明我的问题。如果你 运行 它 with/without:
你会得到不同的结果
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 创建
与
的区别介绍
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,]
But if I choose to use pairwise.complete.obs instead of complete.obs, then i'll keep the code as it was? And further, since i have NA values different places in my data then i will benefite from "pairwise" rather than "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])
…
我在计算相关性时遇到了问题。 我有两个数据集： d3：一个有 1259 个观察值和 264 个变量。 d4：一个有 1196 个观察值和 185 个变量——这些只是血液测试。 两个数据都有相同的参与者唯一 ID，因此可以合并它们。
合并时（称为：d）我有 1150 个观察值（因为没有验血的数据已经出来了，验血数据有很多对照样本等等。当删除没有信息 x 的参与者时（因为这是纳入标准）我们最终得到 1144 个观察值。
在 d4 中，我有两个样本在所有变量中都缺少值。他们被包括在1144名参与者中。
所以在 sum datamanagement 之后我想计算信息 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 值略有不同。但还是一样的排名。
希望有人能帮助我。无论是否包括参与者，我都应该得到相同的结果。
我尝试制作以下可重现的小示例来说明我的问题。如果你 运行 它 with/without:
你会得到不同的结果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 创建与
的区别介绍 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,]
But if I choose to use
pairwise.complete.obsinstead ofcomplete.obs, then i'll keep the code as it was? And further, since i haveNAvalues different places in my data then i will benefite from"pairwise"rather than"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])
…