使用 plyr 计算组的样本协方差矩阵答案

【问题标题】：Calculating a Sample Covariance Matrix for Groups with plyr使用 plyr 计算组的样本协方差矩阵
【发布时间】：2010-04-28 07:13:09
【问题描述】：

我将使用来自http://gettinggeneticsdone.blogspot.com/2009/11/split-apply-and-combine-in-r-using-plyr.html 的示例代码作为这个示例。所以，首先，让我们复制他们的示例数据：

mydata=data.frame(X1=rnorm(30), X2=rnorm(30,5,2),
SNP1=c(rep("AA",10), rep("Aa",10), rep("aa",10)),
SNP2=c(rep("BB",10), rep("Bb",10), rep("bb",10)))

在本例中，我将忽略 SNP2，并假设 SNP1 中的值表示组成员身份。那么，我可能想要一些关于 SNP1 中每个组的汇总统计信息：“AA”、“Aa”、“aa”。

然后，如果我想计算每个变量的均值，使用（稍微修改他们的代码）是有意义的：

> ddply(mydata, c("SNP1"), function(df)
data.frame(meanX1=mean(df$X1), meanX2=mean(df$X2)))
  SNP1      meanX1   meanX2
1   aa  0.05178028 4.812302
2   Aa  0.30586206 4.820739
3   AA -0.26862500 4.856006

但是如果我想要每个组的样本协方差矩阵怎么办？理想情况下，我想要一个 3D 数组，其中我有每个组的协方差矩阵，第三维表示相应的组。我尝试了之前代码的修改版本，得到了以下结果，这让我确信我做错了什么。

> daply(mydata, c("SNP1"), function(df) cov(cbind(df$X1, df$X2)))
, ,  = 1


SNP1         1          2
  aa 1.4961210 -0.9496134
  Aa 0.8833190 -0.1640711
  AA 0.9942357 -0.9955837

, ,  = 2


SNP1          1        2
  aa -0.9496134 2.881515
  Aa -0.1640711 2.466105
  AA -0.9955837 4.938320

我在想第 3 维的 dim() 应该是 3，但实际上是 2。实际上，这是每个组的协方差矩阵的切片版本。如果我们手动计算 aa 的样本协方差矩阵，我们得到：

           [,1]       [,2]
[1,]  1.4961210 -0.9496134
[2,] -0.9496134  2.8815146

使用 plyr，下面给出了我想要的 list() 形式：

> dlply(mydata, c("SNP1"), function(df) cov(cbind(df$X1, df$X2)))
$aa
           [,1]       [,2]
[1,]  1.4961210 -0.9496134
[2,] -0.9496134  2.8815146

$Aa
           [,1]       [,2]
[1,]  0.8833190 -0.1640711
[2,] -0.1640711  2.4661046

$AA
           [,1]       [,2]
[1,]  0.9942357 -0.9955837
[2,] -0.9955837  4.9383196

attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
  SNP1
1   aa
2   Aa
3   AA

但就像我之前所说的，我真的很喜欢这个 3D 数组。关于我在 daply() 或建议上哪里出错的任何想法？当然，我可以将 dlply() 中的列表类型转换为 3D 数组，但我宁愿不这样做，因为我将在模拟中多次重复此过程。

作为旁注，我找到了一种方法 (http://www.mail-archive.com/r-help@r-project.org/msg86328.html)，它为每个组提供样本协方差矩阵，但输出的对象是臃肿的。

提前致谢。

【问题讨论】：

标签： r

【解决方案1】：

daply 使拆分变量成为数组中的第一个维。

a <- daply(mydata, c("SNP1"), function(df) cov(cbind(df$X1, df$X2)))
l <- dlply(mydata, c("SNP1"), function(df) cov(cbind(df$X1, df$X2)))

这样a[1, , ] 和l[[1]] 对应于相同的输出。正如 wkmor1 所建议的那样，您可以使用 aperm 重新排列尺寸，但我想了解更多关于为什么初始形式不适合您的需求。

【讨论】：

Hadley，我只是对上面显示的输出感到困惑。我没有意识到我正确使用了 daply() 。我没有正确完成的是使用 a[1, , ] 而不是 a[, , 1]。我期待使用 [ , , 1]。

【解决方案2】：

怎么样...

aperm(daply(mydata, c("SNP1"), function(df) cov(cbind(df$X1, df$X2))),perm=c(2,3,1))

'apem' 是数组，就像 't' 是矩阵。 perm 参数指定了 dims 应该改变的方式。

【讨论】：