Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Plot vertical density of normal distribution in R and ggplot2

I want to plot vertical normal distributions in ggplot2 into an existing linear regression model, in order to visualize homoscedasticity.
I currently have the following code and plot:

x <- runif(100, 0, 15)
y <- 1000 + 200*x + rnorm(100, 0, 300)
df <- data.frame(x, y)

lm_fit <- lm(y ~ x, data = df)


#with regression line
ggplot(df, mapping = aes(x=x, y=y)) + geom_point(color="blue") + geom_smooth(method='lm', se=FALSE, color="red") 

sample regression

And I would like to insert density curves like this (just facing the other way round):
sample wanted result

like image 287
Doflaminhgo Avatar asked Mar 02 '23 11:03

Doflaminhgo


1 Answers

library(ggplot2)

x <- runif(100, 0, 15)
y <- 1000 + 200*x + rnorm(100, 0, 300)
df <- data.frame(x, y)
lm_fit <- lm(y ~ x, data = df)

k <- 2.5
sigma <- sigma(lm_fit)
ab <- coef(lm_fit); a <- ab[1]; b <- ab[2]

x <- seq(-k*sigma, k*sigma, length.out = 50)
y <- dnorm(x, 0, sigma)/dnorm(0, 0, sigma) * 3

x0 <- 0
y0 <- a+b*x0
path1 <- data.frame(x = y + x0, y = x + y0)
segment1 <- data.frame(x = x0, y = y0 - k*sigma, xend = x0, yend = y0 + k*sigma)
x0 <- 5
y0 <- a+b*x0
path2 <- data.frame(x = y + x0, y = x + y0)
segment2 <- data.frame(x = x0, y = y0 - k*sigma, xend = x0, yend = y0 + k*sigma)
x0 <- 10
y0 <- a+b*x0
path3 <- data.frame(x = y + x0, y = x + y0)
segment3 <- data.frame(x = x0, y = y0 - k*sigma, xend = x0, yend = y0 + k*sigma)

ggplot(df, mapping = aes(x=x, y=y)) + geom_point(color="blue") + 
  geom_smooth(method='lm', se=FALSE, color="red") + 
  geom_path(aes(x,y), data = path1, color = "green") + 
  geom_segment(aes(x=x,y=y,xend=xend,yend=yend), data = segment1) +
  geom_path(aes(x,y), data = path2, color = "green") + 
  geom_segment(aes(x=x,y=y,xend=xend,yend=yend), data = segment2) +
  geom_path(aes(x,y), data = path3, color = "green") + 
  geom_segment(aes(x=x,y=y,xend=xend,yend=yend), data = segment3)

enter image description here

like image 64
Stéphane Laurent Avatar answered Mar 05 '23 16:03

Stéphane Laurent