#r #dplyr #mass
Вопрос:
У меня есть данные, которые выглядят так:
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
Но было бы здорово, если бы существовал способ избавиться от спецификации строк и просто запустить функцию по строкам набора данных и сохранить выходные данные в новом столбце набора данных.
Я экспериментировал с циклами и применял функции, но не добился никакого успеха, в основном потому, что функции matrix() и mvrnorm() принимают значения, а не векторы.
Ответ №1:
Мы можем использовать 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