Rethinking Measurement Error - models using the rethinking package that incorporate measurement in both the response and predictor variables.

rethinking_measurement_errors.R

library(rethinking)
library(dplyr)
data(WaffleDivorce)


mod <- alist(
  #likelihood
  div_est ~ dnorm(mu, sigma),
  
  #model
  mu <- a + bA*age_obs + bR*mar_est[i],
  
  #error in both response and predictor
  div_obs ~ dnorm(div_est, div_sd),
  mar_obs ~ dnorm(mar_est, mar_sd),
  
  #priors
  a ~ dnorm(0,10),
  bA ~ dnorm(0,10),
  bR ~ dnorm(0,10),
  sigma ~ dcauchy(0, 2.5)
)

#do a little manipulation of the data frame
#for fitting
d <- WaffleDivorce %>%
  dplyr::select(Marriage,
         Marriage.SE,
         MedianAgeMarriage,
         Divorce,
         Divorce.SE) %>%
  rename(age_obs = MedianAgeMarriage,
         mar_obs = Marriage,
         mar_sd = Marriage.SE,
         div_obs = Divorce,
         div_sd = Divorce.SE)


#Fit using STAN
fit <- map2stan(mod,
                data = d,
                start = list(div_est = d$div_obs,
                             mar_est = d$mar_obs),
                             WAIC=FALSE,
                             iter=5000, 
                             warmup = 1000,
                             chains = 3,
                             cores = 3,
                             control = list(adapt_delta=0.95))

plot(fit, pars="bR", window=c(1000, 5000))
plot(fit, pars="bA", window=c(1000, 5000))


#Can we do a model with measurement error in both predictors and response?
set.seed(2016)
d <- d %>%
  mutate(age_sd = abs(rnorm(n(), sd(age_obs), sd = 0.5))) #abs to make sure sd's positive

mod_all_err <- alist(
  #likelihood
  div_est ~ dnorm(mu, sigma),
  
  #model
  mu <- a + bA*age_est[i] + bR*mar_est[i],
  
  #error in both response and predictor
  age_obs ~ dnorm(age_est, age_sd),
  div_obs ~ dnorm(div_est, div_sd),
  mar_obs ~ dnorm(mar_est, mar_sd),
  
  #priors
  a ~ dnorm(0,10),
  bA ~ dnorm(0,10),
  bR ~ dnorm(0,10),
  sigma ~ dcauchy(0, 2.5)
)


fit_all_err <- map2stan(mod_all_err,
                data = d,
                start = list(div_est = d$div_obs,
                             mar_est = d$mar_obs,
                             age_est = d$age_obs),
                WAIC=FALSE,
                iter=5000, 
                warmup = 1000,
                chains = 3,
                cores = 3,
                control = list(adapt_delta=0.95))

plot(fit_all_err, pars="bR", window=c(1000, 5000))
plot(fit_all_err, pars="bA", window=c(1000, 5000))

#Can we do a model with measurement error in both predictors and response?
#With an interaction?

mod_int <- alist(
  #likelihood
  div_est ~ dnorm(mu, sigma),
  
  #model
  mu <- a + bA*age_est[i] + bR*mar_est[i] + bI*age_est[i]*mar_est[i],
  
  #error in both response and predictor
  age_obs ~ dnorm(age_est, age_sd),
  div_obs ~ dnorm(div_est, div_sd),
  mar_obs ~ dnorm(mar_est, mar_sd),
  
  #priors
  a ~ dnorm(0,10),
  bA ~ dnorm(0,10),
  bR ~ dnorm(0,10),
  bI ~ dnorm(0,10),
  sigma ~ dcauchy(0, 2.5)
)



fit_all_int <- map2stan(mod_int,
                        data = d,
                        start = list(div_est = d$div_obs,
                                     mar_est = d$mar_obs,
                                     age_est = d$age_obs),
                        WAIC=FALSE,
                        iter=5000, 
                        warmup = 1000,
                        chains = 3,
                        cores = 3,
                        control = list(adapt_delta=0.95))


plot(fit_all_int, pars="bR", window=c(1000, 5000))
plot(fit_all_int, pars="bA", window=c(1000, 5000))
plot(fit_all_int, pars="bI", window=c(1000, 5000))
plot(fit_all_int, pars="sigma", window=c(1000, 5000))
plot(fit_all_int, window=c(1000, 5000))

precis(fit_all_int)