# January 22, 2024

library(tidyverse)

d <- read.csv("https://whlevine.hosted.uark.edu/psyc5143/ch7.csv", header = T)

# scatterplot with linear relation visualized
ggplot(data = d, aes(x = MILES, y = TIME)) +
	geom_point() +
	geom_smooth(method = "lm", se = F)

usualModel <- lm(TIME ~ MILES, d)
summary(usualModel)

# easy residuals plot
plot(usualModel, which = 1)

# scatterplot with relation visualized
ggplot(data = d, aes(x = MILES, y = TIME)) +
	geom_point() +
	geom_smooth(se = F)

# mean-centering miles
d <- d %>% 
	mutate(MILES.c = MILES - mean(MILES))

# adding miles-squared to the data
d <- d %>% 
	mutate(M2 = MILES.c^2)

polyModel <- lm(TIME ~ MILES.c + M2, d)
summary(polyModel)

# some stuff for graphing purposes
coef(polyModel)[2] -> slope0 # the slope at mean miles (MILES.c = 0)
coef(polyModel)[1] + coef(polyModel)[2]*(0 - mean(d$MILES)) -> int0 # the y-intercept of the tangent line at mean miles


# scatterplot with some extra info added
plot1 <- ggplot(data = d, aes(x = MILES, y = TIME)) +
	geom_point() +
	geom_vline(xintercept = mean(d$MILES)) + # mean MILES
	stat_smooth(method = "lm", formula = y ~ x + I(x^2), se = F) # the model, visualized

# the slope at MILES.c = 0 (i.e., mean miles)
plot2 <- plot1 + geom_abline(slope = slope0, intercept = int0, col = "red")
plot2

# what's the slope at MILES = 40?
d <- d %>% 
	mutate(M40 = MILES - 40,
				 M40sq = M40^2)

coef(lm(TIME ~ M40 + M40sq, d))

polyModel40 <- lm(TIME ~ M40 + M40sq, d)

# visualizing the slope at MILES = 40

slope40 <- coef(polyModel40)[2]
intercept40 <- coef(polyModel40)[1] + coef(polyModel40)[2]*(-40)
plot2 + geom_abline(slope = slope40, intercept = intercept40, col = "black")