R语言 变量组分块对角协方差矩阵

hrysbysz  于 2023-02-27  发布在  其他
关注(0)|答案(1)|浏览(150)

如果这是我的数据集,3个受试者,在两个时间点t1、t2测量,在每个时间点每个受试者测量两次

ID         Time     Observation    Y
     1          t1       1              3.1
     1          t1       2              4.5
     1          t2       1              4.2
     1          t2       2              1.7

     2          t1       1              2.3
     2          t1       2              3.1
     2          t2       1              1.5
     2          t2       2              2.2

     3          t1       1              4.1
     3          t1       2              4.9
     3          t2       1              3.5
     3          t2       2              3.2

这个数据集的块对角协方差矩阵是什么?我假设这个矩阵有三个块,每个块代表受试者i在t1和t2时的2x2方差-协方差矩阵。

j8ag8udp

j8ag8udp1#

下面的代码计算每个ID的协方差矩阵,然后用它们创建一个分块矩阵。

make_block_matrix <- function(x, y, fill = NA) {
  u <- matrix(fill, nrow = dim(x)[1], ncol = dim(y)[2])
  v <- matrix(fill, nrow = dim(y)[1], ncol = dim(x)[2])
  u <- cbind(x, u)
  v <- cbind(v, y)
  rbind(u, v)
}

# compute covariance matrices for each ID
sp <- split(df1, df1$ID)
cov_list <- lapply(sp, \(x) {
  y <- reshape(x[-1], direction = "wide", idvar = "Time", timevar = "Observation")
  cov(y[-1])
})
# no longer needed, tidy up
rm(sp)

# fill with NA's, simple call
# res <- Reduce(make_block_matrix, cov_list)

# fill with zeros
res <- Reduce(\(x, y) make_block_matrix(x, y, fill = 0), cov_list)
# finally, the results' row and column names
nms <- paste(sapply(cov_list, rownames), names(cov_list), sep = ".")
dimnames(res) <- list(nms, nms)
res
#>        Y.1.1 Y.2.2 Y.1.3 Y.2.1 Y.1.2 Y.2.3
#> Y.1.1  0.605 -1.54  0.00 0.000  0.00 0.000
#> Y.2.2 -1.540  3.92  0.00 0.000  0.00 0.000
#> Y.1.3  0.000  0.00  0.32 0.360  0.00 0.000
#> Y.2.1  0.000  0.00  0.36 0.405  0.00 0.000
#> Y.1.2  0.000  0.00  0.00 0.000  0.18 0.510
#> Y.2.3  0.000  0.00  0.00 0.000  0.51 1.445

创建于2023年2月20日,使用reprex v2.0.2

数据

df1 <- "ID         Time     Observation    Y
     1          t1       1              3.1
     1          t1       2              4.5
     1          t2       1              4.2
     1          t2       2              1.7

     2          t1       1              2.3
     2          t1       2              3.1
     2          t2       1              1.5
     2          t2       2              2.2

     3          t1       1              4.1
     3          t1       2              4.9
     3          t2       1              3.5
     3          t2       2              3.2"
df1 <- read.table(text = df1, header = TRUE)

创建于2023年2月20日,使用reprex v2.0.2

相关问题