###########################
# MS-C1620
# Statistical inference
# Lecture 7


# We use the significance level 0.05 throughout the script



##########################################################
# Examples of data sets with near linear relationships


#####################
# Iris data

# The sepal lengths and widths for Setosa-irises
# An almost-linear relationship
plot(iris[1:50, 1], iris[1:50, 2], xlab = "Sepal length", ylab = "Sepal width")





####################
# Prestige data - linear with outliers

library(car)

data("Prestige")
Prestige
help(Prestige)

# What is the relationship between "education" and "income" (in 1971)?
plot(Prestige$education, Prestige$income, xlab = "Education", ylab = "Income")
# Mostly linear with some outliers in the top right corner




#####################
# Alligator data - linearity after a transform


# A study in central Florida where 15 alligators were captured.
# Two measurements were made on each of the alligators:
# The weight (in pounds),
# The snout vent length (in inches - this is the distance between the back of the head to the end of the nose).
# Source: https://www.r-bloggers.com/simple-linear-regression-2/
alligator = data.frame(
  Length = exp(c(3.87, 3.61, 4.33, 3.43, 3.81, 3.83, 3.46, 3.76,
               3.50, 3.58, 4.19, 3.78, 3.71, 3.73, 3.78)),
  Weight = exp(c(4.87, 3.93, 6.46, 3.33, 4.38, 4.70, 3.50, 4.50,
               3.58, 3.64, 5.90, 4.43, 4.38, 4.42, 4.25))
)

# Is the relationship linear?
plot(alligator$Length, alligator$Weight, xlab = "Length", ylab = "Weight")

# Maybe more so for logarithmic variables
plot(log(alligator$Length), log(alligator$Weight), xlab = "log(Length)", ylab = "log(Weight)")












##########################################################
# Fitting a linear model


# We choose the line so that deviations of the points from it are as small a possible
x <- rnorm(100)
y <- -1*x + rnorm(100, 0.1)
plot(x, y)

# A possible best line?
plot(x, y)
abline(a = -1, b = -0.5)
segments(x, y, x, -1 -0.5*x, col="red")
# Sum of squared deviations
sum((y - (-1 -0.5*x))^2)


# Another candidate?
plot(x, y)
abline(a = 1, b = -0.7)
segments(x, y, x, 1 -0.7*x, col="red")
# Sum of squared deviations
sum((y - (1 -0.7*x))^2)


# The optimal least squares choice:
plot(x, y)
abline(lm(y ~ x))
segments(x, y, x, fitted(lm(y ~ x)), col="red")
sum((y - fitted(lm(y ~ x)))^2)




###########
# Examples

# 1
x <- rnorm(100)
y <- x
plot(x, y)
abline(lm(y ~ x))


# 2
x <- rnorm(50)
y <- -1*x + rnorm(50, 0.1)
plot(x, y)
abline(lm(y ~ x))
lm(y ~ x)$coef

# Computing the coefficients also "by hand"

mx <- mean(x)
my <- mean(y)
sx <- sd(x)
sy <- sd(y)
corxy <- cor(x, y) 

beta_1 <- corxy*sy/sx
beta_1
beta_0 <- my - beta_1*mx
beta_0

# True line in red
abline(a = 0, b = -1, col = "red")




# 3
x <- rnorm(100)
y <- x^2 + x + rnorm(100, 0.1)
plot(x, y)
abline(lm(y ~ x))


# 4
phi <- runif(400, 0, 2*pi)
r <- runif(400, 1, 1.5)
x <- r*cos(phi)
y <- r*sin(phi)
plot(x, y)
abline(lm(y ~ x))


# Iris
plot(iris[1:50, 1], iris[1:50, 2], xlab = "Sepal length", ylab = "Sepal width")
abline(lm(iris[1:50, 2] ~ iris[1:50, 1]))
lm(iris[1:50, 2] ~ iris[1:50, 1])$coef
# 1cm increase in sepal length increases the sepal width by 0.80cm


# Prestige
plot(Prestige$education, Prestige$income, xlab = "Education", ylab = "Income")
abline(lm(Prestige$income ~ Prestige$education))
lm(Prestige$income ~ Prestige$education)$coef
# 1 point increase in education increases income by ~900$


# Alligator
plot(alligator$Length, alligator$Weight, xlab = "Length", ylab = "Weight")
abline(lm(alligator$Weight ~ alligator$Length))
lm(alligator$Weight ~ alligator$Length)$coef

plot(log(alligator$Length), log(alligator$Weight), xlab = "log(Length)", ylab = "log(Weight)")
abline(lm(log(alligator$Weight) ~ log(alligator$Length)))
lm(log(alligator$Weight) ~ log(alligator$Length))$coef



