José Cláudio Faria
2004-May-31 19:42 UTC
[R] Doubts on anova and use of contrasts in multcomp package
Dear list, I have been studying R and I would like the aid of more experienced to solve the problems of the analysis below: r = gl(3, 8, label = c('r1', 'r2', 'r3')) e = rep(gl(2, 4, label = c('e1', 'e2')), 3) y = c(26.2, 26.0, 25.0, 25.4, 24.8, 24.6, 26.7, 25.2, 25.7, 26.3, 25.1, 26.4, 19.6, 21.1, 19.0, 18.6, 22.8, 19.4, 18.8, 19.2, 19.8, 21.4, 22.8, 21.3) df = data.frame(r, e, y) aux = sort(rep(letters[1:6], 4)) #auxiliary variable df = data.frame(df, aux) attach(df) par(mfrow = c(2, 1)) interaction.plot(r, e, y, col = 'blue', ylab = 'y', xlab = 'r') interaction.plot(e, r, y, col = 'blue', ylab = 'y', xlab = 'r') av1 = aov(y ~ r*e) av2 = aov(y ~ r/e) efR_E = summary(av2, split = list('r:e' = list( 'e1 vs e2/r1' = 1, 'e1 vs e2/r2' = 2, 'e1 vs e2/r3' = 3))) av3 = aov(y ~ e/r) efE_R = summary(av3, split = list('e:r' = list( 'r/e1' = c(1,3), 'r/e2' = c(2,4)))) mds = model.tables(av1, ty = 'means') detach(df) cat('\nData:'); cat('\n') print(df) cat('\nMeans:'); cat('\n') print(mds) cat('\nANOVA:'); cat('\n') print(summary(av1)); cat('\n') cat('\nANOVA - E effect in R levels:'); cat('\n') print(efR_E); cat('\n') cat('\nANOVA - R effect in E levels:'); cat('\n') print(efE_R); cat('\n') #===Below: my original intention (post in this list and still not answered...) # ANOVA - R effect in E levels: #---------------------------------------------- # Df # e 1 # e:r 4 # e:r: r/e1 2 # r1 vs (r2,r3)/e1 1 ?... # r2 vs r3/e1 1 ?... # e:r: r/e2 2 # r1 vs (r2,r3)/e1 1 ?... # r2 vs r3/e2 1 ?... #---------------------------------------------- #Residuals 18 #---------------------------------------------- #===Below: alternative using multcomp # (with auxiliary variable - aux) to study R effect in E levels: # a: r1/e1 # c: r2/e1 # e: r3/e1 # b: r1/e2 # d: r2/e2 # f: r3/e2 #a b c d e f C1 = c(2, 0, -1, 0, -1, 0) # r1 vs (r2,r3)/e1 C2 = c(0, 0, 1, 0, -1, 0) # r2 vs r3/e1 C3 = c(0, 2, 0, -1, 0, -1) # r1 vs (r2,r3)/e2 C4 = c(0, 0, 0, 1, 0, -1) # r2 vs r3/e2 C = rbind(C1, C2, C3, C4) row.names(C) = c('r1 vs (r2,r3)/e1', 'r2 vs r3/e1', 'r1 vs (r2,r3)/e2', 'r2 vs r3/e2') lim1 = lm(y ~ aux, data = df) print(anova(lim1)) tc1 = simtest(y ~ aux, data = df, conf.level = 0.9, alternative = 'less', eps = 1e-04, cmatrix = C) print(summary(tc1)) #===Below: verifying E effect in R levels (already analized in av2) # (with auxiliary variable - aux) # a: e1/r1 # c: e1/r2 # e: e1/r3 # b: e2/r1 # d: e2/r2 # f: e2/r3 #a b c d e f C1 = c(1, -1, 0, 0, 0, 0) # e1 vs e2/r1 C2 = c(0, 0, 1, -1, 0, 0) # e1 vs e2/r2 C3 = c(0, 0, 0, 0, 1, -1) # e1 vs e2/r3 C = rbind(C1, C2, C3) row.names(C) = c('e1 vs e2/r1', 'e1 vs e2/r2', 'e1 vs e2/r3') lim2 = lm(y ~ aux, data = df) print(anova(lim2)) tc2 = simtest(y ~ aux, data = df, conf.level = 0.9, alternative = 'less', eps = 1e-04, cmatrix = C) print(summary(tc2)) #===My Questions: # a) Is possible the resolution of the original intention? How? # b) Why p-values of soluctions av2 and lim2 dont agree? # c) Are there another better way to lead of this analysis? #================================================================== Best regards, Jos?? Cl??udio Faria UESC/DCET Brasil 73-634.2779 joseclaudio.faria at terra.com.br jc_faria at uol.com.br