Measurement Error
Measurement Error
Measurement error is the difference between a measured quantity and its true value. It can be due to
- systematic bias (e.g., a scale is miscalibrated by 1 pound for everyone)
- random error (e.g., some people take off their shoes, others are wearing coats, some may be dehydrated or have just eaten) that may be naturally occurring and may occur with any experiment
Measurement error is often countered by tactics like taking the mean of multiple measurements (e.g., research quality blood pressure measures take the mean of three values) or standardizing experimental conditions. However, sometimes substantial sources of error are unavoidable.
Example: Divorce and Marriage Rate
McElreath (2016) considers the relationship among divorce rate, marriage rate, and median age at marriage based on state-level data.
Note that data from Nevada are not included.
Is divorce associated with marriage? Well, yes!
However, does a high marriage rate imply a high divorce rate? How does median age at marriage affect divorce rates?
One issue analyzing these data is that we have error involved in the measurement of both marriage rate and divorce rate. First, we’ll explore measurement error of our outcome, divorce rate.
library(rethinking)
data(WaffleDivorce)
d=WaffleDivorce
plot(d$Divorce~d$MedianAgeMarriage,ylim=c(4,15),
xlab="Median age at marriage",ylab="Divorce rate per 1000 population")
#add interval of 1 SE in each direction
for (i in 1:nrow(d)) {
ci=d$Divorce[i]+c(-1,1)*d$Divorce.SE[i]
x=d$MedianAgeMarriage[i]
lines(c(x,x),ci)
}
There is substantial variability in the certainty in the estimated divorce rates – why?
A hunch is that the size of the state’s population may be involved.
plot(d$Divorce~log(d$Population),ylim=c(4,15),
xlab="Log(population)",ylab="Divorce rate per 1000 population")
#add interval of 1 SE in each direction
for (i in 1:nrow(d)) {
ci=d$Divorce[i]+c(-1,1)*d$Divorce.SE[i]
x=log(d$Population[i])
lines(c(x,x),ci)
}
Yes, there is a relationship between population size and certainty in the estimated rate!
We also see this in marriage rates!
Handling Measurement Error
First, we focus on measurement error in our response, the divorce rate. One reasonable approach is to use a hierarchical model.
- Define the parameter \(D_{TRUE,i}\) to be the true (unknown) divorce rate for state \(i\)
- Define our observed outcome (subject to measurement error) as \(D_{OBS,i}\) and its associated standard error (provided in the data) as \(D_{SE,i}\)
- Model \(D_{OBS,i} \sim N\left(D_{TRUE,i},D_{SE,i}^2\right)\)
- Here the observed divorce rates are centered on the true rates with the estimated measurement error treated as known
- Define as covariates the median age at marriage \(A_i\) and marriage rate \(R_i\)
Now we can specify our desired model, for the true divorce rates, as follows.
\[D_{TRUE,i} \sim N(\mu_i,\sigma^2)\]
\[\mu_i=\beta_0+\beta_1A_i+\beta_2R_i\]
\[D_{OBS,i} \sim N\left(D_{TRUE,i},D_{SE,i}^2\right)\]
\[\beta_j \sim N(0,100)\] \[\sigma \sim \text{HalfCauchy}(0,2.5)\]
First, we fit the model with no adjustment for measurement error, so that the outcome is just the observed (with error) divorce rate.
library(brms)
#put data into a list
dlist <- list(
div_obs = d$Divorce,
div_sd = d$Divorce.SE,
R = d$Marriage,
A = d$MedianAgeMarriage-mean(d$MedianAgeMarriage))
m1 <-
brm(data = dlist, family = gaussian,
div_obs ~ 0+ Intercept + R+A,
#had to rescale intercept prior
prior = c(prior(normal(0,50),class=b,coef=Intercept),
prior(normal(0, 10), class = b),
prior(cauchy(0, 2.5), class = sigma)),
iter = 5000, warmup = 1000, chains = 4, cores = 4,
seed = 14,control=list(adapt_delta=0.95))m1## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: div_obs ~ 0 + Intercept + R + A
## Data: dlist (Number of observations: 50)
## Samples: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 10.80 1.67 7.52 14.08 1.00 4600 5663
## R -0.06 0.08 -0.22 0.11 1.00 4611 5433
## A -0.99 0.25 -1.49 -0.50 1.00 5407 6482
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.52 0.16 1.24 1.87 1.00 6590 6625
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The interpretation of this model is that while marriage rate is not associated with divorce rate conditional on median age at marriage, conditional on the marriage rate, a one-year higher median age at marriage is associated with an expected 0.99 fewer divorces per 1000 population (95% CrI=(0.50,1.49)). However, we may be concerned because of the error in determination of our outcome.
What does \(D_{OBS,i} \sim N\left(D_{TRUE,i},D_{SE,i}^2\right)\) imply for each state?
d %>%
mutate(Divorce_distribution = str_c("Divorce ~ Normal(", Divorce, ", ", Divorce.SE, ")")) %>%
select(Loc, Divorce_distribution) %>%
head()## Loc Divorce_distribution
## 1 AL Divorce ~ Normal(12.7, 0.79)
## 2 AK Divorce ~ Normal(12.5, 2.05)
## 3 AZ Divorce ~ Normal(10.8, 0.74)
## 4 AR Divorce ~ Normal(13.5, 1.22)
## 5 CA Divorce ~ Normal(8, 0.24)
## 6 CO Divorce ~ Normal(11.6, 0.94)You can see big differences in the error for CA (very large state) and AK (small state in terms of population).
# here we specify the initial (i.e., starting) values
inits <- list(Yl = dlist$div_obs)
inits_list <- list(inits, inits)
#want 0+intercept notation if data not mean-centered
m2 <-
brm(data = dlist, family = gaussian,
div_obs | mi(div_sd) ~ 0 + intercept + R + A,
prior = c(prior(normal(0, 10), class = b),
prior(cauchy(0, 2.5), class = sigma)),
iter = 5000, warmup = 1000, cores = 2, chains = 2,
seed = 14,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
save_mevars = TRUE, # note this line for the `mi()` model
inits = inits_list)m2## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: div_obs | mi(div_sd) ~ 0 + intercept + R + A
## Data: dlist (Number of observations: 50)
## Samples: 2 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 8000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept 9.29 1.74 5.88 12.71 1.00 2688 4380
## R 0.01 0.09 -0.16 0.19 1.00 2700 4341
## A -0.97 0.25 -1.47 -0.47 1.00 3296 4859
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.07 0.20 0.72 1.48 1.00 2554 3868
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The interpretation of this model is similar to what we saw before, though our estimate of \(\sigma\) is now lower.
Measurement Error in Exposure
Measurement error in the exposure variable, here marriage rate, can have an effect on estimation as well. Here we allow the marriage rate to be measured with error as well by fitting the following model.
\[D_{TRUE,i} \sim N(\mu_i,\sigma^2)\]
\[\mu_i=\beta_0+\beta_1A_i+\beta_2R_{TRUE,i}\]
\[D_{OBS,i} \sim N\left(D_{TRUE,i},D_{SE,i}^2\right)\]
\[R_{OBS,i} \sim N\left(R_{TRUE,i},R_{SE,i}^2\right)\]
\[\beta_j \sim N(0,100)\] \[\sigma \sim \text{HalfCauchy}(0,2.5)\]
What does \(R_{OBS,i} \sim N\left(R_{TRUE,i},R_{SE,i}^2\right)\) imply for each state?
d %>%
mutate(Marriage_distribution = str_c("Marriage ~ Normal(", Marriage, ", ", Marriage.SE, ")")) %>%
select(Loc, Marriage_distribution) %>%
head()## Loc Marriage_distribution
## 1 AL Marriage ~ Normal(20.2, 1.27)
## 2 AK Marriage ~ Normal(26, 2.93)
## 3 AZ Marriage ~ Normal(20.3, 0.98)
## 4 AR Marriage ~ Normal(26.4, 1.7)
## 5 CA Marriage ~ Normal(19.1, 0.39)
## 6 CO Marriage ~ Normal(23.5, 1.24)dlist <- list(
div_obs = d$Divorce,
div_sd = d$Divorce.SE,
mar_obs = d$Marriage,
mar_sd = d$Marriage.SE,
A = d$MedianAgeMarriage)
# the `inits`
inits <- list(Yl = dlist$div_obs)
inits_list <- list(inits, inits)
# the model
m3 <-
brm(data = dlist, family = gaussian,
div_obs | mi(div_sd) ~ 0 + intercept + me(mar_obs, mar_sd) + A,
prior = c(prior(normal(0, 10), class = b),
prior(cauchy(0, 2.5), class = sigma)),
iter = 5000, warmup = 1000, cores = 2, chains = 2,
seed = 1235,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
save_mevars = TRUE,
inits = inits_list)## [[1]]
## Stan model 'cf10f52a9e78d534fbc2d2a42e429046' does not contain samples.## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: div_obs | mi(div_sd) ~ 0 + intercept + me(mar_obs, mar_sd) + A
## Data: dlist (Number of observations: 50)
## Samples: 1 chains, each with iter = 5000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## intercept 15.88 6.92 1.89 29.17 1.00 1911
## A -0.44 0.21 -0.84 -0.03 1.00 2113
## memar_obsmar_sd 0.27 0.11 0.06 0.49 1.00 1743
## Tail_ESS
## intercept 2280
## A 2390
## memar_obsmar_sd 2587
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.01 0.21 0.63 1.49 1.00 924 1587
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Now that we’ve accounted for measurement error in the exposure and outcome, we see substantial changes in effect estimates. The interpretation of this model is that conditional on the marriage rate, a one-year higher median age at marriage is associated with an expected 0.44 fewer divorces per 1000 population (95% CrI=(0.04,0.83)). Conditional on the median age at marriage, an increase of the marriage rate by 1 per 1000 is associated with an expected increase in the divorce rate of 0.27 per 1000 (95% CrI=(0.07, 0.49)).
Moral of the Story
The moral of this story is that when you have error associated with a predictor or response (i.e., a distribution of responses), reducing the response to a single value – discarding uncertainty – can lead to spurious inference.