1

我的数据如下所示：

data = data.frame(a.coef = c(.14, .15, .16), 
                  b.coef = c(.4, .5, .6), 
                  a.var = c(0.0937, 0.0934, 0.0945), 
                  b.var = c(0.00453, 0.00564, 0.00624), 
                  ab.cov = c(0.000747, 0.000747, 0.000747))

我想在数据集的每一行上运行以下代码（来源： http ://www.quantpsy.org/medmc/medmc.htm）。

require(MASS)
a = data$a.coef
b = data$b.coef
rep = 10000
conf = 95
pest = c(a, b)
acov <- matrix(c(data$a.var, data$ab.cov, 
                 data$ab.cov, data$b.var), 2, 2)
mcmc <- mvrnorm(rep, pest, acov, empirical = FALSE)
ab <- mcmc[ , 1] * mcmc[ , 2]
low = (1 - conf / 100) / 2
upp = ((1 - conf / 100) / 2) + (conf / 100)
LL = quantile(ab, low)
UL = quantile(ab, upp)
LL4 = format(LL, digits = 4)
UL4 = format(UL, digits = 4)

我创建了一个相对简单的函数，它将数据和行号作为输入：

MCMAM <- function(data_input, row_number) {
  data = data_input[row_number, ]
  a = data[["a.coef"]]
  b = data[["b.coef"]]
  rep = 10000
  conf = 95
  pest = c(a, b)
  acov <- matrix(c(data[["a.var"]], data[["ab.cov"]], 
                   data[["ab.cov"]], data[["b.var"]]), 2, 2)
  require(MASS)
  mcmc <- mvrnorm(rep, pest, acov, empirical = FALSE)
  ab <- mcmc[, 1] * mcmc[, 2]
  low = (1 - conf / 100) / 2
  upp = ((1 - conf / 100) / 2) + (conf / 100)
  LL = quantile(ab, low)
  UL = quantile(ab, upp)
  return(c(LL, UL))
}

MCMAM(data, 1)

    2.5%      97.5% 
-0.1901272  0.3104614

但是，如果有一种方法可以摆脱行规范，让函数逐行遍历数据集并将输出保存到数据集中的新列，那就太好了。

我一直在试验 for 循环和应用函数，但没有取得任何成功，主要是因为 matrix() 和 mvrnorm() 函数都采用值而不是向量。

4

1 回答 1

3

我们可以用lapply

do.call(rbind, lapply(seq_len(nrow(data)), MCMAM, data_input = data))

-输出

    2.5%     97.5%
[1,] -0.1832449 0.3098362
[2,] -0.2260856 0.3856575
[3,] -0.2521126 0.4666583

或使用rowwise

library(dplyr)
library(tidyr)
data %>% 
     rowwise %>% 
     mutate(new = list(MCMAM(cur_data(), 1))) %>%
     unnest_wider(new)
# A tibble: 3 x 7
#  a.coef b.coef  a.var   b.var   ab.cov `2.5%` `97.5%`
#   <dbl>  <dbl>  <dbl>   <dbl>    <dbl>  <dbl>   <dbl>
#1   0.14    0.4 0.0937 0.00453 0.000747 -0.185   0.309
#2   0.15    0.5 0.0934 0.00564 0.000747 -0.219   0.396
#3   0.16    0.6 0.0945 0.00624 0.000747 -0.259   0.472
于 2021-06-01T22:04:57.030 回答