########################### ### Bang & Olufsen data ### ########################### # 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) ############################### ### Input and plot toy data ### ############################### # 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) 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))) # Box plots by treatment plot(treatm, y, xlab = "Treatment", ylab = "y") # Box plots by block plot(block, y, xlab = "Block", ylab="y") #################################################### ### Compute estimates of parameters in the model ### #################################################### # Sample mean (mu_hat <- mean(y)) # Sample mean deviation for each treatment (alpha_hat <- tapply(y, treatm, mean) - mu_hat) # Sample mean deviation for each block (beta_hat <- tapply(y, block, mean) - mu_hat) ######################################################## ### Compute the total variation, sum of squares SST ### ######################################################## # SST for the example data (SST <- sum((y - mu_hat)^2)) ##################################################################### ### Compute variance explained by the treatment part of the model ### ##################################################################### # Sum of squares of treatment SS(Tr) for the example data (SSTr <- l * sum(alpha_hat^2)) ################################################################# ### Compute variance explained by the block part of the model ### ################################################################# # Sum of squares for blocks SS(Bl) for the example data (SSBl <- k * sum(beta_hat^2)) ##################################################### ### Compute residual variance after the model fit ### ##################################################### # Residual sum of squares SSE (SSE <- SST - SSTr - SSBl) ########################################## ### Hypothesis test: Treatment effect? ### ########################################## # Plot density of relevant F-distribution. Remember that this is "under H0" # (computed as if H0 were true) xseq <- seq(0, 10, by = 0.1) plot(xseq, df(xseq, df1 = k-1, df2 = (k-1)*(l-1)), type = "l") # Show critical value (5% signif. level) for test of treatment hypothesis critical_value <- qf(0.95, df1 = k-1, df2 = (k-1)*(l-1)) abline(v = critical_value, col = "red") # Compute value of the test statistic (FTr <- (SSTr/(k-1)) / (SSE/((k-1)*(l-1)))) # Compute p-value for the test 1 - pf(FTr, df1 = k-1, df2 = (k-1)*(l-1)) ###################################### ### Hypothesis test: Block effect? ### ###################################### # Plot density of relevant F-distribution. Remember that this is "under H0" # (computed as if H0 were true) xseq <- seq(0, 10, by = 0.1) plot(xseq, df(xseq, df1 = l-1, df2 = (k-1)*(l-1)), type = "l") # Show critical value (5% signif. level) for test of treatment hypothesis critical_value <- qf(0.95, df1 = l-1, df2 = (k-1)*(l-1)) abline(v = critical_value, col = "red") # Compute value of the test statistic (FBl <- (SSBl/(l-1)) / (SSE/((k-1)*(l-1)))) # Compute p-value for the test 1 - pf(FBl, df1 = l-1, df2 = (k-1)*(l-1)) ################################################### ### Two-way ANOVA output using anova() and lm() ### ################################################### anova(lm(y ~ treatm + block)) ############################ ### Model validation: ### ### Variance homogeneity ### ############################ # 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") ######################### ### Model validation: ### ### Normality ### ######################### # Normal QQ-plot of the residuals qqnorm(fit\$residuals) qqline(fit\$residuals)