#------------------------------------------------------------------------------#
# Ex4.1 Two-way classification model                                           # 
#
# Last revised: Spring 2019
#
# Copyright © 2018 by Rodrigo Labouriau                                               #
#------------------------------------------------------------------------------#

# Basic preparation 

# Loading and attaching the required data file (please insert the address
# of the directory where the file referred above is placed in your computer)
load("DataFramesStatisticalModelling.RData")

ls()

str(Ex4.1)

attach(Ex4.1)

# Item a)

plot(week.day, n.cars, col=week, pch=19, cex=1.5,
     xlab="Day of the week", ylab="Number of cars")

plot(week.day+runif(length(week))/2, n.cars, col=week, pch=19, cex=1.1,
     xlab="Day of the week", ylab="Number of cars")


Means <- tapply(n.cars, week.day, mean)
Means

points(1:5, Means, type="b", cex=2, lwd=3, pch=19)

table(week, week.day)

week <- factor(week)

week.day <- factor(week.day)

Full <- glm(n.cars ~ week.day + week + week.day:week, family=poisson(link="log"))
summary(Full)

length(week.day)

length(coef(Full))

deviance(Full)

pchisq(deviance(Full), df=150, lower.tail=F)

Multiplicative <- glm(n.cars ~ week.day + week , 
                      family=poisson(link="log"))

anova(Multiplicative, Full, test="Chisq")

OnlyDay <- glm(n.cars ~ week.day , 
                      family=poisson(link="log"))

anova(OnlyDay, Multiplicative, test="Chisq")

Null <- glm(n.cars ~ 1 , 
              family=poisson(link="log"))

anova(Null, OnlyDay, test="Chisq")

summary(OnlyDay)

RawResiduals <- residuals(Full, "response")

Fitted <- fitted(Full)

cbind(n.cars, week, week.day, Fitted, RawResiduals)

plot(Fitted, RawResiduals)

PearsonResiduals <- residuals(Full, "pearson")

plot(Fitted, PearsonResiduals)

detach(Ex4.1)
################################################################################

# END !