# =============================================
# Advertising data
# =============================================
#
# Summary:
# These are records of sales and advertising budgets
# from 200 different markets
#
# Data description:
# The data set consists of the following 4 variables:
# Response variable:
#  Sales: sales in thousands of units
# Covariates:
#  TV: budget for advertising on TV in thousands of dollars
#  Radio: budget for advertising on radio in thousands of dollars
#  Newspaper: budget for advertising in newspapers in thousands of dollars
#
#======================================================

# Read the data
advertising <- read.csv("http://www.uio.no/studier/emner/matnat/math/STK2100/data/advertising.csv")

# Fit a model with all 3 covariates
fit.lm <- lm(Sales~TV+Radio+Newspaper,data=advertising)
summary(fit.lm)

# T-test for the effect of Newspaper
N <- nrow(advertising)
p <- 3
beta.hat <- coef(fit.lm)[4]
sigma.hat <- sqrt(sum(fit.lm$residuals^2)/(N-p-1))
sigma.hat
summary(fit.lm)$sigma
X <- model.matrix(fit.lm)
v <- solve(t(X)%*%X)
se.beta.hat <- sigma.hat*sqrt(v[4,4])
se.beta.hat
summary(fit.lm)$coef[4,2]
tobs <- beta.hat/se.beta.hat #observed t-statistic
tobs
df <- N-p-1
qt(0.975,df) #rejection level
2*(1-pt(abs(tobs),df)) #P-value

# F-test for the effects of all 3 covariates
fit.null <- lm(Sales~1,data=advertising) #fit null model
RSS <- sum(fit.lm$residuals^2)
RSS0 <- sum(fit.null$residuals^2)
q <- 3
Fobs <- ((RSS0-RSS)/q)/(RSS/(N-p-1))
Fobs
qf(0.95,q,N-p-1) #rejection level
1-pf(Fobs,q,N-p-1) #P-value

anova(fit.null,fit.lm)

# F-test for the effect of Newspaper
fit.lm2 <- lm(Sales~TV+Radio,data=advertising)
RSS <- sum(fit.lm$residuals^2)
RSS0 <- sum(fit.lm2$residuals^2)
q <- 1
Fobs <- ((RSS0-RSS)/q)/(RSS/(N-p-1))
Fobs
qf(0.95,q,N-p-1) #rejection level
1-pf(Fobs,q,N-p-1) #P-value

anova(fit.lm2,fit.lm)

# F-test for the effect of Newspaper and Radio
fit.lm3 <- lm(Sales~TV,data=advertising)
RSS <- sum(fit.lm$residuals^2)
RSS0 <- sum(fit.lm3$residuals^2)
q <- 2
Fobs <- ((RSS0-RSS)/q)/(RSS/(N-p-1))
Fobs
qf(0.95,q,N-p-1) #rejection level
1-pf(Fobs,q,N-p-1) #P-value

anova(fit.lm3,fit.lm)

#Confidence interval for the betas
beta.hat+se.beta.hat*qt(c(0.025,0.975),N-p-1)

confint(fit.lm,level=0.95)

#Confidence interval for E[Y|x] for a new x
newdata <- data.frame(TV=100,Radio=20,Newspaper=50)
y.hat <- predict(fit.lm,newdata)
y.hat
x <- c(1,100,20,50)
se.y.hat <- as.numeric(sigma.hat*sqrt(t(x)%*%v%*%x))
y.hat+se.y.hat*qt(c(0.025,0.975),N-p-1)
predict(fit.lm,newdata,interval="confidence")

#Prediction interval for Y for a new x
y.hat+sqrt(sigma.hat^2+se.y.hat^2)*qt(c(0.025,0.975),N-p-1)
predict(fit.lm,newdata,interval="predict")

# =============================================
# Credit data
# =============================================
#
# Summary:
# These are records of the average credit debt and
# for 400 customers, as well as diverse information
# about the customers, such as age, gender, credit
# rating, and so on
#
# Data description:
# The data set consists of the following 10 variables
# (the first colums contains the number of the individual):
# Response variable:
#  Balance: average credit debt
# Covariates:
#  Income: income in thousands of dollars
#  Limit: credit limit
#  Rating: credit rating
#  Cards: number of credit cards
#  Age: age in years
#  Education: years of education
#  Gender: gender
#  Student: is a student (Yes/No)
#  Married: is married (Yes/No)
#  Ethnicity: ethnicity (Caucasian/African American/Asian)
#
#======================================================

# Read the data
credit <- read.csv("http://www.uio.no/studier/emner/matnat/math/STK2100/data/credit.csv")

# Box plots of Balance for the qualitative covariates Gender, Student, Status and Ethnicity
par(mfrow=c(2,2))
boxplot(Balance~factor(Gender),data=credit)
boxplot(Balance~factor(Student),data=credit)
boxplot(Balance~factor(Married),data=credit)
boxplot(Balance~factor(Ethnicity),data=credit)
par(mfrow=c(1,1))

# Fit individual models with Income and each of the qualitative
# covariates Student and Ethnicity
fit.income <- lm(Balance~Income,data=credit)
summary(fit.income)
fit.student <- lm(Balance~Income+factor(Student),data=credit)
summary(fit.student)
plot(credit$Income,credit$Balance)
abline(fit.income$coef[1],fit.income$coef[2],col='green',lwd=2)
abline(fit.student$coef[1],fit.student$coef[2],col='red',lwd=2)
abline(fit.student$coef[1]+fit.student$coef[3],fit.income$coef[2],col='blue',lwd=2)

fit.ethnicity <- lm(Balance~Income+factor(Ethnicity),data=credit)
summary(fit.ethnicity)
plot(credit$Income,credit$Balance)
abline(fit.income$coef[1],fit.income$coef[2],col='green',lwd=2)
abline(fit.ethnicity$coef[1],fit.ethnicity$coef[2],col='red',lwd=2)
abline(fit.ethnicity$coef[1]+fit.ethnicity$coef[3],fit.ethnicity$coef[2],col='blue',lwd=2)
abline(fit.ethnicity$coef[1]+fit.ethnicity$coef[4],fit.ethnicity$coef[2],col='purple',lwd=2)

# Zoomm in to see the differences
plot(credit$Income,credit$Balance,xlim=c(40,50),ylim=c(488,550))
abline(fit.income$coef[1],fit.income$coef[2],col='green',lwd=2)
abline(fit.ethnicity$coef[1],fit.ethnicity$coef[2],col='red',lwd=2)
abline(fit.ethnicity$coef[1]+fit.ethnicity$coef[3],fit.ethnicity$coef[2],col='blue',lwd=2)
abline(fit.ethnicity$coef[1]+fit.ethnicity$coef[4],fit.ethnicity$coef[2],col='purple',lwd=2)