# R code to analyze data from Latin Square Design # We use the bank teller productivity data from class and chapter 28 ############################ # # Some code related to designing a Latin Square: # # A random permutation of the numbers 1, 2, 3, 4, 5: sample(1:5,replace=F) # Another random permutation of the numbers 1, 2, 3, 4, 5: sample(1:5,replace=F) # ############################# # Save the data file into a directory and # use the full path name: bankteller.data <- read.table(file = "z:/My Documents/teaching/stat_705/banktellerdata.txt", header=FALSE, col.names = c('productivity', 'week', 'day', 'music')) # attaching the data frame: attach(bankteller.data) # Making "week", "day" and "music" factors: week <- factor(week) day <- factor(day) music <- factor(music) # Fitting the model and viewing the ANOVA table: bank.fit <- lm(productivity ~ week + day + music) anova(bank.fit) # Sample mean productivity for each level of music: music.means <- tapply(productivity, music, mean); print(music.means) # Main effects plot to visually examine differences # in mean response across treatments: plot(as.numeric(levels(music)), music.means ,type='p', main="Main Effects Plot, Kenton Data Set") lines(as.numeric(levels(music)), music.means ,lty=2) abline(h=mean(productivity),lty=3) # Note we are calculating and plotting the mean productivity # at each level of "music" here. The 20.6 value # is the OVERALL sample mean productivity. # ############################################################ # Tukey CIs for pairwise treatment mean differences: TukeyHSD(aov(bank.fit),conf.level=0.90)$music # NOTE: The CIs which do NOT contain zero indicate the treatment means # that are significantly different at (here) the 0.10 family significance level. # We see from the ANOVA output that there are significant # differences among the music types in terms of mean productivity. # (F* = 10.58, P-value = 0.0007) The Tukey output shows which # pairs of music types are significantly different from each other # in terms of mean productivity. # ############### Checking Model Assumptions: ############### # Residual Plots and Q-Q plots: # Residual plots for checking ANOVA model assumptions: # plotting residuals versus fitted values: plot(fitted.values(bank.fit), residuals(bank.fit), xlab="Fitted Values", ylab="Residuals"); abline(h=0) # A normal Q-Q plot of the residuals: qqnorm(residuals(bank.fit))