#############
## 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("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))

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")

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

model.0 <- lm(Admission~1, data=fullmoon)
model.1 <- lm(Admission~Moon, data=fullmoon)

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

mu <- mean(fullmoon$Admission)
sst <- sum((fullmoon$Admission-mu)^2)
mu.before <- mean(fullmoon$Admission[fullmoon$Moon=="Before"])
mu.during <- mean(fullmoon$Admission[fullmoon$Moon=="During"])
mu.after <- mean(fullmoon$Admission[fullmoon$Moon=="After"])
ssa <- 12*(mu.before-mu)^2+12*(mu.during-mu)^2+12*(mu.after-mu)^2
sse <- sum((fullmoon$Admission[fullmoon$Moon=="Before"]-mu.before)^2)+sum((fullmoon$Admission[fullmoon$Moon=="During"]-mu.during)^2)+sum((fullmoon$Admission[fullmoon$Moon=="After"]-mu.after)^2)

mu
mu.before
mu.during
mu.after
ssa
sse
sst
ssa/2
sse/33
sst/35
sqrt(sse/33)
sqrt(sst/35)

summary(model.0)
mu
sqrt(sst/35)
t.test(fullmoon$Admission)

summary(model.1)
mu.before
mu.during-mu.before
mu.after-mu.before
sqrt(sse/33)
1-sse/sst
1-(sse/33)/(sst/35)
(ssa/2)/(sse/33)
1-pf((ssa/2)/(sse/33), df1=2, df2=33)

anova(model.1)
ssa
ssa/2
sse
sse/33
(ssa/2)/(sse/33)
1-pf((ssa/2)/(sse/33), df1=2, df2=33)

anova(model.0, model.1)
sst
sse 
ssa
(ssa/2)/(sse/33)
1-pf((ssa/2)/(sse/33), df1=2, df2=33)

aov(Admission~Moon, data=fullmoon)
summary(aov(Admission~Moon, data=fullmoon))
sse 
ssa
ssa/2
sse/33
(ssa/2)/(sse/33)
1-pf((ssa/2)/(sse/33), df1=2, df2=33)
?aov

oneway.test(Admission~Moon, data=fullmoon, var.equal=TRUE)
(ssa/2)/(sse/33)
1-pf((ssa/2)/(sse/33), df1=2, df2=33)
?oneway.test

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

A <- rbind(c(1, 0, 0), c(1, 1, 0), c(1, 0, 1))
tab.lm <- summary(model.1)$coefficients
tab.anova <- data.frame(row.names=c("Before", "During", "After"), Estimate=rep(NA, 3), StdError=NA, Tvalue=NA, Pvalue=NA, PvaluePos=NA, PvalueNeg=NA)
tab.anova[, 1] <- A%*%tab.lm[, 1]
tab.anova[, 2] <- sqrt(diag(A%*%vcov(model.1)%*%t(A)))
tab.anova[, 3] <- tab.anova[, 1]/tab.anova[, 2]
tab.anova[, 4] <- 2*pt(-abs(tab.anova[, 3]), df=summary(model.1)$df[2])
tab.anova[, 5] <- 1-pt(tab.anova[, 3], df=summary(model.1)$df[2])
tab.anova[, 6] <- pt(tab.anova[, 3], df=summary(model.1)$df[2])
tab.anova

xx <- seq(-4, 4, length=501)
yy <- dt(xx, df=33)
plot(yy~xx, type="l", main="Density of t(33)", xlab="x", ylab="f(x)")

xx <- seq(0, 5, length=501)
yy <- df(xx, df1=2, df2=33)
plot(yy~xx, type="l", main="Density of F(2,33)", xlab="x", ylab="f(x)")

###########
# Cviceni 4
###########

model.matrix(model.1)

model.2 <- lm(Admission~-1+Moon, data=fullmoon)
model.matrix(model.2)
summary(model.2)
tab.anova
summary(model.1)

model.help <- lm(Admission-mean(Admission)~-1+Moon, data=fullmoon)
model.matrix(model.help)
summary(model.help)
summary(model.2)
summary(model.1)

###########
# Cviceni 6
###########

summary(model.1)
model.matrix(model.1)
summary(model.1)
contr.treatment(3)

levels(fullmoon$Moon)

c1 <- c(1/3, 1/3, 1/3)
c2 <- c(-1, 1, 0)
c3 <- c(0, -1, 1)
C <- rbind(c1, c2, c3)
C.solve <- solve(C)
C.cont <- C.solve[, -1]
model.3 <- lm(Admission~Moon, data=fullmoon, contrasts=list(Moon=C.cont))
summary(model.3)
summary(model.1)
model.matrix(model.3)

mu
-mu.before+mu.during
-mu.during+mu.after

model.4 <- lm(Admission~relevel(Moon, ref="During"), data=fullmoon)
summary(model.4)
summary(model.1)
model.matrix(model.4)

model.5 <- lm(Admission~Moon, data=fullmoon, contrasts=list(Moon=contr.sum(3)))
summary(model.5)
summary(model.1)
model.matrix(model.5)

?contr.sum