brms 包中存在 x 和 y 误差的线性回归

Question

我正在尝试对我的数据进行线性回归，以计算出海平面变化的速率。但是，简单的线性回归将不起作用，因为我同时有 x（年龄）和 y（RSL）错误，例如：

RSL	RSL error	Age	Age error
-0.31	0.05	1815	1
-0.29	0.07	1880	5
-0.29	0.05	1895	5
-0.2	0.05	1935	1

我一直在做一些研究，看起来变量误差方法或贝叶斯测量模型都可以工作https://www.r-bloggers.com/2021/04/how-to-estimate-models-with-measurement-error-for-our-covid-19-indices/

我决定从贝叶斯测量模型开始，因为作者将其描述为更有利且更易于实施的模型。

我试图用我自己的数据复制他们的示例，但是我收到以下错误Error: The following variables can neither be found in 'data' nor in 'data2': 'Wap'

有谁知道我哪里出了问题以及如何让模型运行？

注意在我的数据框中，我有 Ageupper 和 Agelower 以及 RSLupper RSLlower，但它们是高斯​​的，所以我只在代码中使用 Ageupper RSLupper 等。

谢谢

## Load packages

library(brms)

### Load csv

Wap<-read.csv("Wapengo.csv",header=TRUE)

### Set errors

brms_formula<-bf(Wap~
me(RSL,RSLupper)+
me(Age,Ageupper),
center=TRUE)+
set_mecor(FALSE)

### Run model


model <- brm(brms_formula,
                data=Wap,
                silent=0,
                chains=1,save_pars = save_pars(),
                iter=500,
                warmup=250,
                backend='rstan')
## Wap Data

structure(list(Site = structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Label = "Wapengo", class = "factor"), 
    RSL = c(-0.068238463, -0.073155693, -0.02581141, -0.017379805, 
    -0.014178649, 0.026706959), RSLupper = c(0.16795545, 0.168146638, 
    0.16916378, 0.16951953, 0.168921232, 0.168238356), RSLlower = c(0.16795545, 
    0.168146638, 0.16916378, 0.16951953, 0.168921232, 0.168238356
    ), Age = c(1832L, 1860L, 1881L, 1894L, 1906L, 1913L), Ageupper = c(14.09253495, 
    13.7156267, 12.99997671, 12.25404364, 10.13081851, 10.19587526
    ), Agelower = c(14.09253495, 13.7156267, 12.99997671, 12.25404364, 
    10.13081851, 10.19587526), Rate = c(-0.037244426, -0.174854911, 
    2.332632731, 0.61776332, 0.279128313, 5.371978495)), row.names = c(NA, 
6L), class = "data.frame")

### Session info

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.16

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] brms_2.15.0     Rcpp_1.0.6      patchwork_1.1.1 dplyr_1.0.7     tidypaleo_0.1.1
[6] ggplot2_3.3.5  

loaded via a namespace (and not attached):
  [1] nlme_3.1-144         matrixStats_0.60.0   xts_0.12.1           threejs_0.3.3       
  [5] rstan_2.21.2         backports_1.1.7      tools_3.6.3          utf8_1.2.1          
  [9] R6_2.5.0             DT_0.13              DBI_1.1.0            mgcv_1.8-36         
 [13] projpred_2.0.2       colorspace_2.0-2     withr_2.4.2          prettyunits_1.1.1   
 [17] processx_3.4.5       tidyselect_1.1.1     gridExtra_2.3        Brobdingnag_1.2-6   
 [21] curl_4.3             compiler_3.6.3       cli_2.5.0            shinyjs_2.0.0       
 [25] labeling_0.4.2       colourpicker_1.1.0   scales_1.1.1         dygraphs_1.1.1.6    
 [29] mvtnorm_1.1-1        callr_3.5.1          ggridges_0.5.2       StanHeaders_2.21.0-7
 [33] stringr_1.4.0        digest_0.6.27        minqa_1.2.4          base64enc_0.1-3     
 [37] pkgconfig_2.0.3      htmltools_0.5.1.1    sessioninfo_1.1.1    lme4_1.1-23         
 [41] fastmap_1.1.0        htmlwidgets_1.5.3    rlang_0.4.11         rstudioapi_0.13     
 [45] shiny_1.6.0          farver_2.1.0         generics_0.1.0       jsonlite_1.7.2      
 [49] zoo_1.8-8            crosstalk_1.1.0.1    gtools_3.9.2         inline_0.3.19       
 [53] magrittr_2.0.1       loo_2.4.1            bayesplot_1.8.1      Matrix_1.2-18       
 [57] munsell_0.5.0        fansi_0.5.0          abind_1.4-5          lifecycle_1.0.0     
 [61] stringi_1.6.2        MASS_7.3-51.5        pkgbuild_1.2.0       plyr_1.8.6          
 [65] ggstance_0.3.4       grid_3.6.3           blob_1.2.1           parallel_3.6.3      
 [69] promises_1.1.0       crayon_1.4.1         miniUI_0.1.1.1       lattice_0.20-38     
 [73] splines_3.6.3        ps_1.5.0             pillar_1.6.1         igraph_1.2.6        
 [77] boot_1.3-24          markdown_1.1         shinystan_2.5.0      codetools_0.2-16    
 [81] reshape2_1.4.4       stats4_3.6.3         rstantools_2.1.1     glue_1.4.2          
 [85] V8_3.4.2             RcppParallel_5.1.4   vctrs_0.3.8          nloptr_1.2.2.1      
 [89] httpuv_1.6.1         gtable_0.3.0         purrr_0.3.4          tidyr_1.1.3         
 [93] assertthat_0.2.1     mime_0.9             xtable_1.8-4         coda_0.19-4         
 [97] later_1.0.0          rsconnect_0.8.18     tibble_3.1.2         shinythemes_1.2.0   
[101] gamm4_0.2-6          statmod_1.4.34       ellipsis_0.3.2       bridgesampling_1.1-2

Answer 1

我认为您需要的信息可以在这里找到：https://github.com/paul-buerkner/brms/issues/643

为了确保其他一切正常，我使用下面的代码拟合模型。我遇到了许多不同的转换，以及一些超过最大树深度的转换。我不确定您发布的数据是否是您的所有数据，但如果您遇到相同的问题，brms 和 Stan 警告消息通常非常有助于为您指明正确的方向。我也会考虑设置比这更好的先验，但模型本身的公式应该保持不变。

m1 <- brm(RSL | mi(RSLupper) ~ me(Age, Ageupper), data = Wap, save_mevars = TRUE)