[R] Doubts on anova and use of contrasts in multcomp package

[José Cláudio Faria]

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