Матрица дисперсии-ковариации: найти ковариацию для уникальных комбинаций переменных * переменных - PullRequest
Новая
       94

Матрица дисперсии-ковариации: найти ковариацию для уникальных комбинаций переменных * переменных

0 голосов
/ 31 мая 2018

Цель

Используйте эти данные: 

dat<-structure(list(study = c(1, 1, 2, 2, 3, 4, 4, 5, 5), nrate = c(1,                                                                1, 1, 2, 1, 1, 1, 1, 2), trt = c(1, 2, 1, 1, 1, 1, 2, 1, 2),                n2i = c(25, 25, 40, 40, 50, 30, 30, 20, 30), Ni = c(75, 75,                                                                    80, 80, 100, 90, 90, 40, 60), yi = structure(c(1.75557336268135,                                                                                                                   1.16269114535263, 2.25236533601502, 1.65098691534697, 1.93238812372334,                                                                                                                   2.80537854506277, 2.47373334918987, 1.36964712768673, 1.18135471573816                                                                   ), measure = "ROM", ni = c(50, 50, 80, 80, 100, 60, 60, 40,                                                                                               60)), vi = c(0.0972473617680551, 0.10417464101422, 0.0525739144226032,                                                                                                            0.0135660003587117, 0.036197209164285, 0.341666364303935,                                                                                                            0.342935708755073, 0.0303744729767536, 0.00416144452369287                                                                                              )), .Names = c("study", "nrate", "trt", "n2i", "Ni", "yi",                                                                                                              "vi"), row.names = c(NA, -9L), class = c("escalc", "data.frame"                                                                                                             ), yi.names = "yi", vi.names = "vi", digits = 4)
 dat<-data.frame(dat)

Чтобы получить эту матрицу дисперсии-ковариации: 

      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7] [,8]  [,9]
 [1,] 0.097 0.054 0.000 0.000 0.000 0.000 0.000 0.00 0.000
 [2,] 0.054 0.104 0.000 0.000 0.000 0.000 0.000 0.00 0.000
 [3,] 0.000 0.000 0.053 0.000 0.000 0.000 0.000 0.00 0.000
 [4,] 0.000 0.000 0.000 0.014 0.000 0.000 0.000 0.00 0.000
 [5,] 0.000 0.000 0.000 0.000 0.036 0.000 0.000 0.00 0.000
 [6,] 0.000 0.000 0.000 0.000 0.000 0.342 0.072 0.00 0.000
 [7,] 0.000 0.000 0.000 0.000 0.000 0.072 0.343 0.00 0.000
 [8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.03 0.000
 [9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.004

Я хочу использоватьэта функция для вычисления ковариаций между "yi's" по уровню ("study" * "nrate") и оставления "vi" дисперсий, которые у меня есть, на диагоналях: 

 library(metafor) #bldiag comes from here
calc.v <- function(x) {
  v <- matrix(1/x$n2i[1] + outer(x$yi, x$yi, "*")/(2*x$Ni[1]), nrow=nrow(x), ncol=nrow(x))
  diag(v) <- x$vi
  v
}

А применение следующего к данным дает матрицу var-cov, которая почти есть *, если бы я только мог применить ее на уникальных уровнях «изучения» * «обучения», а не на уровне «изучения» и получить обратно матрицу.

V <- bldiag(lapply(split(dat, dat[,c("study")]), calc.v))

Проблема

Я пытался: 

V <- bldiag(lapply(split(dat, dat[,c("study","nrate")]), calc.v))

и 

V <- bldiag(lapply(unique(dat[,c("study","nrate")]), calc.v))

, которые дают ошибки Error in bldiag(lapply(split(dat, dat[, c("study", "nrate")]), calc.v)) : replacement has length zeroи Error in x$n2i : $ operator is invalid for atomic vectors о bldiag, а затем о моей calc.v функции соответственно.

Сноски

* Почти там матрица (сравните с вышеупомянутой): 

       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7] [,8]  [,9]
 [1,] 0.097 0.054 0.000 0.000 0.000 0.000 0.000 0.00 0.000
 [2,] 0.054 0.104 0.000 0.000 0.000 0.000 0.000 0.00 0.000
 [3,] 0.000 0.000 0.053 0.048 0.000 0.000 0.000 0.00 0.000
 [4,] 0.000 0.000 0.048 0.014 0.000 0.000 0.000 0.00 0.000
 [5,] 0.000 0.000 0.000 0.000 0.036 0.000 0.000 0.00 0.000
 [6,] 0.000 0.000 0.000 0.000 0.000 0.342 0.072 0.00 0.000
 [7,] 0.000 0.000 0.000 0.000 0.000 0.072 0.343 0.00 0.000
 [8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.03 0.070
 [9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.07 0.004

Также, очевидно, есть функция vcov для встроенных моделей, но я не понимаю, как это мне помогает.

1 Ответ

0 голосов
/ 31 мая 2018

Проблема в том, что мы, когда делаем split, есть все комбинации разделения, даже может быть элемент list с 0 строками.Чтобы удалить эти элементы, нам нужно использовать drop = TRUE

library(metafor)
V <- bldiag(lapply(split(dat, dat[,c("study","nrate")], drop = TRUE), calc.v))
dim(V)
#[1] 9 9

round(V, 3)
#       [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8]  [,9]
# [1,] 0.097 0.054 0.000 0.000 0.000 0.000 0.00 0.000 0.000
# [2,] 0.054 0.104 0.000 0.000 0.000 0.000 0.00 0.000 0.000
# [3,] 0.000 0.000 0.053 0.000 0.000 0.000 0.00 0.000 0.000
# [4,] 0.000 0.000 0.000 0.036 0.000 0.000 0.00 0.000 0.000
# [5,] 0.000 0.000 0.000 0.000 0.342 0.072 0.00 0.000 0.000
# [6,] 0.000 0.000 0.000 0.000 0.072 0.343 0.00 0.000 0.000
# [7,] 0.000 0.000 0.000 0.000 0.000 0.000 0.03 0.000 0.000
# [8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.014 0.000
# [9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.000 0.004
...