[banner]

An R Companion for the Handbook of Biological Statistics

Salvatore S. Mangiafico

Two-way Anova with Robust Estimation

A two-way anova using robust estimators can be performed with the WRS2 package.  Options for estimators are M-estimators, trimmed means, and medians.  This type of analysis is resistant to deviations from the assumptions of the traditional ordinary-least-squares anova, and are robust to outliers.  However, it may not be appropriate for data that deviate too widely from parametric assumptions.

 

The main analysis using M-estimators for a two-way anova is conducted with the pbad2way function in the WRS2 package.   Post-hoc tests can be performed with the mcp2a function in the WRS2 package or with my custom functions pairwise.robust.test and pairwise.robust.matrix, which rely on the pb2gen function in WRS2.

 

My custom function groupwise.huber uses the HuberM function in the DescTools package  to determine the Huber M-estimators across groups in a data frame.

 

For more information on robust tests available in the WRS2 package, see:

 

help(package="WRS2")

 

Consult the chapter on Two-way Anova for general consideration about conducting analysis of variance.

 

### --------------------------------------------------------------
### Two-way anova with robust estimators, hypothetical data
###   Using WRS2 package
### --------------------------------------------------------------

Input = (

"Factor.A  Factor.B  Response
 l         x          0.9
 l         y          1.4
 l         x          1.3
 l         y          2.0
 l         x          1.6
 l         y          2.6
 m         x          2.4
 m         y          3.6
 m         x          2.8
 m         y          3.7
 m         x          3.2
 m         y          3.0
 n         x          1.6
 n         y          1.2
 n         x          2.0
 n         y          1.9
 n         x          2.7
 n         y          0.9
")

Data = read.table(textConnection(Input),header=TRUE)

 

Produce Huber M-estimators and confidence intervals by group

 

library(rcompanion)

Sum = groupwiseHuber(data = Data,
                     group = c("Factor.A", "Factor.B"),
                     var = "Response",
                     conf.level=0.95,
                     conf.type="wald")

Sum

  Factor.A Factor.B n  M.Huber  lower.ci upper.ci

1        l        x 3 1.266667 0.9421910 1.591142

2        l        y 3 2.000000 1.4456385 2.554362

3        m        x 3 2.800000 2.4304256 3.169574

4        m        y 3 3.538805 3.2630383 3.814572

5        n        x 3 2.100000 1.5855743 2.614426

6        n        y 3 1.333333 0.8592063 1.807460

 

 

 

Interaction plot using summary statistics

 

library(ggplot2)

pd = position_dodge(.2)

ggplot(Sum, aes(x=Factor.A,
                y=M.Huber,
                color=Factor.B)) +
    geom_errorbar(aes(ymin=lower.ci,
                      ymax=upper.ci),
                   width=.2, size=0.7, position=pd) +
    geom_point(shape=15, size=4, position=pd) +
    theme_bw() +
    theme(
          axis.title.y = element_text(vjust= 1.8),
          axis.title.x = element_text(vjust= -0.5),
          axis.title = element_text(face = "bold")) +
    scale_color_manual(values = c("black", "blue"))

 

 

 

 

Two-way analysis of variance for M-estimators

The est = "mom" option uses a modified M-estimator for the analysis.  To analyze using medians, use the est= "median" option in the pbad2way function in the WRS2 package.  To analyze using trimmed means, use the t2way function in the WRS2 package.


library(WRS2)

pbad2way(Response ~ Factor.A + Factor.B + Factor.A:Factor.B,
         data = Data,
         est = "mom",    # modified M-estimator
         nboot = 5000)   # number of bootstrap samples
                         # a higher number will take longer to compute

 

pbad2way(formula = Response ~ Factor.A + Factor.B + Factor.A:Factor.B,

    data = Data, est = "mom", nboot = 3000)

 

                  p.value

Factor.A           0.0000

Factor.B           0.3403

Factor.A:Factor.B  0.0460

 

 

Produce post-hoc tests for main effects with mcp2a

 

post = mcp2a(Response ~ Factor.A + Factor.B + Factor.A:Factor.B,
             data = Data,
             est = "mom",    # M-estimator
             nboot = 5000)   # number of bootstrap samples

post$contrasts

post

 

                                             Factor.A1:  Factor.A2:  Factor.A3:                       

    Factor.A1 Factor.A2 Factor.A3 Factor.B1  Factor.B1   Factor.B1   Factor.B1

l_x         1         1         0         1          1           1           0

l_y         1         1         0        -1         -1          -1           0

m_x        -1         0         1         1         -1           0           1

m_y        -1         0         1        -1          1           0          -1

n_x         0        -1        -1         1          0          -1          -1

n_y         0        -1        -1        -1          0           1           1

 

 

                          V1 ci.lower ci.upper p-value

Factor.A1           -3.18333 -4.20000 -1.60000 0.00000

Factor.A2           -0.16667 -1.70000  1.36667 0.40233

Factor.A3            3.01667  1.40000  4.05000 0.00000

Factor.B1           -0.81667 -2.28333  1.00000 0.22233

Factor.A1:Factor.B1  0.11667 -1.50000  1.16667 0.48033

Factor.A2:Factor.B1 -1.50000 -3.10000  0.00000 0.01767

Factor.A3:Factor.B1 -1.61667 -2.80000  0.00000 0.01433

 

### The Factor.A1 contrast compares l to m; since it is significant,

###   l is significantly different than m.

### The Factor.A2 contrast compares l to n; since it is not significant,

###   l is not significantly different than n.

 

 

Produce post-hoc tests for main effects with pairwiseRobustTest or pairwiseRobustMatrix

 

Table output with pairwiseRobustTest

 

library(rcompanion)

PT = pairwiseRobustTest(
     Data$Response,
     Data$Factor.A,
     est="mom",
     nboot=5000,
     method="fdr")      # adjust p-values; see ?p.adjust for options

PT

 

  Comparison Statistic p.value p.adjust

1  l - m = 0    -1.483   8e-04   0.0012

2  l - n = 0  -0.08333  0.7226   0.7226

3  m - n = 0       1.4   6e-04   0.0012

### p-values may differ

 

 

Compact letter display output with pairwiseRobustMatrix

 

library(rcompanion)

PM = pairwiseRobustMatrix(
     Data$Response,
     Data$Factor.A,
     est="mom",
     nboot=5000,
     method="fdr")      # adjust p-values; see ?p.adjust for options

PM$Adjusted

 

PM$Adjusted

       l     m      n

l 1.0000 9e-04 0.7284

m 0.0009 1e+00 0.0009

n 0.7284 9e-04 1.0000

### p-values may differ

 


library(multcompView)

multcompLetters(PM$Adjusted,
                compare="<",
                threshold=0.05,
                Letters=letters,
                reversed = FALSE)

 

  l   m   n

"a" "b" "a"

### Note, means are not ordered from largest to smallest

 

 

Produce post-hoc tests for interaction effect

 

Data$Factor.int = interaction (Data$Factor.A, Data$Factor.B)

   ### Create a factor which is the interaction of Factor.A and Factor.B

library(FSA)

headtail(Data)

 

   Factor.A Factor.B Response Factor.int

1         l        x      0.9        l.x

2         l        y      1.4        l.y

3         l        x      1.3        l.x

16        n        y      1.9        n.y

17        n        x      2.7        n.x

18        n        y      0.9        n.y

 


Table output with pairwiseRobustTest


library(rcompanion)

PT = pairwise.robust.test(
     Data$Response,
     Data$Factor.int,
     est="mom",
     nboot=5000,
     method="fdr")      # adjust p-values; see ?p.adjust for options

PT

 

      Comparison Statistic p.value p.adjust

1  l.x - l.y = 0   -0.7333  0.1342   0.1629

2  l.x - m.x = 0    -1.533       0   0.0000

3  l.x - m.y = 0    -2.383       0   0.0000

4  l.x - n.x = 0   -0.8333  0.0687   0.1472

5  l.x - n.y = 0  -0.06667  0.9457   0.9457

6  l.y - m.x = 0       0.8  0.1298   0.1629

7  l.y - m.y = 0     -1.65       0   0.0000

8  l.y - n.x = 0       0.1  0.7681   0.8230

9  l.y - n.y = 0    0.6667  0.1368   0.1629

10 m.x - m.y = 0     -0.85  0.1314   0.1629

11 m.x - n.x = 0       0.7  0.1374   0.1629

12 m.x - n.y = 0     1.467       0   0.0000

13 m.y - n.x = 0     -1.55       0   0.0000

14 m.y - n.y = 0     2.317       0   0.0000

15 n.x - n.y = 0    0.7667  0.1412   0.1629

### p-values may differ

 

 

Compact letter display output with pairwiseRobustMatrix


library(rcompanion)

PM = pairwiseRobustMatrix(
     Data$Response,
     Data$Factor.int,
     est="mom",
     nboot=5000,
     method="fdr")      # adjust p-values; see ?p.adjust for options

PM

 

$Unadjusted

    l.x    l.y    m.x    m.y    n.x    n.y

l.x  NA 0.1322 0.0000 0.0000 0.0666 0.9515

l.y  NA     NA 0.1366 0.0000 0.7465 0.1324

m.x  NA     NA     NA 0.1418 0.1242 0.0000

m.y  NA     NA     NA     NA 0.0000 0.0000

n.x  NA     NA     NA     NA     NA 0.1396

n.y  NA     NA     NA     NA     NA     NA

 

$Method

[1] "fdr"

 

$Adjusted

       l.x    l.y    m.x    m.y    n.x    n.y

l.x 1.0000 0.1636 0.0000 0.0000 0.1427 0.9515

l.y 0.1636 1.0000 0.1636 0.0000 0.7998 0.1636

m.x 0.0000 0.1636 1.0000 0.1636 0.1636 0.0000

m.y 0.0000 0.0000 0.1636 1.0000 0.0000 0.0000

n.x 0.1427 0.7998 0.1636 0.0000 1.0000 0.1636

n.y 0.9515 0.1636 0.0000 0.0000 0.1636 1.0000

### p-values may differ

 

 

library(multcompView)

multcompLetters(PM$Adjusted,
                compare="<",
                threshold=0.05,
                Letters=letters,
                reversed = FALSE)

 

l.x  l.y  m.x  m.y  n.x  n.y

 "a" "ab" "bc"  "c" "ab"  "a"

### Note, means are not ordered from largest to smallest

 

#     #     #