# Ruback and Juieng # This page starts with a comparison of orthogonal and nonorthogonal, crossed # with L or Linv. # change last 2 items from 44.05 and 66.46 to 94.05 and 120.46 for tails time <- c(33.98, 37.50, 20.24, 29.80, 58.22, 33.60, 57.80, 42.15, 48.63, 47.92, 26.03, 55.03, 31.77, 3.85, 14.24, 57.65, 47.66, 48.16, 21.90, 20.88, 0.10, 36.27, 24.38, 29.59, 39.32, 48.72, 26.36, 10.51, 29.48, 33.86, 5.91, -5.51, 28.22, 52.89, 36.84, 29.05, 28.50, 82.78, 28.51, 53.77, 16.63, 46.77, 25.51, 37.66, 36.08, 51.89, 60.22, 17.56, 7.42, 21.68, 49.70, 52.99, 57.51, 49.77, 45.27, 33.64, 42.43, 47.53, 33.66, 52.20,70.48, 48.33, 45.69, 52.39, 48.40, 27.27, 64.67, 62.58, 37.51, 33.87, 59.99, 57.67, 43.61, 44.05, 66.46 ) group <- c(rep ("NoWait", 25), rep("Wait",25), rep("Honk",25)) group <- factor(group) nreps <- 1000 F.resamp <- numeric(nreps) F.resampC1 <- numeric(nreps) F.resampC2 <- numeric(nreps) model1 <- aov(time~group) model2 <- summary(model1, split = list(group = list("First and third vs middle"=1, "first vs last" = 2))) Fobt <- model2[[1]][[4]][1] # The overall F for groups par(mfrow = c(2,1)) for (i in 1:nreps) { resamp.time <- sample(time, length(time), replace = FALSE) model.resamp <- aov(resamp.time~group) F.resamp[i] <- summary(model.resamp, split = list(group = list("First and third vs middle"=1, "first vs last" = 2))) [[1]][[4]][1] #Distrib of overall F for groups } probF <- length(F.resamp[F.resamp >= Fobt])/nreps Fobt <- round(Fobt, digits = 2) probF <- round(probF, digits = 3) hist(F.resamp, mai.1) ; ht <- df(dd,2,72) dd <- c(0,dd) ; ht <- c(0,ht) polygon(dd,ht) n = "Distrib. of Overall F Under Null", freq = FALSE, xlab = "Resampled F", xlim = c(0,12), breaks = 50 ) dd <- seq(0,12, text(4, .4, paste("F = ",Fobt)) text(6, .4, paste("p = ",probF)) cat("The probability for overall F under the null is = ", probF, '\n') # The following pulls out the 1st contrast FobtC1 <- model2[[1]][[4]][2] # The First F for Contrasts for (i in 1:nreps) { resamp.time <- sample(time, length(time), replace = FALSE) model.resamp <- aov(resamp.time~group) F.resampC1[i] <- summary(model.resamp, split = list(group = list("First and third vs middle"=1, "first vs last" = 2))) [[1]][[4]][2] } probFC1 <- length(F.resamp[F.resamp >= FobtC1])/nreps FobtC1 <- round(FobtC1, digits = 2) probF <- round(probF, digits = 3) hist(F.resampC1, main = "Distrib. of First F Contrast Under Null", freq = FALSE, xlab = "Resampled F Contrasts#1", xlim = c(0,12), breaks = 50 ) text(4, .4, paste("F = ",FobtC1)) text(6, .4, paste("p = ",probFC1)) cat("The probability for contrast 1 under the null is = ", probFC1, '\n') ### At this point I could look at other contrasts ######################################### ### The following stuff is useful, but extra. contrasts(group) <- matrix(c(-2,1,1, 0,-1,1), ncol=2) mod <- lm(time ~ group, contrasts = contrasts(group)) # This is orthog using inverse L <- matrix(c(1/3, 1/3, 1/3, 1/2, -1, 1/2, 1, 0, -1), ncol = 3) LInv <- solve(t(L)) orth.Contrasts <- LInv[,2:3] contrasts(group) = orth.Contrasts model5 <- lm(time~group) summary(model5,split=list(Group=list("1 & 3 vs 2"=1, "1 vs 3" = 2))) # This is orthog without inverse L <- matrix(c(1/3, 1/3, 1/3, 1/2, -1, 1/2, 1, 0, -1), ncol = 3) orth.Contrasts <- L[,2:3] contrasts(group) = orth.Contrasts model6 <- lm(time~group) summary(model6,split=list(group=list("1&3 vs 2"=1, "1 vs 3" = 2))) #Now make them nonorthogonal # With inverse L <- matrix(c(1/3, 1/3, 1/3, 1/2, -1, 1/2, 1, -1, 0), ncol = 3) LInv <- solve(t(L)) orth.Contrasts <- LInv[,2:3] contrasts(group) = orth.Contrasts model7 <- lm(time~group) summary(model7,split=list(Group=list("1 & 3 vs 2"=1, "1 vs 3" = 2))) #Nonorthogonal without inverse L <- matrix(c(1/3, 1/3, 1/3, 1/2, -1, 1/2, 1, -1, 0), ncol = 3) Nonorth.Contrasts <- L[,2:3] contrasts(group) = Nonorth.Contrasts model8 <- lm(time~group) summary(model8,split=list(group=list("1&3 vs 2"=1, "1 vs 3" = 2)))