机器算法验证 - ICC 与同一组中两个随机抽取的单元之间的预期相关性 - 吾爱随笔录

ICC 与同一组中两个随机抽取的单元之间的预期相关性

机器算法验证相关性混合模式类内相关

2022-03-09 18:09:59

在多级建模中，通常从随机效应方差分析中计算出类内相关性

y_{i j} = γ_{00} + u_{j} + e_{i j}

$y_{ij} = \gamma_{00} + u_j + e_{ij}$

其中是二级残差，是一级残差。然后我们分别得到和方差的和，并将它们代入以下方程： $u_j$ $e_{ij}$ $\hat{\sigma}_u^2$ $\hat{\sigma}_e^2$ $u_j$ $e_{ij}$

ρ = \frac{{\hat{σ}}_{u}^{2}}{{\hat{σ}}_{u}^{2} + {\hat{σ}}_{e}^{2}}

$ρ = \frac{\hat{\sigma}_u^2}{\hat{\sigma}_u^2 +\hat{\sigma}_e^2}$

Hox (2002)在 p15上写道

类内相关性 ρ 也可以解释为同一组中两个随机抽取的单元之间的预期相关性

这里有一个问题，它提出了一个高级问题（为什么它完全等于这个而不是近似等于）并得到一个高级答案。

但是，我想问一个更简单的问题。

问题：谈论同一组中两个随机抽取的单位之间的相关性甚至意味着什么？

我对类内相关性适用于组而不是配对数据这一事实有一个基本的了解。但是，如果我们只有来自同一组的两个随机抽取的单元，我仍然不明白如何计算相关性。例如，如果我查看ICC 维基百科页面上的点图，我们有多个组和每个组内的多个点。

2个回答

如果您考虑每组只有两个人的情况，则可能最容易看到等价性。所以，让我们看一个具体的例子（我会用 R 来做这个）：

