library(ggplot2)
pearson <- read.csv("https://turtlegraphics.org/stat2300/data/pearson-father-son.csv")
xbar = mean(pearson$father)
sx = sd(pearson$father)
ybar = mean(pearson$son)
sy = sd(pearson$son)

gp <- ggplot(data=pearson, aes(x=father,y=son)) + 
  geom_hline(yintercept = c(ybar+sy,ybar,ybar-sy),color="blue") +
  geom_vline(xintercept = c(xbar+sx,xbar,xbar-sx),color="blue") +
  geom_abline(slope=sy/sx, intercept=ybar-(sy/sx)*xbar,color="green") +
  geom_point(size=0.5) +
  geom_smooth(method="lm", se=FALSE, color="red")
print(gp)