Linear regression with both x and y errors in package brms

Question

I am trying to run linear regressions on my data to work out the rate of sea-level change. However, a simple linear regression will not work as I have both x (Age) and y (RSL) errors for example:

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

I have been doing some research and it looks like either an error-in-variables approach or a Bayesian measurement model will work https://www.r-bloggers.com/2021/04/how-to-estimate-models-with-measurement-error-for-our-covid-19-indices/

I decided to start with the Bayesian measurement model as the authors describe this as the more advantageous and easier to implement model.

I tried to replicate their example with my own data however I get the following error Error: The following variables can neither be found in 'data' nor in 'data2': 'Wap'

Does anyone know where I am going wrong and how I can get the model to run?

N.B. In my dataframe I have Ageupper and Agelower and RSLupper RSLlower but they are gaussian so I just use Ageupper RSLupper etc in the code.

Thanks


## 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

I do not use this package, but in the data you're sharing there is not a `Wap` variable, that are you calling in the `brms_formula`, which is called in `brm()` . Maybe the error is there but I am guessing. — s__, Jul 29 '21 at 13:11
Hi S_ yes I thought it could be this but when I put ```brms_formula``` in it wouldn't work because data= needs to be a data frame rather than a list which makes me think it must need to come from the original data frame — Sophie Williams, Jul 29 '21 at 13:14
What is the response you're trying to predict? :) Is it `Rate`? If so, I'd replace `Wap ~` with `Rate ~` — Jonny Phelps, Jul 29 '21 at 13:20
Hi Jonny thanks for your comment. Yes i am trying to predict the rate using the RSL and Age for example Rate=(RSL1-RSL2)/(AGE1/AGE2) if that makes sense? — Sophie Williams, Jul 29 '21 at 13:29
Sorry I should mention that the ```rate``` column is just exactly that and doesn't take into account the errors - it is just to see a roundabout estimate of what it should be. — Sophie Williams, Jul 29 '21 at 13:32
I believe that the following from package ```deming``` may do the same thing: library(deming) fit <- deming(RSL ~ Age, data=Wap, xstd=Ageupper, ystd=RSLupper) print(fit) — Sophie Williams, Jul 29 '21 at 14:01

sjp · Answer 1 · 2021-08-03T10:41:03.120

0

I think the information you need can be found here: https://github.com/paul-buerkner/brms/issues/643

Just to make sure everything else worked, I fit the model using the code below. I ran into a number of divergent transitions, as well as a number of transitions that exceeded that maximum treedepth. I'm unsure if the data you posted was all of your data or not, but if you run into the same problems the brms and Stan warning messages are normally very helpful in pointing you in the right direction. I'd also look into setting better priors than this, but the formula for the model itself should stay the same.

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

edited Aug 03 '21 at 10:41

answered Jul 29 '21 at 22:00

sjp

820
1
5
10

Hi sjp, thanks for your help. Unfortunately i am still getting the error Error: The following variables can neither be found in 'data' nor in 'data2': 'RSLUpper' when I run your line of code – Sophie Williams Aug 02 '21 at 07:53
Sorry, there was a typo in my code. I had `RSLUpper` when it should have been `RSLupper` – sjp Aug 03 '21 at 10:40
Hi there, thanks for this. Running this gives me a new error though! Any idea what this error is? ```SAMPLING FOR MODEL ‘d7e503b7b07e1f8b09c966d1d3ccec06’ NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 2.5e-05 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup) [1] “Error in sampler$call_sampler(args_list[[i]]) : " " c++ exception (unknown reason)” error occurred during calling the sampler; sampling not done``` – Sophie Williams Aug 04 '21 at 09:40
Can you give the `sessionInfo` output? This in an issue that has been solved for others: https://github.com/stan-dev/rstan/issues/716 – sjp Aug 05 '21 at 10:03
Hi there I have included ```sessionInfo``` into the original comment now - thanks for attaching that link, I will have a read now – Sophie Williams Aug 05 '21 at 13:25
I want to *tentatively* say this has been solved by installing both the latest version of R and R studio and reinstalling stan from CRAN - new error messages but it at least seems to run now – Sophie Williams Aug 05 '21 at 14:14

Linear regression with both x and y errors in package brms

1 Answers1