# Linear regression with prior

Im trying to estimate the linear curve (y~x) where I know intercept must be normally distributed around -100, and slope always positive and normally distributed around 2 (blue continous line in plots below).

The example.data.1 below is "clean" and the linear regression (red dashed line) is ok. The resulting red dashed line is what I want.

But example.data.2 has many measurement errors so the red dashed line becomes unrealistic. The resulting line should be parallel with the blue line, but lower.

How can I assign a strong prior similar to the blue line in the plots, so that I get a posterior reasonably similar to the blue line?

example.data.1 <- structure(list(x = c(1.36, 2.22, 2.53, 3.09, 3.44, 3.25, 3.15,
3.21, 3.57, 3.63, 3.51, 2.85, 2.56, 2.25, 1.61, 1.35, 1, 1.6,
1.92, 1.9, 2.3, 2.61, 3.9, 3.74, 3.77, 3.77, 3.49, 3.37, 3.35,
2.79, 2.31, 1.88, 1.5, 1.18, 1.83, 2.32, 3.06, 3.37, 3.77, 3.82,
3.75, 3.72, 3.53, 3.35, 3.67, 3.18, 3.11, 2.43, 1.9, 1.39, 1.17,
1.48, 2.05, 2.62, 3.08, 3.65, 3.92, 4.08, 4.1, 3.47, 3.84, 3.45,
2.87, 2.83, 2.49, 1.87, 2.06, 2.49, 1.78, 2.33, 2.95, 3.73, 3.64,
3.62, 4.1, 3.85, 4.06, 3.67, 3.3, 2.86, 2.46, 2.32, 2.08, 1.64,
1.96), y = c(-101.04, -99.42, -98.33, -96.88, -95.22, -91.89,
-91.63, -90.19, -92.98, -95.58, -95.69, -96.32, -96.94, -98.25,
-100.11, -100.81, -101.87, -99.72, -99.94, -100.87, -100.38,
-98.64, -93.38, -92.98, -93.39, -93.76, -93.25, -93.12, -94.46,
-96.45, -97.46, -99.75, -100.09, -101.62, -101.1, -97.8, -96.33,
-96.21, -94.37, -93.18, -93.32, -93.73, -94.13, -94.4, -94.63,
-94.83, -96.29, -98.11, -100.2, -100.82, -101.56, -101.35, -100.61,
-98.65, -97.37, -95.36, -95.45, -95.33, -95.63, -95.26, -97.08,
-97.1, -97.14, -96.9, -98.17, -99.47, -100.17, -100.58, -100.55,
-99.94, -99.02, -97.3, -96.25, -95.44, -95.69, -95.21, -95.87,
-95.87, -97.71, -96.91, -97.62, -97.94, -98.9, -99.79, -99.88
)), .Names = c("x", "y"), row.names = c(NA, -85L), class = "data.frame")

example.data.2 <- structure(list(x = c(3.11, 3.46, 3.42, 3.34, 3.3, 2.45, 4, 4.2,
4.08, 3.57, 1.97, 1.83, 1.07, 0.68, 0.54, 0.47, 0.63, 3.19, 3.52,
3.49, 3.47, 3.36, 2.76, 3.42, 3.17, 3.54, 2.56, 1.06, 1.09, 0.84,
0.64, 0.61, 0.74, 0.49, 3.49, 3.56, 3.46, 3.25, 3.72, 3.57, 3.58,
2.62, 1.99, 1.85, 1.04, 1.06, 0.62, 0.49, 0.48, 0.68, 0.5, 3.63,
3.71, 3.75, 3.67, 3.78, 3.52, 3.04, 2.26, 1, 1.17, 1.01, 0.92,
0.65, 0.54, 0.36, 0.38, 0.3, 3.08, 3.79, 3.9, 3.5, 3.4, 2.57,
3.03, 1.93, 2.02, 1.5, 0.67, 0.63, 0.72, 0.6, 0.67, 0.63, 0.53
), y = c(-105.28, -104.1, -104.81, -104.34, -104.37, -105.31,
-103.59, -103.32, -102.66, -103.57, -103.73, -104.47, -97.69,
-92.56, -95.9, -95.72, -107.6, -104.39, -105.12, -104.18, -104.46,
-102.19, -103.59, -103.38, -103.48, -102.84, -96.52, -88.54,
-90.36, -93.7, -85.21, -89.68, -99.47, -91.92, -104.58, -103.91,
-104.47, -104.49, -104.41, -104.41, -102.6, -98.65, -87.98, -89.23,
-86.34, -94.21, -91.57, -84.62, -84.14, -95.33, -102.14, -104.18,
-103.8, -102.47, -101.75, -101.73, -102.84, -97.49, -92.67, -91.72,
-80.45, -80.97, -84.94, -80.2, -81.05, -77.84, -82.72, -91.75,
-105.19, -104.66, -104.36, -104.31, -103.57, -102.68, -98.4,
-89.48, -85.92, -84.59, -84.49, -81.13, -83.28, -83.12, -85.62,
-85.89, -90.07)), .Names = c("x", "y"), row.names = c(NA, -85L
), class = "data.frame")

lm.1 <- coefficients(lm(example.data.1$y~example.data.1$x))
lm.2 <- coefficients(lm(example.data.2$y~example.data.2$x))

library(ggplot2)
p <- ggplot(example.data.1, aes(x=x, y=y))
p <- p + geom_point()
p <- p + geom_abline(intercept=(-100), slope=2, color="blue")
p <- p + geom_abline(intercept=lm.1[1], slope=lm.1[2], color="red", linetype="dashed")
p <- p + xlim(0, 10)
p <- p + ylim(-110, -50)
p

p <- ggplot(example.data.2, aes(x=x, y=y))
p <- p + geom_point()
p <- p + geom_abline(intercept=(-100), slope=2, color="blue")
p <- p + geom_abline(intercept=lm.2[1], slope=lm.2[2], color="red", linetype="dashed")
p <- p + xlim(0, 10)
p <- p + ylim(-110, -50)
p


I need to do this for tens of thousands of data-groups like each example above, so a fast algorithm is important. I have tried to use Stan, Jags, and arm-package, but don't understand how to tell those functions what I want.

My limited knowledge about statistics lead me to think that a bayesian approach is best, but I could be wrong.

• Did you forget the plots? – user30490 Oct 9 '15 at 2:28
• They are in the code, not sure how to paste them into the question. Sorry about that. – Chris Oct 9 '15 at 2:30
• Are you certain that what you are trying to do is the best solution to the problem? The bottom line is, that under the assumptions of linear regression, your data are telling a very different story than the theory. Using a very strong prior because you don't like the results isn't logically consistent. Rather, if you believe the problem is caused by measurement error, you should apply an appropriate technique to correct for it, such as total least squares, instrumental variables, etc. That way, your modeling assumptions match the situation your modeling. – Zachary Blumenfeld Oct 9 '15 at 5:54
• Another option could be to assume that the error term in the regression model is not normal, but e.g. a t distribution with 3 DFs (or a mixture of a normal and something long-tailed). The long-tailed part of the error terms would then reflect the potential for measurement errors. – Björn Oct 9 '15 at 6:13
• @Bjorn You actually have to be careful with that one, the presence of measurement error still violates the assumptions of such model(s) . Moreover, doing what you recommend could further compound the problem as the introduction of non-linearities and distributions outside of common location-scale families will likely further bias the regression when measurement error exists, even if such distributions seem to provide a better predictive fit. My feeling is that it would be better to attend to the measurement error phenomenon first before playing around with more sophisticated distributions. – Zachary Blumenfeld Oct 9 '15 at 9:54