#############
## Cviceni 1
#############

fullmoon <- read.table("fullmoon.txt", header=TRUE)

summary(fullmoon)

fullmoon$Moon <- factor(fullmoon$Moon, levels=c("Before", "During", "After"))

fullmoon$Month <- factor(fullmoon$Month, levels=c("Aug", "Sep", "Oct", "Nov", "Dec", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul"))

fullmoon$Fullmoon <- factor(ifelse(fullmoon$Moon=="During", "Fullmoon", "Before/After"), levels=c("Fullmoon", "Before/After"))


table(fullmoon$Fullmoon, fullmoon$Moon)

summary(fullmoon)

plot(fullmoon$Admission~fullmoon$Month)
plot(fullmoon$Admission~fullmoon$Moon)

plot(fullmoon$Admission~as.numeric(fullmoon$Month), ylab="Admissions", xlab="Month", main="Admissions by month and moon")
points(fullmoon$Admission[fullmoon$Moon=="Before"]~as.numeric(fullmoon$Month)[fullmoon$Moon=="Before"], pch=19, col="orange")
points(fullmoon$Admission[fullmoon$Moon=="During"]~as.numeric(fullmoon$Month)[fullmoon$Moon=="During"], pch=19, col="firebrick")
points(fullmoon$Admission[fullmoon$Moon=="After"]~as.numeric(fullmoon$Month)[fullmoon$Moon=="After"], pch=19, col="springgreen")
abline(h=mean(fullmoon$Admission), col=2, lty=2, lwd=2)

legend(x="topleft", legend=levels(fullmoon$Moon), col=c("orange", "firebrick", "springgreen"), pch=19)

#############
## Cviceni 2
#############

# (a)
# ---
model.month.and.moon <- lm(Admission ~ Month + Moon, data=fullmoon)
summary(model.month.and.moon)
anova(model.month.and.moon)

anova(lm(Admission ~ Moon + Month, data=fullmoon))

library(car)
?Anova
Anova(model.month.and.moon, type="II")
Anova(model.month.and.moon, type="III")

# (b)
# ---
model.month <- lm(Admission ~ Month, data=fullmoon)
model.moon <- lm(Admission ~ Moon, data=fullmoon)
anova(model.month, model.month.and.moon)
anova(model.moon, model.month.and.moon)

model.0 <- lm(Admission~1, data=fullmoon)
anova(model.0, model.moon)
anova(model.0, model.month)

# (c)
# ---
model.month.times.moon <- lm(Admission ~ Month * Moon, data=fullmoon)
summary(model.month.times.moon)

#############
## Cviceni 3
#############

# (a)
# ---
model.month.times.fullmoon <- lm(Admission ~ Month * Fullmoon, data=fullmoon)
summary(model.month.times.fullmoon)

anova(model.month.times.fullmoon)
Anova(model.month.times.fullmoon, type="II")
Anova(model.month.times.fullmoon, type="III")

points(fitted(model.month.times.fullmoon)[1:12]~c(1:12), pch=11)
lines(fitted(model.month.times.fullmoon)[1:12]~c(1:12), pch=11, lty=2)

points(fitted(model.month.times.fullmoon)[13:24]~c(1:12), pch=8)
lines(fitted(model.month.times.fullmoon)[13:24]~c(1:12), pch=11, lty=2)

# vsimneme si, ze fitovane hodnoty behem uplnku jsou rovny pozorovanym hodnotam 
# nas model je "prefitovany" a neni rozumne jej pozit
# ve zbytku cviceni budeme dale pracovat s timto modelem za ucelem ilustrace
# pro realnou analyzu bychom sloucili kategorie i v promenne Month (napr. uvazovat rocni obdobi),
# abychom dostali vic nez jedno pozorovani na kombinaci kategorii

# (b)
# ---
model.month.and.fullmoon <- lm(Admission ~ Month + Fullmoon, data=fullmoon)
model.month <- lm(Admission ~ Month, data=fullmoon)
model.fullmoon <- lm(Admission ~ Fullmoon, data=fullmoon)

anova(model.month.and.fullmoon, model.month.times.fullmoon)
anova(model.month, model.month.and.fullmoon)
# ...

# (c)
# ---
model.month.times.fullmoon <- lm(Admission ~ Fullmoon * Month, data=fullmoon)
summary(model.month.times.fullmoon)
anova(model.month.times.fullmoon)
Anova(model.month.times.fullmoon, type="II")
Anova(model.month.times.fullmoon, type="III")

model.month.times.fullmoon <- lm(Admission ~ Month * Fullmoon, data=fullmoon, contrasts = list(Month=contr.sum))
summary(model.month.times.fullmoon)
anova(model.month.times.fullmoon)
Anova(model.month.times.fullmoon, type="II")
Anova(model.month.times.fullmoon, type="III")

model.month.times.fullmoon <- lm(Admission ~ Month * Fullmoon, data=fullmoon)
summary(model.month.times.fullmoon)
anova(model.month.times.fullmoon)
Anova(model.month.times.fullmoon, type="II")
Anova(model.month.times.fullmoon, type="III")

# (d)
# ---
fullmoon.1 <- fullmoon[-1, ]

model.month.times.fullmoon.1 <- lm(Admission ~ Month * Fullmoon, data=fullmoon.1)
model.month.times.fullmoon.2 <- lm(Admission ~ Fullmoon * Month, data=fullmoon.1)

anova(model.month.times.fullmoon.1)
anova(model.month.times.fullmoon.2)
Anova(model.month.times.fullmoon.1, type="II")
Anova(model.month.times.fullmoon.1, type="III")

Anova(model.month.times.fullmoon.2, type="II")
Anova(model.month.times.fullmoon.2, type="III")

#############
## Cviceni 5
#############

# (a)
# ---
TukeyHSD(aov(Admission~Month*Fullmoon, data=fullmoon))

# (c)
# ---
?pairwise.t.test
monthfullmoon <- factor(paste(as.character(fullmoon$Month), as.character(fullmoon$Fullmoon)))
pairwise.t.test(x = fullmoon$Admission, g = monthfullmoon, p.adjust.method="bonferroni")
# nedava vysledky, protoze mame jenom jednom pozorovani na kategorii pro nektere kategorie 
# => nelze odhadnout rozptyl pro tyto skupiny
# pro rozumny model (dostatek pozorovani na kategorii) by tento kod daval zadane vysledky