Я собираюсь использовать пример кода из 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
Но что, если мне нужна выборочная ковариационная матрица для каждой группы? В идеале я хотел бы получить трехмерный массив, где у меня есть ковариационная матрица для каждой группы, а третье измерение обозначает соответствующую группу. Я попробовал модифицированную версию предыдущего кода и получил следующие результаты, которые убедили меня, что я делаю что-то не так.
> 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
Я думал, что dim () 3-го измерения будет равен 3, но вместо этого равен 2. На самом деле это разрезанная версия ковариационной матрицы для каждой группы. Если мы вручную вычисляем образец ковариационной матрицы для аа, мы получим:
[,1] [,2]
[1,] 1.4961210 -0.9496134
[2,] -0.9496134 2.8815146
Используя plyr, следующее дает мне то, что я хочу в форме списка ():
> 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 () к трехмерному массиву, но я бы предпочел этого не делать, потому что я буду повторять этот процесс много раз в симуляции.
В качестве дополнительного примечания я нашел один метод (http://www.mail-archive.com/r-help@r-project.org/msg86328.html), который предоставляет образец ковариационной матрицы для каждой группы , но выводимый объект раздутый.
Заранее спасибо.