##########################################################
# Assessing model fit


##############################
# Residuals in iris data
x <- iris[1:50, 1]
y <- iris[1:50, 2]
plot(x, y, xlab = "Sepal length", ylab = "Sepal width")
lm_iris <- lm(iris[1:50, 2] ~ iris[1:50, 1])
abline(lm_iris)

res_iris <- signif(residuals(lm_iris), 5)
pred_iris <- predict(lm_iris) # plot distances between points and the regression line
segments(x, y, x, pred_iris, col="red")

# add labels (res values) to points
library(calibrate)
textxy(x, y, res_iris, cex=0.8)





##############################
# Coefficients of determination for the fits


# 1
x <- rnorm(100)
y <- x
plot(x, y)
abline(lm(y ~ x))
summary(lm(y ~ x))$r.squared


# 2
x <- rnorm(100)
y <- -1*x + rnorm(100, 0.1)
plot(x, y)
abline(lm(y ~ x))
summary(lm(y ~ x))$r.squared


# 3
x <- rnorm(100)
y <- x^2 + x + rnorm(100, 0.1)
plot(x, y)
abline(lm(y ~ x))
summary(lm(y ~ x))$r.squared


# 4
phi <- runif(400, 0, 2*pi)
r <- runif(400, 1, 1.5)
x <- r*cos(phi)
y <- r*sin(phi)
plot(x, y)
abline(lm(y ~ x))
summary(lm(y ~ x))$r.squared


# Iris
plot(iris[1:50, 1], iris[1:50, 2], xlab = "Sepal length", ylab = "Sepal width")
abline(lm(iris[1:50, 2] ~ iris[1:50, 1]))
lm(iris[1:50, 2] ~ iris[1:50, 1])$coef
summary(lm(iris[1:50, 2] ~ iris[1:50, 1]))$r.squared


# Prestige
plot(Prestige$education, Prestige$income, xlab = "Education", ylab = "Income")
abline(lm(Prestige$income ~ Prestige$education))
lm(Prestige$income ~ Prestige$education)$coef
summary(lm(Prestige$income ~ Prestige$education))$r.squared


# Alligator
plot(alligator$Length, alligator$Weight, xlab = "Length", ylab = "Weight")
abline(lm(alligator$Weight ~ alligator$Length))
lm(alligator$Weight ~ alligator$Length)$coef
summary(lm(alligator$Weight ~ alligator$Length))$r.squared
# Large even though line is not the "correct" choice!


plot(log(alligator$Length), log(alligator$Weight), xlab = "log(Length)", ylab = "log(Weight)")
abline(lm(log(alligator$Weight) ~ log(alligator$Length)))
lm(log(alligator$Weight) ~ log(alligator$Length))$coef
summary(lm(log(alligator$Weight) ~ log(alligator$Length)))$r.squared




#########################################################
# Inference for linear models


# Iris
iris_short <- iris[1:50, 1:2]
iris_lm <- lm(Sepal.Width ~ Sepal.Length, data = iris_short)
summary(iris_lm)
# The final column in the "Coefficients" table tells the p-value of the hypothesis tests of the form
# H0: beta = 0
# H1: beta != 0

# Small p-value for Sepal length -> its true value is not zero
# -> Sepal length has an effect on Sepal length (there is a true relationship between the two)

# Confidence intervals:
confint(iris_lm)

# How the interval is calculated:
iris_coef <- summary(iris_lm)$coefficients
iris_coef["Sepal.Length", "Estimate"] - qt(0.975, 48)*iris_coef["Sepal.Length", "Std. Error"]
iris_coef["Sepal.Length", "Estimate"] + qt(0.975, 48)*iris_coef["Sepal.Length", "Std. Error"]




# Prestige
prestige_short <- Prestige[, 1:2]
prestige_lm <- lm(income ~ education, data = prestige_short)
summary(prestige_lm)
# Education and income have a statistically significant relationship between them

# Confidence intervals:
confint(prestige_lm)



# Alligator
alligator_lm <- lm(log(alligator$Weight) ~ log(alligator$Length), data = alligator)
summary(alligator_lm)
# Length and weight (their logarithms) have a statistically significant relationship between them

# Confidence intervals:
confint(alligator_lm)



# Circle data
phi <- runif(400, 0, 2*pi)
r <- runif(400, 1, 1.5)
x <- r*cos(phi)
y <- r*sin(phi)
circle_lm <- lm(y ~ x)
summary(circle_lm)
# No *linear* relationship between x and y!

confint(circle_lm)