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