Measurement error models using me()

Interfaces brms
,
1

Hi,
I’m dealing with a very simple model that has measurement error in the DV and IV. I’m almost sure that this was working before. But now it’s throwing divergent errors and the estimates don’t make sense.

This is the data:

df_VOTmandarin <- structure(list(subject = c("F01", "F02", "F03", "F04", "F05", 
"F06", "F07", "F08", "F09", "F10", "M11", "M12", "M13", "M14", 
"M15", "M16", "M17", "M18", "M19", "M20"), meanVOT = c(105.7, 
86.7, 97.8, 84.9, 84.6, 98.6, 82, 108, 91, 84.2, 66.9, 72.1, 
79.4, 81.4, 79.4, 92, 62.3, 79.1, 83, 91.5), meanvdur = c(160.240500790632, 
177.088580577028, 166.431480766486, 191.743223448853, 164.443296364447, 
209.664984526471, 164.963891458981, 178.764459628451, 164.695243878764, 
165.454699023242, 145.064740184494, 156.549546592339, 157.902925276459, 
186.697506183006, 180.052870601759, 137.171787564198, 146.799319226475, 
145.136248276608, 144.52767543895, 167.64967057907), sevdur = c(11.6375246514992, 
13.6695923443727, 14.7210641393673, 13.3598869761091, 13.033972044188, 
13.8418943147466, 13.3182550042139, 16.8821184127746, 17.1090382793098, 
13.931131552853, 11.9578088505664, 11.6801297648966, 12.9319905990126, 
11.1184014826944, 17.9909366691651, 9.27778876792764, 15.0041274537309, 
13.996813051598, 13.2763776071296, 15.3970331017383), seVOT = c(3.79195053882417, 
4.15812457725836, 4.62313025268955, 4.67724990423504, 4.49246282368839, 
4.09932243723819, 3.89301368550838, 8.34132949701532, 5.16397779494322, 
5.98479555020702, 4.69858134617957, 4.16719996587103, 3.82738669184195, 
3.81284379142101, 4.70744091837593, 7.21572357194112, 2.08193286261696, 
3.49745939536433, 5.81950742474539, 2.54405625374562), c_meanvdur = c(-5.31163172870365, 
11.5364480576924, 0.879348247150347, 26.1910909295173, -1.10883615488865, 
44.1128520071353, -0.588241060354648, 13.2123271091153, -0.856888640571668, 
-0.0974334960936574, -20.4873923348416, -9.00258592699666, -7.64920724287666, 
21.1453736636703, 14.5007380824233, -28.3803449551376, -18.7528132928607, 
-20.4158842427277, -21.0244570803856, 2.09753805973435)), row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20"), class = "data.frame")

And this is the model with its priors.

priors_me <- c(
  prior(normal(0, 200), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(normal(100, 50), class = meanme),
  prior(normal(0, 50), class = sdme),
  prior(normal(0, 50), class = sigma)
)
fit_mvotme <- brm(meanVOT | resp_se(seVOT, sigma = TRUE) ~
                    me(meanvdur, sevdur),
                  data = df_VOTmandarin,
                  family = gaussian(),
                  prior = priors_me,
                  control = list(adapt_delta = .999)
                  )

Did something change in the last (couple of) brms versions? Am I doing something wrong? Should I just be using mi()?

Thanks

  • Operating System: Ubuntu 20.04.2 LTS
  • rstan_2.21.3
  • StanHeaders_2.21.0-7
  • brms_2.14.4
2

I am sorry you are experiencing these problems. I don’t think something has changed recently from the brms side. If set up correctly mi() will be equivalent to me() in this case, but that his not going to help you I suppose.

3

Ok, but do you spot an error? Why might be that this simple model is not working?

4

I don’t see an error unfortunately but i didn’t check the appropriateness of the prior.

5

Along with Paul’s comments, the priors seem very broad. Have you done prior predictive checks? Also, you used adapt_delta=0.999 (perhaps already before you tested this again?) This indicates that the sampler is struggling a bit. Have you tried this with cmdstanr and the latest cmdstan? brms allows for using cmdstan as a backend.

If I use slightly more modest priors and the latest cmdstan then it samples with adapt_delta = 0.95.

In short, prior predictive checks is probably what needs to be done :)

[EDIT]: Here are the priors I used

priors_me <- c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(normal(100, 10), class = meanme),
    prior(normal(0, 10), class = sdme),
    prior(normal(0, 10), class = sigma)
)
6

But the true intercept is around 50.

Ok I’ll debug a bit more and report back…
I thought that something changed in brms because I was sure that the model converged a couple of years ago…

7

ok, if I increase the max_treedepth it works fine… I also rewrote the code in Stan and I get more or less the same output.
I’m still puzzled by how come the intercept goes to 0 with a huge CrI, but it must be the lack of data…

8

the intercept behavior is likely because of the scale and mean of the predictor in combination with the amount of measurement error.