dat <- read.table(header=TRUE, text = "
group person   y
1     1        5
1     2        6
2     1        3
2     2        2
3     1        7
3     2        9
4     1        2
4     2        2
5     1        3
5     2        5
6     1        6
6     2        9
7     1        4
7     2        2
8     1        8
8     2        7")

所以，我们有 8 组，每组 2 个人。现在让我们拟合随机效应方差分析模型：

library(nlme)
res <- lme(y ~ 1, random = ~ 1 | group, data=dat, method="ML")

最后，让我们计算 ICC：

getVarCov(res)[1] / (getVarCov(res)[1] + res$sigma^2)

这产生：（0.7500003准确地说是 0.75，但这里的估计过程中有一些轻微的数字印象）。

现在让我们将数据从长格式重塑为宽格式：

dat <- as.matrix(reshape(dat, direction="wide", v.names="y", idvar="group", timevar="person"))

现在看起来像这样：

   group y.1 y.2
1      1   5   6
3      2   3   2
5      3   7   9
7      4   2   2
9      5   3   5
11     6   6   9
13     7   4   2
15     8   8   7

现在计算和之间的相关y.1性y.2：

cor(dat[,2], dat[,3])

这产生：0.8161138

等等，什么？这里发生了什么？不应该是0.75吗？不完全的！我上面计算的不是 ICC（类内相关系数），而是常规的 Pearson 积矩相关系数，它是一个类间相关系数。请注意，在长格式数据中，谁是人 1 和谁是人 2 完全是任意的——这些对是无序的。您可以重新排列组内的数据，您会得到相同的结果。y.1但在宽格式数据中，谁列在下面，谁列在下面并不是任意的y.2。如果你要切换一些个体，你会得到不同的相关性（除非你要切换所有个体——那么这相当于cor(dat[,3], dat[,2])这当然仍然给你0.8161138）。

费舍尔指出的是用宽格式数据获取 ICC 的一个小技巧。将每对按两个顺序包含两次，然后计算相关性：

dat <- rbind(dat, dat[,c(1,3,2)])
cor(dat[,2], dat[,3])

这产生：0.75.

所以，正如你所看到的，ICC 实际上是一个相关系数——对于来自同一组的两个人的“未配对”数据。

如果每个组有两个以上的人，您仍然可以以这种方式考虑 ICC，除了在组内创建个人配对的方式会更多。然后，ICC 是所有可能配对之间的相关性（再次以无序的方式）。

@Wolfgang 已经给出了很好的答案。我想稍微扩展一下，以表明您还可以通过逐字地实现随机选择多对值的直观算法——其中每对的成员来自同一组——然后简单地计算它们的相关性。然后这个相同的过程可以很容易地应用于具有任何大小的组的数据集，正如我还将展示的那样。 $y$

首先我们加载@Wolfgang 的数据集（此处未显示）。现在让我们定义一个简单的 R 函数，它接受一个 data.frame 并返回来自同一组的一对随机选择的观察值：

get_random_pair <- function(df){
  # select a random row
  i <- sample(nrow(df), 1)
  # select a random other row from the same group
  # (the call to rep() here is admittedly odd, but it's to avoid unwanted
  # behavior when the first argument to sample() has length 1)
  j <- sample(rep(setdiff(which(dat$group==dat[i,"group"]), i), 2), 1)
  # return the pair of y-values
  c(df[i,"y"], df[j,"y"])
}

这是一个例子，如果我们在@Wolfgang 的数据集上调用这个函数 10 次，我们会得到什么：

test <- replicate(10, get_random_pair(dat))
t(test)
#       [,1] [,2]
#  [1,]    9    6
#  [2,]    2    2
#  [3,]    2    4
#  [4,]    3    5
#  [5,]    3    2
#  [6,]    2    4
#  [7,]    7    9
#  [8,]    5    3
#  [9,]    5    3
# [10,]    3    2

现在要估计 ICC，我们只需多次调用此函数，然后计算两列之间的相关性。

random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
#           [,1]      [,2]
# [1,] 1.0000000 0.7493072
# [2,] 0.7493072 1.0000000

同样的过程可以应用到具有任何大小组的数据集，而无需任何修改。例如，让我们创建一个由 100 组每组 100 个观察值组成的数据集，真实的 ICC 设置为 0.75，如 @Wolfgang 的示例所示。

set.seed(12345)
group_effects <- scale(rnorm(100))*sqrt(4.5)
errors <- scale(rnorm(100*100))*sqrt(1.5)
dat <- data.frame(group = rep(1:100, each=100),
                  person = rep(1:100, times=100),
                  y = rep(group_effects, each=100) + errors)

stripchart(y ~ group, data=dat, pch=20, col=rgb(0,0,0,.1), ylab="group")

根据混合模型的方差分量估计 ICC，我们得到：

library("lme4")
mod <- lmer(y ~ 1 + (1|group), data=dat, REML=FALSE)
summary(mod)
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  group    (Intercept) 4.502    2.122   
#  Residual             1.497    1.223   
# Number of obs: 10000, groups:  group, 100

4.502/(4.502 + 1.497)
# 0.7504584

如果我们应用随机配对过程，我们得到

random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
#           [,1]      [,2]
# [1,] 1.0000000 0.7503004
# [2,] 0.7503004 1.0000000

这与方差分量估计非常吻合。

请注意，虽然随机配对过程有点直观，并且在教学上很有用，但@Wolfgang 说明的方法实际上要聪明得多。对于像这样大小为 100*100 的数据集，组内唯一配对的数量（不包括自配对）为 505,000——一个很大但不是天文数字——因此我们完全有可能计算相关性完全耗尽所有可能的配对集合，而不是需要从数据集中随机抽样。这是一个函数，用于检索具有任意大小组的一般情况下的所有可能配对：

get_all_pairs <- function(df){
  # do this for every group and combine the results into a matrix
  do.call(rbind, by(df, df$group, function(group_df){
    # get all possible pairs of indices
    i <- expand.grid(seq(nrow(group_df)), seq(nrow(group_df)))
    # remove self-pairings
    i <- i[i[,1] != i[,2],]
    # return a 2-column matrix of the corresponding y-values
    cbind(group_df[i[,1], "y"], group_df[i[,2], "y"])
  }))
}

现在，如果我们将此函数应用于 100*100 数据集并计算相关性，我们得到：

cor(get_all_pairs(dat))
#           [,1]      [,2]
# [1,] 1.0000000 0.7504817
# [2,] 0.7504817 1.0000000

这与其他两个估计非常吻合，并且与随机配对过程相比，计算速度要快得多，并且在方差较小的意义上也应该是更有效的估计。

其它你可能感兴趣的问题

上一篇指数分布族背后的基本原理是什么？下一篇贝叶斯究竟如何定义（或解释？）概率？