[R] Doubts on anova and use of contrasts in multcomp package
José Cláudio Faria
joseclaudio.faria at terra.com.br
Mon May 31 21:42:03 CEST 2004
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
More information about the R-help
mailing list