install.packages("lmerTest")

# Get the B&O data from the lmerTest-package
library(lmerTest)
data(TVbo)
# Each of 8 assessors scored each of 12 combinations 2 times.
# Take a look at the sharpness scores for one single picture
# and one of the two repetitions
TVbo_sub <- subset(TVbo, Picture == 1 & Repeat == 1)[, c(1, 2, 9)]
sharp <- matrix(TVbo_sub$Sharpness, nrow = 8, byrow = T)
colnames(sharp) <- c("TV3", "TV2", "TV1")
rownames(sharp) <- c("Person 1", "Person 2", "Person 3",
                     "Person 4", "Person 5", "Person 6",
                     "Person 7", "Person 8")
library(xtable)
xtable(sharp) ## Outputter i latex

sharp

## Envejs anova
anova(lm(Sharpness ~ TVset, data = TVbo_sub))


## Tovejs anova
anova(lm(Sharpness ~ TVset + Assessor, data = TVbo_sub))

## samme konklusion i dette tilfælde, med lidt forskellige statistiske resultater

######

# Observations
y <- c(2.8, 3.6, 3.4, 2.3,
       5.5, 6.3, 6.1, 5.7,
       5.8, 8.3, 6.9, 6.1)
# Treatments (groups, varieties)
## HUSK: treatm skal være factor
treatm <- factor(c(1, 1, 1, 1,
                   2, 2, 2, 2,
                   3, 3, 3, 3))
# Blocks (persons, fields)
block <- factor(c(1, 2, 3, 4,
                  1, 2, 3, 4,
                  1, 2, 3, 4))
# No. of treatments and no. of blocks (for later formulas)
(k <- length(unique(treatm)))
(l <- length(unique(block)))
par(mfrow = c(1,2))
# Box plots by treatment
plot(treatm, y, xlab = "Treatment", ylab = "y")
# Box plots by block
plot(block, y, xlab = "Block", ylab="y")

#### Estimater
# Sample mean
(mu_hat <- mean(y))

tapply(y, treatm, mean) ### giver gruppegennemsnit
(alpha_hat <- tapply(y, treatm, mean) - mu_hat)

# Sample mean deviation for each block
tapply(y, block, mean) ### giver gruppegennemsnit
(beta_hat <- tapply(y, block, mean) - mu_hat)

k <- 3
l <- 4

SST <- sum((y - mu_hat)^2)
SST

SS_tr <- sum(alpha_hat^2) * l
SS_tr

SS_bl <- sum(beta_hat^2) * k
SS_bl

SSE <- SST - SS_tr - SS_bl
SSE

# Compute value of the test statistic
(FTr <- (SS_tr/(k-1)) / (SSE/((k-1)*(l-1))))
# Compute p-value for the test
1 - pf(FTr, df1 = k-1, df2 = (k-1)*(l-1))

## Tilsvarende for block:
# Compute value of the test statistic
(FBl <- (SS_bl/(l-1)) / (SSE/((k-1)*(l-1))))
# Compute p-value for the test
1 - pf(FBl, df1 = l-1, df2 = (k-1)*(l-1))

## tovejs ANOVA i R 
anova(lm(y ~  treatm + block))


#### Modelkontrol
# Save the fitted model
fit <- lm(y ~ treatm + block)
# Make box plots of residuals
par(mfrow = c(1,2))
plot(treatm, fit$residuals, xlab = "Treatment")
plot(block, fit$residuals, xlab = "Block")

qqnorm(fit$residuals)
qqline(fit$residuals)

## residualplot


plot(fitted(fit), residuals(fit), xlab = "Treatment")

MESS::wallyplot(fit) ## Det ser ok ud. 











D <- data.frame(
  y=c(22.5, 24.3, 24.9, 22.4,
      21.5, 21.3, 23.9, 18.4,
      22.2, 21.9, 21.7, 17.9),
  car=c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4),  ## FORKERT, vi skal ikke lave (multipel) lineær regression
  tire=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3)
)


D <- data.frame(
  y=c(22.5, 24.3, 24.9, 22.4,
      21.5, 21.3, 23.9, 18.4,
      22.2, 21.9, 21.7, 17.9),
  car=factor(c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4)), ## VIGTIGT at det er factor
  tire=factor(c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3))
)

D ## Tjek at faktorer og y-værdier stemmer overens


par(mfrow=c(1,2))
plot(D$tire, D$y, xlab="Tire", ylab="y")
plot(D$car, D$y, xlab="Car", ylab="y")


fit <- lm(y ~ car + tire, data=D)
anova(fit)
## Konklusion: begge variable signifikante

## Modelkontrol
qqnorm(fit$residuals)
qqline(fit$residuals)

par(mfrow=c(1,2))
plot(D$car, fit$residuals, xlab="Car", ylab="Residuals")
plot(D$tire, fit$residuals, xlab="Tire", ylab="Residuals")



### Post-Hoc
## Ingen a priori valgte sml, af dæk, så vi tager alle tre. 


tapply(D$y, D$tire, mean) ## Kører længst med dæk nr. 1

M <- 3
alpha <- 0.05
alpha / M
1 - alpha / M / 2

k <- 3
l <- 4
(k-1)*(l-1)


LSD_bonf <- qt(1-0.05/6, df=6) * sqrt(2*1.19/4)
LSD_bonf