如何将矩阵数据输入到brms公式中?

How to input matrix data into brms formula?

我正在尝试将矩阵数据输入 brm() 函数以 运行 进行信号回归。 brm 来自 brms 包,它提供了一个使用 Stan 拟合贝叶斯模型的接口。信号回归是指您在更大的模型中使用另一个协变量对一个协变量建模,并使用 by 参数,如下所示:model <- brm(response ~ s(matrix1, by = matrix2) + ..., data = Data)。问题是,我无法使用 'data' 参数输入我的矩阵,因为它只允许输入一个 data.frame 对象。

这是我的代码以及我在尝试绕过该约束时遇到的错误...

首先,我的模型构建前的可重现代码:

library(brms)
#100 rows, 4 columns. Each cell contains a number between 1 and 10
Data <- data.frame(runif(100,1,10),runif(100,1,10),runif(100,1,10),runif(100,1,10))
#Assign names to the columns
names(Data) <- c("d0_10","d0_100","d0_1000","d0_10000")
Data$Density <- as.matrix(Data)%*%c(-1,10,5,1)
#the coefficients we are modelling
d <- c(-1,10,5,1) 
#Made a matrix with 4 columns with values 10, 100, 1000, 10000 which are evaluation points. Rows are repeats of the same column numbers
Bins <- 10^matrix(rep(1:4,times = dim(Data)[1]),ncol = 4,byrow =T)
Bins

如上所述,由于'data'只允许输入一个data.frame对象,我尝试了其他输入矩阵数据的方法。这些方法包括:

1) 使用 as.matrix()

在 brm() 函数中创建矩阵
signalregression.brms <- brm(Density ~ s(Bins,by=as.matrix(Data[,c(c("d0_10","d0_100","d0_1000","d0_10000"))])),data = Data)
#Error in is(sexpr, "try-error") : 
  argument "sexpr" is missing, with no default

2) 在公式外创建矩阵,将其存储在一个变量中,然后在 brm() 函数中调用该变量

Donuts <- as.matrix(Data[,c(c("d0_10","d0_100","d0_1000","d0_10000"))])
signalregression.brms <- brm(Density ~ s(Bins,by=Donuts),data = Data)
#Error: The following variables can neither be found in 'data' nor in 'data2':
'Bins', 'Donuts'

3) 使用 'data2' 参数输入包含矩阵的列表

signalregression.brms <- brm(Density ~ s(Bins,by=donuts),data = Data,data2=list(Bins = 10^matrix(rep(1:4,times = dim(Data)[1]),ncol = 4,byrow =T),donuts=as.matrix(Data[,c(c("d0_10","d0_100","d0_1000","d0_10000"))])))
#Error in names(dat) <- object$term : 
  'names' attribute [1] must be the same length as the vector [0]
以上的

None 有效;每个人都有自己的错误,很难对它们进行故障排除,因为我无法在 brms 的上下文中找到具有相似性质的在线答案或示例。

我能够在 mgcv 程序包中对 gam() 使用上述技术——您不必使用 'data' 定义 data.frame,您可以调用变量在 gam() 公式之外定义,您可以在 gam() 函数本身内部创建矩阵。见下文:

library(mgcv)
signalregression2 <- gam(Data$Density ~ s(Bins,by = as.matrix(Data[,c("d0_10","d0_100","d0_1000","d0_10000")]),k=3))
#Works!

似乎 brms 不太灵活...:(

我的问题:有人对如何使我的 brm() 函数有任何建议吗运行?

非常感谢!

我对信号回归的理解非常有限,我不相信这是正确的,但我认为这至少是朝着正确方向迈出的一步。问题似乎是 brm() 期望其公式中的所有内容都是 data 中的一列。因此,我们可以通过确保我们想要的所有内容都存在于 data:

中来编译模型
library(tidyverse)
signalregression.brms = brm(Density ~
                              s(cbind(d0_10_bin, d0_100_bin, d0_1000_bin, d0_10000_bin),
                                by = cbind(d0_10, d0_100, d0_1000, d0_10000),
                                k = 3),
                            data = Data %>%
                              mutate(d0_10_bin = 10,
                                     d0_100_bin = 100,
                                     d0_1000_bin = 1000,
                                     d0_10000_bin = 10000))

手写每一列有点烦人;我相信还有更通用的解决方案。

作为参考,这是我安装的软件包版本:

map_chr(unname(unlist(pacman::p_depends(brms)[c("Depends", "Imports")])), ~ paste(., ": ", pacman::p_version(.), sep = ""))
 [1] "Rcpp: 1.0.6"           "methods: 4.0.3"        "rstan: 2.21.2"         "ggplot2: 3.3.3"       
 [5] "loo: 2.4.1"            "Matrix: 1.2.18"        "mgcv: 1.8.33"          "rstantools: 2.1.1"    
 [9] "bayesplot: 1.8.0"      "shinystan: 2.5.0"      "projpred: 2.0.2"       "bridgesampling: 1.1.2"
[13] "glue: 1.4.2"           "future: 1.21.0"        "matrixStats: 0.58.0"   "nleqslv: 3.3.2"       
[17] "nlme: 3.1.149"         "coda: 0.19.4"          "abind: 1.4.5"          "stats: 4.0.3"         
[21] "utils: 4.0.3"          "parallel: 4.0.3"       "grDevices: 4.0.3"      "backports: 1.2.1"