# Remember to set the wd 
# Some useful functions are 
# getwd() # to see the folder where you are
# setwd() # to set the wd by hand (otherwise right click on the R file)


# Read the data and explore the structure
Advertising <- read.csv("Advertising.csv")

dim(Advertising)
names(Advertising)
str(Advertising)

# (Minimal) Exploratory Analysis
summary(Advertising[,-1])

par(mfrow=c(1,4))
for (i in 2:5) hist(Advertising[,i], breaks = 30, main = names(Advertising)[i],
                    xlab = names(Advertising)[i])

cor(Advertising[,-1])

pairs(Advertising[,-1], cex = 0.5)

# Regressing Sales on TV (comments)
fitTV <- lm(Sales ~ TV, data = Advertising)
summary(fitTV)$coefficients


with(Advertising, plot(TV, Sales, cex = 0.5))
abline(fitTV$coefficients, col = "red")
points(Advertising$TV, predict(fitTV), col = "blue", cex = 0.5)

par(mfrow=c(1,2))
with(Advertising, plot(TV, Sales, cex = 0.5)); abline(fitTV$coefficients, col = "red")
points(Advertising$TV, predict(fitTV), col = "blue", cex = 0.5)
with(Advertising, segments(x0 = TV, y0 = Sales, x1 = TV, y1 = predict(fitTV), lty = "dashed"))
plot(fitTV, which = 1)


confint(fitTV)

summary(fitTV)$coefficients[2,]
2 * pt(-17.66763, df = 198)


summary(fitTV)$sigma

summary(fitTV)$r.squared

# Regress sales on Radio 
fitRadio <- lm(Sales ~ Radio, data = Advertising)
summary(fitRadio)$coefficients

# Regress sales on Newspper
fitNewspaper <- lm(Sales ~ Newspaper, data = Advertising)
summary(fitNewspaper)$coefficients


round(cor(Advertising[,-1]),3)

# Fit the model including all the predictors
fitLM <- lm(Sales ~ TV + Radio + Newspaper, data = Advertising)
round(summary(fitLM)$coefficients,4)
summary(fitLM)$sigma
summary(fitLM)$r.squared


summary(fitLM)
anova(lm(Sales ~ 1, data = Advertising), fitLM)

# Comparing R square
summary(fitTV)$r.squared
summary(fitRadio)$r.squared
summary(lm(Sales ~ TV + Radio, data = Advertising))$r.squared
summary(fitLM)$r.squared # Model with all the variables

# Comparing residual standard error
summary(fitTV)$sigma
summary(fitRadio)$sigma
summary(lm(Sales ~ TV + Radio, data = Advertising))$sigma
summary(fitLM)$sigma # Model with all the variables

# Analyze single p-values
round(summary(fitLM)$coefficients,3)

# Confidence intervals 
confint(fitLM)
 
# Take a look to the coefficents table  
summary(fitTV)$coefficients
summary(fitRadio)$coefficients
summary(fitNewspaper)$coefficients


# Accuracy of the prediction 
plot(Advertising$Sales, predict(fitLM))
abline(0,1)

# Confidence interval 
predict(fitLM, se.fit = TRUE, interval = "confidence",
        newdata = data.frame(TV = 100, Radio = 20, Newspaper = 0))$fit

# Prediction interval 
predict(fitLM, se.fit = TRUE, interval = "predict",
        newdata = data.frame(TV = 100, Radio = 20, Newspaper = 0))$fit


par(mfrow=c(2,2))
plot(fitLM)


# Intercation/Sinergy
# Let's consider the model including only TV and Radio 
fitTVRadio <- lm(Sales ~ TV + Radio, data = Advertising)
tv.seq <- seq(min(Advertising$TV), max(Advertising$TV), length.out = 30)
radio.seq <- seq(min(Advertising$Radio), max(Advertising$Radio), length.out = 30)
grid <- expand.grid(TV = tv.seq, Radio = radio.seq)

grid$Sales <- predict(fitTVRadio, newdata = grid)
Z <- matrix(grid$Sales, nrow = length(tv.seq), ncol = length(radio.seq))
par(mfrow=c(1,2))
p <- persp(tv.seq, radio.seq, Z, 
           theta = 45, phi = 15, expand = 0.6,
           col = "lightblue", shade = 0.1,
           xlab = "TV", ylab = "Radio", zlab = "Sales",
           ticktype = "detailed",
           main = "Sales ~ TV + Radio")

pts <- trans3d(Advertising$TV, Advertising$Radio, Advertising$Sales, p)
predicted <- predict(fitTVRadio)
visible <- Advertising$Sales > predicted
points(pts, col = ifelse(visible, "red", "blue"), pch = 19)

p <- persp(tv.seq, radio.seq, Z,
           theta = 135, phi = 15, expand = 0.6,
           col = "lightblue", shade = 0.1,
           xlab = "TV", ylab = "Radio", zlab = "Sales",
           ticktype = "detailed",
           main = "Sales ~ TV + Radio")

pts <- trans3d(Advertising$TV, Advertising$Radio, Advertising$Sales, p)
predicted <- predict(fitTVRadio)
visible <- Advertising$Sales > predicted
points(pts, col = ifelse(visible, "red", "blue"), pch = 19)

# Let's fit the model including the interaction term 
fitLMint <- lm(Sales ~ TV*Radio, data = Advertising)
# Equivalently
fitLMint <- lm(Sales ~ TV + Radio + TV:Radio, data = Advertising)
round(summary(fitLMint)$coefficients,4)
summary(fitLMint)$r.squared

par(mfrow=c(2,2))
plot(fitLMint)

# What if we transform the outcome
fitLMtrans <- lm(I(log(Sales)) ~ TV + I(TV^2) +  Radio + TV:Radio, data = Advertising)
summary(fitLMtrans)

par(mfrow=c(2,2))
plot(fitLMtrans)

par(mfrow = c(1,2))
plot(log(Advertising$Sales), predict(fitLMtrans))
abline(0,1)
plot(Advertising$Sales, exp(predict(fitLMtrans)))
abline(0,1)