This is all super-simple; still, it might be useful. In class today a student asked for some intuition as to why, when you’re regressing y on x, measurement error on x biases the coefficient estimate but measurement error on y does not.
I gave the following quick explanation:
– You’re already starting with the model, y_i = a + bx_i + e_i. If you add measurement error to y, call it y*_i = y_i + eta_i, and then you regress y* on x, you can write y* = a + bx_i + e_i + eta_i, and as long as eta is independent of e, you can just combine them into a single error term.
– When you have measurement error in x, two things happen to attenuate b—that is, to pull the regression coefficient toward zero. First, if you spreading out x but keep y unchanged, this will reduce the slope of y on x. Second, when you add noise to x you’re changing the ordering of the data, which will reduce the strength of the relationship.
But that’s all words (and some math). It’s simpler and clearer to do a live simulation, which I did right then and there in class!
Here’s the R code:
# simulation for measurement error
library("arm")
set.seed(123)
n <- 1000
x <- runif(n, 0, 10)
a <- 0.2
b <- 0.3
sigma <- 0.5
y <- rnorm(n, a + b*x, sigma)
fake <- data.frame(x,y)
fit_1 <- lm(y ~ x, data=fake)
display(fit_1)
sigma_y <- 1
fake$y_star <- rnorm(n, fake$y, sigma_y)
sigma_x <- 4
fake$x_star <- rnorm(n, fake$x, sigma_x)
fit_2 <- lm(y_star ~ x, data=fake)
display(fit_2)
fit_3 <- lm(y ~ x_star, data=fake)
display(fit_3)
fit_4 <- lm(y_star ~ x_star, data=fake)
display(fit_4)
x_range <- range(fake$x, fake$x_star)
y_range <- range(fake$y, fake$y_star)
par(mfrow=c(2,2), mar=c(3,3,1,1), mgp=c(1.5,.5,0), tck=-.01)
plot(fake$x, fake$y, xlim=x_range, ylim=y_range, bty="l", pch=20, cex=.5)
abline(coef(fit_1), col="red", main="No measurement error")
plot(fake$x, fake$y_star, xlim=x_range, ylim=y_range, bty="l", pch=20, cex=.5)
abline(coef(fit_2), col="red", main="Measurement error on y")
plot(fake$x_star, fake$y, xlim=x_range, ylim=y_range, bty="l", pch=20, cex=.5)
abline(coef(fit_3), col="red", main="Measurement error on x")
plot(fake$x_star, fake$y_star, xlim=x_range, ylim=y_range, bty="l", pch=20, cex=.5)
abline(coef(fit_4), col="red", main="Measurement error on x and y")
The resulting plot is at the top of this post.
I like this simulation for three reasons:
1. You can look at the graph and see how the slope changes with measurement error in x but not in y.
2. This exercise shows the benefits of clear graphics, including little things like making the dots small, adding the regression lines in red, labeling the individual plots, and using a common axis range for all four graphs.
3. It was fast! I did it live in class, and this is an example of how students, or anyone, can answer this sort of statistical question directly, with a lot more confidence and understanding than would come from a textbook and some formulas.
P.S. As Eric Loken and I discuss in this 2017 article, everything gets more complicated if you condition on "statistical significance."
P.P.S. Yes, I know my R code is ugly. Think of this as an inspiration: even if, like me, you’re a sloppy coder, you can still code up these examples for teaching and learning.
I don’t understand this sentence: “In class today a student asked for some intuition as to why, when you’re regressing y on x, measurement error on x biases the coefficient estimate by measurement error on y does not.” Should “by” be “but”?
Should “estimate by measurement” read “estimate but measurement . . .”?
I don’t think that is such ugly code. I can mostly understand it and I don’t use R.
Bob76
Typo fixed; thanks.
P.S. The code is fine for my purposes, but I think it’s ugly to users of ggplot2 and tidyverse in R.
I’ll bet you didn’t mean to include a 12.5 mB image in your post! I’m surprised, but it clobbered the Feedly app on my iPad.
A useful heuristic in many situations is to take things to the extreme. It’s easy to see that measurement error in x when exploded to infinite standard deviation, will lead to a line with zero slope. I remember a discussion on this on the blog a year ago ish and I didn’t get the point of the question until Carlos made me look harder at it, I finally used this heuristic of blowing up the x errors to very large size and it immediately made sense.
Worth it to remember this trick in many applications.
Nice simulation by Andrew, and cool idea Daniel. I will create a shiny app for class where one can dial up the uncertainty on x and on y. This will help students a lot to understand why regressions on averaged values need to be seen with some skepticism (psycholinguistics is littered with such plots, ignoring measurement error but reporting significant correlations, even in top-notch, best-in-category Nature-ish type journals).
The simulation code is very readable. I like that about R’s being purpose-built for stats. Python code for this gets cluttered up with namespace issues (a convention that would help R avoid the inevitable package variable/function naming conflicts).
I think you can remove “display()” when writing about plotting code. I looked it up, and it’s just a replacement for the default lm print in the arm package that prints fewer statistics at lower arithmetic precision.
The plotting code is opaque to non R users because there’s a function called “plot()”—if this were called “scatter_plot(x = , y = )” it’d be easier to understand what it would do without looking at the output. For people like me who don’t use base graphics in R much, I have no idea what all that cut-and-paste args are doing, like “bty=”l”, pch=20, cex=.5″ and “mar=c(3,3,1,1), mgp=c(1.5,.5,0), tck=-.01”. The way to clean this up is to write a function encapsulating all the redundancy and binding the repeated arguments.
x_range <- ... y_range <- ... lr_plot <- function(fit, x, y, fit) { plot(x, y, xlim=x_range, ylim=y_range, bty="l", pch=20, cex=.5) abline(coef(fit_1), col="red", main=title) } lr_plot(fit, fake$x, fake$y, "No measurement error") lr_plot(fit_err_y, fake$x, fake$y_star, "Measurement error on y") lr_plot(fit_err_x, fake$x_star, fake$y, "Measurement error on x") lr_plot(fit_err_xy, fake$x_star, fake$y_star, "Measurement error on x and y")This makes it clear the same thing is happening 4 times and it makes it easeir to see where the plots vary. It also makes it easier to check that the right titles are being attached to the right graphs. I went further and renamed the fits to make it easier to see the alignment. Even better would be to pull the data out of the fit object---I'm pretty sure you can do that because of the everything-and-the-kitchen-sink design of fit objects. You could go even further and name those terse args if you want to. Something like "tiny_dot_plotting_character <- 20" and then "pch=tiny_dot_plotting_character" (OK, I did know this one) or "tick_adjustment <- -0.1" and "tck=tick_adjustment" (in the previous code)---when you name variables this way, you no longer need doc.
I'd be curious what best practices in the Tidyverse would look like for this. I suspect that either there are built ins for lm, or you'd have to wrestle the data into a dataframe, set the faceting, extract the coefficients for the lines, etc. I find the ggplot code pretty readable (more so than base graphics), but I find all the data frame manipulation and piping very opaque.
Here’s one way:
library(tidyverse)
set.seed(123)
n <- 1000
a <- 0.2
b <- 0.3
sigma <- 0.5
fake %
mutate(y_star = rnorm(n, y, sigma_y),
x_star = rnorm(n, x, sigma_x))
bind_rows(
tibble(x=fake$x, y=fake$y, name=”No measurement error”),
tibble(x=fake$x, y=fake$y_star, name=”Measurement error on y”),
tibble(x=fake$x_star, y=fake$y, name=”Measurement error on x”),
tibble(x=fake$x_star, y=fake$y_star, name=”Measurement error on x and y”)
) %>%
mutate(name = fct_inorder(name)) %>%
ggplot(aes(x,y)) +
geom_point() +
geom_smooth(method=”lm”, fullrange=TRUE) +
facet_wrap(~name)
Couple of notes: geom_smooth() has a (imo) nice default of not extrapolating the fitted line beyond the range of the data, hence the fullrange option to match the original. By default the plots would have appeared in alphabetical order by title, hence the slightly esoteric fct_inorder().
An alternative to faceting that might be more readable for non-ggplot users is to make individual plots and compose them using the (excellent!) patchwork.
Bob:
The call to display() takes up no more space than the call to print(), and I prefer the display() output, so that’s what I use. If you use print(), you don’t get standard errors for the coefficient estimates. If you use summary(), you get all sorts of extra stuff like t-statistics and p-values. I wrote display() to give me just what I want.
> The plotting code is opaque to non R users because there’s a function called “plot()”—if this were called “scatter_plot(x = , y = )” it’d be easier to understand what it would do without looking at the output.
It doesn’t seem that difficult to guess what the following may do:
plot(fake$x, fake$y, xlim=x_range, ylim=y_range, …)
Anyway it doesn’t matter much – the objective of the code is not to make easy to understand what may happen when you run it but to make that happen when you run it.
Interestingly you called your auxiliary function lr_plot (which is invalid because of a repeated formal argument by the way) and not scatter_plot or lr_scatter_plot even though that would make it easier to understand what it does without looking at the output.
Any code done live in front of students gets an automatic aesthetic pass. Plus, I agree with Bob76, this seems totally fine for what it is.
Dmitri:
It’s partly a rhetorical trick, that I badmouth my own code and then the commenters can reassure me how great it is. On the other hand, sometimes I badmouth my own code and the commenters agree!
Yeah, I know that one, I use it on my family whenever I can. “You are going to hate this suggestion, but …” It leverages the listener’s inherent contrarianism in your favor. Like, they’re stuck between doing what you predict and doing what you want!
I also use this trick in writing papers. I put in something that any contrarian (i.e., 100% of academics) would object to but is not fatal to the paper. Then the reviewer inevitably objects, and feels happy about themselves and all is good. I call it a honeypot; it helps manage the psychological state of the reviewer, it assuages their strong desire to shoot down the paper just because it’s there.
This trick also works on accreditation self-studies. Put in some obvious things that are “deficient” such as inadequate measurement of learning objectives – the accreditation team feels like they did their job by pointing out the need for improvement. If you give a few things like that, then they might not look any deeper. It is games like these that bore me with accreditation and publishing papers.
Reminds me of what R.A. Heinlein had a writer character say in “Stranger In a Strange Land”: “Always put in something for the editor to take out. After they piss in it, they like the flavor better.”
Thanks for posting this, I knew and trusted this to be true. However fleshing it out like that is a killer demo IMHO.
Jean:
Yes, indeed, there’s nothing like a fake-data demonstration. I’ve been doing more and more of these. Sometimes they can be tricky to set up. I’d like to write an article, something about how to construct a fake-data demonstration and how it can be difficult to do.
Horizontal clusters of systematic errors due to DIE policy of no child left behind auto grade advancement at EA12 and EA16.
https://twitter.com/dux_ie/status/1775706592811585688
@Andrew, do you teach Deming regression in the case of measurement error in both X & Y?
Nice post. What would be really cool would be an easy way to make *animated* versions in R. It would start with the no-errors data and then slowly show the points as you added noise in the x dimension, with the regression line slowly flattening. Then a second animation would return to the original data and add noise to the y values, with the regression line remaining rock-solid in its original position. Maybe there already is an R package that enables such animated graphics, but not to my knowledge.
You can do something like this:
x = (1:100)/10
y = x+rnorm(100, 0, 5)
noisex = rnorm(100, 0, 10)
noisey = rnorm(100, 0, 10)
plotWithNoise = function(kx=0, ky=0){
x <- x+kx*noisex
y <- y+ky*noisey
plot(y~x, xlim=c(-15,25), ylim=c(-15, 25), las=1, bty="n", xlab="", ylab="")
abline(lm(y~x))
}
animation::saveGIF({
for (kx in (0:50)/50) plotWithNoise(kx=kx)
for (kx in (50:0)/50) plotWithNoise(kx=kx)
for (ky in (0:50)/50) plotWithNoise(ky=ky)
for (ky in (50:0)/50) plotWithNoise(ky=ky)
},
movie.name = "animation.gif",
interval=0.1
)
Cool. Impressive display of R’s capabilities. It had struck me that such an animation would be easier in Mathematica than in R. But I doubt that it would be much simpler.
Bob76
Output: https://imgur.com/a/41rBwuH
Nice!! Animated GIFs are not just little toys.
Nice demo.
Regarding the code aesthetics – way too much white space ;)
This is half a nice example. It shows the effect of measurement on bias of the coefficients, but the effect on precision is largely implicit.
This can be made whole simply by adding data (concentration) ellipses to the plot.
library(car)
demo_plot <- function(x, y, fit, title) {
dataEllipse(x, y,
xlim=x_range, ylim=y_range,
pch=20, levels = 0.9,
main = title)
abline(coef(fit), col="red", lwd=2)
}
op <- par(mfrow=c(2,2),
mar=c(3,3,1,1),
mgp=c(1.5,.5,0), tck=-.01)
demo_plot(fake$x, fake$y, fit_1, "No measurement error")
demo_plot(fake$x, fake$y_star, fit_2, "Measurement error on y")
demo_plot(fake$x_star, fake$y, fit_3, "Measurement error on x")
demo_plot(fake$x_star, fake$y_star, fit_4, "Measurement error on x and y")
par(op)
Love this, Andrew! Such a simple yet elegant way to show how measurement error impacts regression models. Below is Stata code in case any Stata-centric folks run across this post and want to implement in their stats program of choice.
clear
version 16.1
set seed 5029
set obs 1000
gen x = 0+int((10-0+1)*runiform())
gen e = rnormal(0,.5)
gen y = .2 + .3*x + e
regress y x
graph twoway (scatter y x, xlabel(-10(5)20) ylabel(-2(2)6)) ///
(lfit y x, range(-10 20)), saving(no_msrmt_error) ///
legend(off) title(“No measurement error”) scheme(plottig)
*Add measerment error to y
gen y_star = y + rnormal(0,1)
regress y_star x
graph twoway (scatter y_star x, xlabel(-10(5)20) ylabel(-2(2)6)) ///
(lfit y_star x, range(-10 20)), saving(y_msrmt_error) ///
legend(off) title(“Measurement error in y”) scheme(plottig)
*Add measurement error to x
gen x_star = x + rnormal(0,4)
regress y x_star
graph twoway (scatter y x_star, xlabel(-10(5)20) ylabel(-2(2)6)) ///
(lfit y x_star, range(-10 20)), saving(x_msrmt_error) ///
legend(off) title(“Measurement error in x”) scheme(plottig)
*Regress error-prone y on x
regress y_star x_star
graph twoway (scatter y_star x_star, xlabel(-10(5)20) ylabel(-2(2)6)) ///
(lfit y_star x_star, range(-10 20)), saving(x_y_msrmt_error) ///
legend(off) title(“Measurement error in x and y”) scheme(plottig)
graph combine no_msrmt_error.gph y_msrmt_error.gph ///
x_msrmt_error.gph x_y_msrmt_error.gph, col(2) scheme(plottig)
I had also just re-created the example in Stata to use in my course next year. Complementary code ;)
set obs 1000
gen x = runiform(0,10)
gen y = 0.2 + 0.3*x + rnormal(0,0.5)
gen y_star = y + rnormal(0,1)
gen x_star = x + rnormal(0,4)
graph twoway (scatter y x) (lfit y x), ylabel(-2(2)6) xlabel(-10(5)20) legend(off) saving(plot1) title(“No measurement error”) ytitle(“y”)
graph twoway (scatter y_star x) (lfit y_star x), ylabel(-2(2)6) xlabel(-10(5)20) legend(off) saving(plot2) title(“Measurement error in y”) ytitle(“y_star”)
graph twoway (scatter y x_star) (lfit y x_star), ylabel(-2(2)6) xlabel(-10(5)20) legend(off) saving(plot3) title(“Measurement error in x”) ytitle(“y”)
graph twoway (scatter y_star x_star) (lfit y_star x_star), ylabel(-2(2)6) xlabel(-10(5)20) legend(off) saving(plot4) title(“Measurement error in both x and y”) ytitle(“y_star”)
graph combine plot1.gph plot2.gph plot3.gph plot4.gph
Adding my voice to the chorus exclaiming that this is fantastic and instructive post!
For those who prefer python, here is Andrew’s code translated:
“`
# import libraries
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
import seaborn as sns
sns.set_theme()
# create dataframe
np.random.seed(123)
n = 1000
x = scipy.stats.uniform.rvs(0, 10, n)
a = 0.2
b = 0.3
sigma = 0.5
y = scipy.stats.norm.rvs(a+b*x, sigma, n)
y_star = scipy.stats.norm.rvs(y, sigma + 0.5, n)
x_star = scipy.stats.norm.rvs(x, sigma + 3.5, n)
df = pd.DataFrame({“x”: x, “x_star”: x_star, “y”: y, “y_star”: y_star})
# plot
palette = sns.color_palette()
scatter_kws = {“alpha”: 0.3, “s”: 2}
line_kws = {“color”: palette[3]}
fig, axs = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(8, 8))
axs[0][0].set_title(“No measurement error”)
sns.regplot(data=df, x=”x”, y=”y”, line_kws=line_kws, scatter_kws=scatter_kws, ax=axs[0][0])
axs[0][1].set_title(“Measurement error in y”)
sns.regplot(data=df, x=”x”, y=”y_star”, line_kws=line_kws, scatter_kws=scatter_kws, ax=axs[0][1])
axs[0][0].set_title(“Measurement error in x”)
sns.regplot(data=df, x=”x_star”, y=”y”, line_kws=line_kws, scatter_kws=scatter_kws, ax=axs[1][0])
axs[0][1].set_title(“Measurement error in both x and y”)
sns.regplot(data=df, x=”x_star”, y=”y_star”, line_kws=line_kws, scatter_kws=scatter_kws, ax=axs[1][1]);
“`
One thing that I find a confusing about this example is that the measurement error added to y is quite a bit smaller than the error that is added to x. For the purpose of comparison, I made another version of Andrew’s plots where the same error is added to x and y. The overall message remains the same (of course): measurement error in x has a much greater effect on the slope than error in y. But, I think this alternate version brings out some of the nuance a bit better.
My code and the plots are available here:
https://github.com/ciyer/intro-data-viz/blob/master/sundry/Regression-Measurement-Error.ipynb
“Second, when you add noise to x you’re changing the ordering of the data, which will reduce the strength of the relationship.”
I’m guessing you mean you have the points (2,2) (3,3) and (4,4), adding noise can give you (2,2) (4,3) and (3,4), swapping the order, is that right? But doesn’t that also happen as you add noise to y? Is it countered by the increased spread in y and just leads to greater uncertainty?
Kj,
Yes, good point. I think I was trying too hard to come up with an intuitive explanation here.
This is very cool! One potential caveat I’d ask about is that in some disciplines (e.g., business, where I am a PhD student) they are using proxies for their y variables, such that they’re not directly “measuring” their y construct of interest. In such cases “measurement error” would seem to be more than an additive problem. And yet you often hear in econometrics classes that “measurement error of y doesn’t matter since it doesn’t bias the coefficients.” My intuition is that we can’t assume a constant additive measurement error when using proxies for y and therefore measurement error in y remains a concern, but I’m not a statistician/econometrician. Thoughts?
My intuition for this is that the slope coefficient is Cov(y,x)/Var(x). Cov(y,x) does not change by this type of measurement error, hence no bias if only y has measurement error (but se goes up because the variance of the error increases). Measurement error in x increase Var(x), hence bias towards zero. Showing this with simulations adds to this intuition.
While this is true if measurement error on y is addictive (y*_i = y_i + eta_i), this is not the case for multiplicative measurement error (y*_i = y_i * eta_i)
Chao:
Yes, if there is ever any question about this sort of thing, my recommended approach is to check it using simulation!