Variance-covariance matrix: find covariance for unique variable*variable combinations

Question

Goal

Use this data:

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)

To get this variance-covariance matrix:

      [,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

I want to use this function to calculate covariances between "yi's" by level of ("study"*"nrate") and leave the "vi" variances I have alone on the diagonals:

 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
}

And applying the following to the data gives a var-cov matrix that is almost there*, if only I could only apply it at unique "study"*"nrate" levels rather than at the "study" level and get back a matrix.

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

Problem

I tried:

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

and

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

which give errors Error in bldiag(lapply(split(dat, dat[, c("study", "nrate")]), calc.v)) : replacement has length zero and Error in x$n2i : $ operator is invalid for atomic vectors about the bldiag and then my calc.v function, respectively.

Footnotes

*Almost there matrix (compare with above):

       [,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

Also apparently there is a vcov function for fitted models but I don't see how that helps me.

Sorry, I thought your input was a matrix by looking at the covariance matrix. It is not. The issue is described in the solution — akrun

akrun akrun · Accepted Answer · 2018-05-31T15:03:20

The issue is that we when do the split, there are all combinations of split happens, even there can be list element with 0 rows. In order to remove those elements, we need to use 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

Variance-covariance matrix: find covariance for unique variable*variable combinations

1 Answers