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

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)

model.year <- lm(Admission~Month, data=fullmoon)
summary(model.year)

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

points(tapply(fullmoon$Admission, fullmoon$Month, mean) ~ c(1:12), pch=8)

# (b)
# ---
model.year.num <- lm(Admission~as.numeric(Month), data=fullmoon)
summary(model.year.num)
abline(coefficients(model.year.num))

anova(model.year.num, model.year)

# (c)
# ---
model.year.num.poly <- lm(Admission~poly(as.numeric(Month), 11), data=fullmoon)
summary(model.year.num.poly)
lines(fitted(model.year.num.poly)[1:12]~as.numeric(fullmoon$Month)[1:12], col=3)
summary(model.year)

model.year.num.poly.5 <- lm(Admission~poly(as.numeric(Month), 5), data=fullmoon)
summary(model.year.num.poly.5)
lines(fitted(model.year.num.poly.5)[1:12]~as.numeric(fullmoon$Month)[1:12], col=4)

anova(model.year.num.poly.5, model.year.num.poly)
anova(model.year.num, model.year.num.poly)

# (d)
# ---
model.year <- lm(Admission~Month, data=fullmoon, contrasts=list(Month=contr.poly))
summary(model.year)
summary(model.year.num.poly)

round(solve(cbind(1, contr.poly(11))), 3)

1-pf((summary(model.year.num.poly)$coefficients[, 3])^2, df1=1, df2=24)
1-pf((summary(model.year.num.poly)$coefficients[, 3])^2/11, df1=11, df2=24)
0.05/11

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

model.year <- lm(Admission~Month, data=fullmoon)

A <- matrix(0, nrow=66, ncol=12)

# August vs others
a <- c(-1, 1, rep(0, 10)) # mu
a <- c(0, 1, rep(0, 10)) # beta

for(i in 1:11){A[i, i+1] <- 1}

# September vs others
a <- c(0, -1, 1, rep(0, 9)) # mu
a <- c(0, -1, 1, rep(0, 9)) # beta

j <- 3
for(i in 12:21){
  A[i, 2] <- -1
  A[i, j] <- 1
  j <- j+1
}

# ...

j <- 4
for(i in 22:30){
  A[i, 3] <- -1
  A[i, j] <- 1
  j <- j+1
}

j <- 5
for(i in 31:38){
  A[i, 4] <- -1
  A[i, j] <- 1
  j <- j+1
}

r.index <- 39
for(c.index in 6:12){
  j <- c.index
  for(i in r.index:(r.index + 12 - c.index)){
    A[i, c.index-1] <- -1
    A[i, j] <- 1
    j <- j+1
  }
  r.index <- r.index + 12 - c.index + 1
}

A

ci.all.t <- data.frame(Lower=rep(NA, 66), Estimate=NA, Upper=NA, StdError=NA, Tval=NA, Pval=NA, Sig=NA)
ci.all.t$Estimate <- A%*%coef(model.year)
ci.all.t$StdError <- sqrt(diag(A%*%vcov(model.year)%*%t(A)))
ci.all.t$Lower <- ci.all.t$Estimate-qt(0.975, df=24)*ci.all.t$StdError
ci.all.t$Upper <- ci.all.t$Estimate+qt(0.975, df=24)*ci.all.t$StdError
ci.all.t$Tval <- ci.all.t$Estimate/ci.all.t$StdError
ci.all.t$Pval <- 2*pt(-abs(ci.all.t$Tval), df=24)
ci.all.t$Sig <- ifelse(ci.all.t$Pval<0.05, "*", "")
ci.all.t

ci.all.t[1:11, ]
summary(model.year)

# (a)
# ---
A[1:11, 1] <- -1

Estimate.LSD <- A%*%tapply(fullmoon$Admission, fullmoon$Month, mean)
StdError.LSD <- summary(model.year)$sigma*sqrt(2/3)

sum(abs(Estimate.LSD-ci.all.t$Estimate))
sum(abs(StdError.LSD-ci.all.t$StdError))

# (b)
# ---
ci.all.S <- data.frame(Lower=NA, Estimate=ci.all.t$Estimate, Upper=NA, StdError=ci.all.t$StdError, Fval=(ci.all.t$Tval)^2/11, Pval=NA, Sig=NA)
ci.all.S$Lower <- ci.all.S$Estimate-sqrt(qf(0.95, df1=11, df2=24)*11)*ci.all.S$StdError
ci.all.S$Upper <- ci.all.S$Estimate+sqrt(qf(0.95, df1=11, df2=24)*11)*ci.all.S$StdError
ci.all.S$Pval <- 1-pf(ci.all.S$Fval, df1=11, df2=24)
ci.all.S$Sig <- ifelse(ci.all.S$Pval<0.05, "*", "")
ci.all.S

# (c)
# ---
ci.all.HSD <- TukeyHSD(aov(Admission~Month, data=fullmoon))
ci.all.HSD
?TukeyHSD

?qtukey

# comparison
# ----------
plot(ci.all.HSD)
for (i in 1:66){
  segments(x0=ci.all.S$Lower[i], x1=ci.all.S$Upper[i], y0=67-i, col="deeppink")
  segments(x0=ci.all.HSD$Month[i, 2], x1=ci.all.HSD$Month[i, 3], y0=67-i, col="firebrick")
  segments(x0=ci.all.t$Lower[i], x1=ci.all.t$Upper[i], y0=67-i, col="orange")
}

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

# (a)
# ---
ci.all.B <- ci.all.t
ci.all.B$Lower <- ci.all.t$Estimate-qt(1-0.05/2/66, df=24)*ci.all.t$StdError
ci.all.B$Upper <- ci.all.t$Estimate+qt(1-0.05/2/66, df=24)*ci.all.t$StdError
ci.all.B$Pval <- ci.all.t$Pval*66
ci.all.B$Sig <- ifelse(ci.all.B$Pval<0.05, "*", "")
ci.all.B

?pairwise.t.test
pairwise.t.test(x=fullmoon$Admission, g=fullmoon$Month, p.adjust.method="bonferroni", pool.sd=TRUE)
round(ci.all.B$Pval[c(7:9, 11, 18, 27, 35, 42, 45)], 5)