Correlated measurement errors in errors-in-variables models

[David West]

Hi Paul,

Is it currently possible to fit errors-in-variables models with correlated measurement error in brms?

Thus for example: brm(y | se(sei, sigma = TRUE) ~ me(bx, sdx), data = dat), where values in me(bx, sdx) are correlated.

Thanks a lot!

David

[Paul Buerkner]

correlated with what?

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+...@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/ebc6012e-449b-4420-8287-b344832a122f%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[David West]

Sorry for being unclear. It is the correlation within the variable that is measured with error. I have tried to include an example below.

#data

d1 <- cbind.data.frame(predictor = c(0.005125100, 0.005435463, 0.021340467, 0.006122980, 0.008509872, 0.009936922),

                       predictor_se = c(0.002682132, 0.003804434, 0.003118820, 0.002060798, 0.002657896, 0.002432118),

                       outcome = c(0.026238331, 0.009947423, 0.038082012, 0.006974313, 0.024211767, 0.002802170),

                       outcome_se = c(0.01271259, 0.01967391, 0.01713570, 0.01199105, 0.01356065, 0.01838512))

rho <- cbind(c(1.00000000, 0.11092781, 0.11570652, -0.21288003, 0.07870986, -0.09393294),

             c(0.01636902, 1.00000000, -0.07730337, 0.07658367, -0.04730668, 0.03229659),

             c(0.08243997, -0.06467834, 1.00000000, -0.28981654, 0.30530699, -0.13025649),

             c(-0.17657762, 0.08264548, -0.32989140, 1.00000000, -0.15878164, 0.46864858),

             c(0.06135205, -0.08437681, 0.36126178, -0.12203443, 1.00000000, -0.27311942),

             c(-0.09783581, -0.02560804, -0.15168974, 0.45968208, -0.30779561, 1.00000000))

#with brms (not accounting for correlation)

library(brms)

fit1 <- brm(outcome | se(outcome_se, sigma = TRUE) ~ 0 + me(predictor, predictor_se),

            data = d1, iter = 1000, chains = 2, cores = 1)

summary(fit1)

#non-bayesian inverse-variance weighted estimate (ignoring correlation)

beta1 <- summary(lm(outcome ~ 0 + predictor, weights = outcome_se^-2, data = d1))$coef[1]

#non-bayesian inverse-variance weighted estimate (accounting for correlation)

Sigma <- d1$outcome_se %o% d1$outcome_se * rho

X1 <- solve(t(chol(Sigma))) %*% d1$predictor

Y1 <- solve(t(chol(Sigma))) %*% d1$outcome

beta2 <- lm(Y1 ~ 0 + X1)$coef[1]

#plot

library(ggplot2)

ggplot(as.data.frame(fit1), aes(bme_mepredictorpredictor_se)) +

  geom_histogram() +

  geom_vline(xintercept = beta1, col = "red") +   # ignoring correlation

  geom_vline(xintercept = beta2, col = "green")   # accounting for correlation

Is it possible to account for the correlation matrix (rho) using brms?

Thanks!

Op vrijdag 13 april 2018 11:26:02 UTC+2 schreef Paul Buerkner:

[Paul Buerkner]

unfortunately this is not possible.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/cdf36a9e-ffb3-41ae-8e28-3d631db0e3de%40googlegroups.com.

[David West]

OK. Thanks. If you think that this is a useful feature for brms I can create an issue.

Op vrijdag 13 april 2018 16:35:34 UTC+2 schreef Paul Buerkner:

[Paul Buerkner]

not yet sure. feel free to open an issue anyway.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/3c2e6d2d-7948-4180-a759-eca9dc0e6df5%40googlegroups.com.