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

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

C <- cbind(1, contr.sum(12))
round(solve(C), digits=3)
1/12
1-1/12
mean(fullmoon$Admission)
tapply(fullmoon$Admission, fullmoon$Month, mean)-mean(fullmoon$Admission)

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

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

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

# (a)
# ---
summary(model.year)

# (b)
# ---
a <- c(rep(-1/5, 5), rep(1/7, 7)) # pro mu
t(a)%*%tapply(fullmoon$Admission, fullmoon$Month, mean)

a <- c(0, rep(-1/5, 4), rep(1/7, 7)) # pro beta
t(a)%*%coef(model.year)

2*pt(-abs(t(a)%*%coef(model.year)/sqrt(t(a)%*%vcov(model.year)%*%a)), df=24)
t(a)%*%coef(model.year)+c(-1, 1)*sqrt(t(a)%*%vcov(model.year)%*%a)

# (c)
# ---
# podzim vs rok
a <- c(-1/12, rep(1/3-1/12, 3), rep(-1/12, 8)) # pro mu
a.p <- c(0, rep(1/3-1/12, 3), rep(-1/12, 8)) # pro beta
t(a.p)%*%coef(model.year)

# zima vs rok
a <- c(rep(-1/12, 4), rep(1/3-1/12, 3), rep(-1/12, 5)) # pro mu
a.z <- c(0, rep(-1/12, 3), rep(1/3-1/12, 3), rep(-1/12, 5)) # pro beta
t(a.z)%*%coef(model.year)

# jaro vs rok
a <- c(rep(-1/12, 7), rep(1/3-1/12, 3), rep(-1/12, 2)) # pro mu
a.j <- c(0, rep(-1/12, 6), rep(1/3-1/12, 3), rep(-1/12, 2)) # pro beta
t(a.j)%*%coef(model.year)

# leto vs rok
a <- c(1/3-1/12, rep(-1/12, 9), rep(1/3-1/12, 2)) # pro mu
a.l <- c(0, rep(-1/12, 9), rep(1/3-1/12, 2)) # pro beta
t(a.l)%*%coef(model.year)

A <- rbind(a.j, a.l, a.p, a.z)

tab.seasons <- data.frame(row.names=c("Spring", "Summer", "Autumn", "Winter"), Estimate=rep(NA, 4), StdError=NA, Tvalue=NA, Pvalue=NA)
tab.seasons[, 1] <- A%*%coef(model.year)
tab.seasons[, 2] <- sqrt(diag(A%*%vcov(model.year)%*%t(A)))
tab.seasons[, 3] <- tab.seasons[, 1]/tab.seasons[, 2]
tab.seasons[, 4] <- 2*pt(-abs(tab.seasons[, 3]), df=summary(model.year)$df[2])
tab.seasons

ci.seasons <- data.frame(row.names=c("Spring", "Summer", "Autumn", "Winter"), Lower=NA, Estimate=tab.seasons$Estimate, Upper=NA)
ci.seasons$Lower <- tab.seasons$Estimate-qt(0.975, df=24)*tab.seasons$StdError
ci.seasons$Upper <- tab.seasons$Estimate+qt(0.975, df=24)*tab.seasons$StdError
ci.seasons

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

# (a)
# ---
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)

# (b)
# ---
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))

# (c)
# ---
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

# (d)
# ---
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")
}