Interpreting the output of sppba(), sppbb(), and sppbi() in R

Question

I'm trying to run a robust mixed ANOVA with a between groups factor (2 levels) and a within groups factor (3 levels). I'd ideally like to be able to set and run planned contrasts, have a full set of relevant summary statistics (presumably I'd report summary stats for the trimmed means rather than the original means..), and then calculate effect sizes for the contrasts.

I'm new to R (I attempted to teach it to myself over the weekend) so if you can explain R language stuff at the most basic level imaginable then that would be appreciated..

Following Field (Discovering Statistics using R, 2012) and Wilcox (Introduction to Robust Estimation and Hypothesis Testing, 2013) I've downloaded Rallfun-v40.txt to my working directory from here and then brought its contents into R using source(file.choose()).

From this I've followed Field's instructions and run the following on my data (converted to wide format such that the columns are betweenL1_withinL1 betweenL1_withinL2 ... betweenL2_withinL3 )

bwtrim(2, 3, data, tr = .2) # using 20% trimmed means

And then I've also run the following three lines of code:

sppba(2, 3, data, est = mom, nboot = 2000) 
sppbb(2, 3, data, est = mom, nboot = 2000) 
sppbi(2, 3, data, est = mom, nboot = 2000) 

The output from bwtrim() is straightforward enough to understand. It gives me test statistics (Q) for both factors and their interaction, each accompanied by a p-value.

I don't really understand what the output from the remaining three lines gives me. The R-documentation tells me "The functions sppba, sppbb, and sppbi compute the main fixed effect, the main within-subjects effect, and the interaction effect only, respectively, using bootstrap."

It seems to return test statistics (psi hat), a p-value, and a table of contrast codes as follows:

sppba(2, 3, data, est = mom, nboot = 2000)

$p.value
[1] 0.112

$psihat
[1] -5.577981

$con
     [,1]
[1,]    1
[2,]   -1

sppbb(2, 3, data, est = mom, nboot = 2000)

$p.value
[1] 0.4562719

$center
[1]  2.9931670 -0.4448510  0.5338472

sppbi(2, 3, data, est = mom, nboot = 2000)

$p.value
[1] 0.58

$psihat
[1] -5.345118 -5.875428  1.012981

$con
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
[4,]   -1    0    0
[5,]    0   -1    0
[6,]    0    0   -1 

Can you explain this output to me? To whatever extent the above code doesn't already do this, can you also give me advice on my initial aim: "to set and run planned contrasts, have a full set of relevant summary statistics (presumably I'd report summary stats for the trimmed means rather than the original means..), and then calculate effect sizes for the contrasts"

Finally, I also ran a robust mixed anova with bootstrapping and a 20% trimmed mean (hopefully I inputted the arguments correctly...) as follows:

bwtrimbt(2, 3, data, tr=0.2,JK=2*3,nboot=2000)

This had the result of making a significant interaction from the original bwtrim() calculation only marginally significant. How would you interpret this? Am I better off using bwtrimbt? The data for each of the columns in my wide format dataset is not normal and the (between) group sizes are n=17 and n=21.

Answer 1

After a few days of reading up on this, I'm going to attempt to answer my own question. However, any corrections to my answer will be appreciated.

bwtrim(), sppba(), sppbb(), and sppbi() seem to give different outputs depending on:

1. whether they are called from the WRS2 package or Rallfun-v40.txt (from https://dornsife.usc.edu/labs/rwilcox/software)
2. whether the arguments are in one of the following two formats:

• bwtrim(outcome ~ between*within, id = id, data = data_long, tr=0.2)

• bwtrim(2, 3, data = data_wide, tr = .2)

From what I can tell, the WRS2 package requires the arguments to be as follows:

bwtrim(outcome ~ between*within, id = id, data = data_long, tr=0.2)

The data should be in long format, i.e., 4 columns:

• Column 1: "id" [observation number (class: integer)]

• Column 2: "between" [between subjects variable (class: factor)]

• Column 3: "within" [within subjects variable (class: factor)]

• Column 4: "outcome" [outcome variable (class: numeric)]

Whereas for the Rallfun-v40.txt version the arguments should (for a 2x3 mixed design (between x within)) be in the following format:

bwtrim(2, 3, data = data_wide, tr = .2)

For this, the data should be in wide format, i.e. (with 'b' for between subjects factor and 'w' for within subjects factor):

• Column 1: b1w1

• Column 2: b1w2

• Column 3: b1w3

• Column 4: b2w1

• Column 5: b2w2

• Column 6: b2w3

The WRS2 version provides more easily interpretable output:

#WRS2 Output for bwtrim(outcome ~ between*within, id = id, data = data_long, tr=0.2)

#library("WRS2")

#                value   df1  df2      p.value
#between         2.6551   1   19.9949  0.1189
#within          1.4066   2   16.3158  0.2732
#between:within  3.2474   2   16.3158  0.0650

vs

#source(file.choose()) #then open Rallfun-v40.txt from working directory

#$Qa #i.e., the test statistic (Q) for the between subjects factor
#[1] 2.655095

#$Qa.p.value
#          [,1]
#[1,] 0.1162786

#$Qb #i.e., the test statistic (Q) for the within subjects factor
#[1] 1.415902

#$Qb.p.value
#          [,1]
#[1,] 0.2664895

#$Qab #i.e., the test statistic (Q) for the interaction
#[1] 3.26877

#$Qab.p.value
#           [,1]
#[1,] 0.05957826

For sppba(), sppbb(), and sppbi() the output is also easier to interpret when called from the WRS2 package:

WRS2 output from sppbb(outcome ~ between*within, id, data = data_long)

#Test statistics:
#                                            Estimate
#within2-within1                            -2.9932
#within2-within3                             0.5338
#within1-within3                            -0.4449

#Test whether the corresponding population parameters are the same:
#p-value: 0.454 

It tests whether the “typical” difference score (as measured by an M-estimator) between any two levels of measurement occasions is 0 (while ignoring the between-subjects groups), resulting in a p-value.

--

WRS2 output from sppba(outcome ~ between*within, id = id, data = data_long, avg = FALSE)

#Test statistics:
#                                      Estimate
#within2      between1-between2          0.368
#within1      between1-between2         -11.250
#within3      between1-between2         -5.852

#Test whether the corresponding population parameters are the same:
#p-value: 0.114 

This performs pairwise group comparisons for the between subjects groups in each of the within subjects conditions, and then tests whether all of the null hypotheses are simultaneously true, resulting in a p-value.

sppba(outcome ~ between*within, id = id, data = data_long)

#Test statistics:
#             Estimate
#names-notes   -5.578

#Test whether the corresponding population parameters are the same:
#p-value: 0.124 

This looks at the difference between the overall averages for between subjects group 1 and between subjects group 2. It's suggested that it's more satisfactory in terms of type I errors.

--

sppbi(outcome ~ between*within, id = id, data = data_long)

#Test statistics:
#                                                         Estimate
#within2-within1               between1-between2            5.345
#within2-within3               between1-between2            1.013
#within1-within3               between1-between2           -5.875

#Test whether the corresponding population parameters are the same:
#p-value: 0.572 

This output is a little confusing. Mair & Wilcox (2020) suggest we should see something more like this:

#Test Statistics:
#                                                         Estimate                                                         
#within1    between1-between2                              5.345
#within2    between1-between2                              1.013
#within3    between1-between2                             -5.875

#Test whether the corresponding population parameters are the same:
#p-value: 0.572 

This computes M-estimators based on measurement occasion differences for each group separately. This would be the null hypothesis: H0 : θ1,2|1 − θ1,3|1 = θ1,4|1 − θ2,3|1 = θ2,4|1 − θ3,4|1 = θ1,2|2 − θ1,3|2 = θ1,4|2 − θ2,3|2 = θ2,4|2 − θ3,4|2 = 0

--

The Rallfun-v40.txt version of all of these provides the same eventual statistics but in a format that is harder to interpret. This being said, it does provide a contrast table for the final analysis:

#     [,1] [,2] [,3]
#[1,]    1    0    0
#[2,]    0    1    0
#[3,]    0    0    1
#[4,]   -1    0    0
#[5,]    0   -1    0
#[6,]    0    0   -1 

Column-wise this tells us that it performs the following comparisons, before conducting a single hypothesis test comprising all of these:

• Contrast 1: b1w1 to b2w1
• Contrast 2: b1w2 to b2w2
• Contrast 3: b1w3 to b2w3

--

Regarding running planned contrasts, I'm not sure if this is an option yet? It would be great if there were some robust way of running e.g., a simple effects analysis, but I'm not sure if there is?

There are a number of ways to run robust multiple comparisons/linear contrasts using the functions in Rallfun-v40.txt though (using data in wide format).

• bwmcp(2,3, data_wide, tr = 0.2)

Described as 'A bootstrap-t for multiple comparisons among for all main effects and interactions in a between-by-within design. The analysis is done by generating bootstrap samples (default n=599) and using an appropriate linear contrast.'

It returns the following:

#$Fac.A
#     con.num    psihat       se      test crit.value    p.value
#[1,]       1 -14.02105 8.604797 -1.629446   1.893387 0.09015025

1 Contrast (between subjects factor): between1 vs between2

#$Fac.B
#     con.num   psihat       se      test crit.value   p.value
#[1,]       1 6.283111 3.783206 1.6607900    2.26226 0.0884808
#[2,]       2 6.952670 5.422252 1.2822477    2.26226 0.1936561
#[3,]       3 0.669559 4.519172 0.1481597    2.26226 0.8714524

3 contrasts (within subjects factor):

1. within1 vs within2
2. within1 vs within3
3. within2 vs within3

#$Fac.AB
#     con.num    psihat       se       test crit.value     p.value
#[1,]       1 -9.451155 3.783206 -2.4981866   2.295026 0.008347245
#[2,]       2 -4.473331 5.422252 -0.8249951   2.295026 0.429048414
#[3,]       3  4.977824 4.519172  1.1014903   2.295026 0.280467446

3 contrasts (interactions):

1. [b1w1 & b2w2] vs [b2w1 & b1w2]
2. [b1w1 & b2w3] vs [b2w1 & b1w3]
3. [b1w2 and b2w3] vs [b2w2 & b1w3]

• bwamcp(2,3, data_wide, tr = 0.2)

• Performs multiple comparisons associated with Factor A using the (bootstrap-t) method BWAMCP

• It returns three sets of results corresponding to Factor A, Factor B, and all interactions. The critical value reported for each of the three set of tests in designed to control the probability of at least one type I error. Note the adjusted p-values.

#$test
#     con.num       test     crit       se       df
#[1,]       1 -2.2886813 2.578012 4.070108 21.58990
#[2,]       2  0.0380117 2.576128 3.577192 21.79337
#[3,]       3 -1.1849082 2.575349 4.086265 21.87870

#$psihat
#     con.num     psihat   ci.lower ci.upper   p.value adj.p.value
#[1,]       1 -9.3151798 -19.807965 1.177605 0.0322647  0.09147072
#[2,]       2  0.1359751  -9.079328 9.351278 0.9700244  0.99997118
#[3,]       3 -4.8418487 -15.365408 5.681710 0.2487656  0.56465996

Contrasts (compares the between groups in every within groups condition):

1. b1w1 vs b2w1
2. b1w2 vs b2w2
3. b1w3 vs b2w3

• bwbmcp(2,3, data_wide, tr = 0.2, con = 0, alpha = 0.05, dif = T,pool=F)

• uses method BWBMCP to compare the levels of Factor B. • Gives all pairwise comparisons among levels of Factor B in a mixed design using trimmed means. The pool option allows you to pool dependent groups across Factor A for each level of Factor B.

• "Critical p-values” are reported in the column headed by p.crit. These indicate how small a p-value must be in order to reject it, given the goal that the family-wise error be equal to some specified α value.

#$TESTS.4.EACH.LEVEL.OF.A[[1]]
#     Group Group        test   p.value    p.crit       se
#[1,]     1     2 -0.35472450 0.7301628 0.02500000 2.135674
#[2,]     1     3 -0.45389406 0.6595966 0.01666667 2.882919
#[3,]     2     3  0.01232317 0.9904102 0.05000000 1.862900

#$TESTS.4.EACH.LEVEL.OF.A[[2]]
#     Group Group        test    p.value     p.crit     se
#[1,]     1     2  1.86769503 0.08641073 0.01666667 2.964151
#[2,]     1     3  0.86431908 0.40435907 0.02500000 3.395994
#[3,]     2     3 -0.08201215 0.93598912 0.05000000 3.054188


#$PSIHAT.4.EACH.LEVEL.OF.A
#$PSIHAT.4.EACH.LEVEL.OF.A[[1]]
#     Group Group      psihat  ci.lower ci.upper
#[1,]     1     2 -0.75757576 -6.887114 5.371963
#[2,]     1     3 -1.30853994 -9.582728 6.965648
#[3,]     2     3  0.02295684 -5.323702 5.369615

#$PSIHAT.4.EACH.LEVEL.OF.A[[2]]
#     Group Group     psihat  ci.lower  ci.upper
#[1,]     1     2  5.5361305 -2.702648 13.774909
#[2,]     1     3  2.9352227 -6.503852 12.374297
#[3,]     2     3 -0.2504805 -8.739514  8.238553

Contrasts: All within group conditions for b1

1. b1w1 vs b1w2
2. b1w1 vs b1w3
3. b1w2 vs b1w3 All within group conditions for b2
4. b2w1 vs b2w2
5. b2w1 vs b2w3
6. b2w2 vs b2w3

If pool is set to true [i.e., bwbmcp(2,3, data_wide, tr = 0.2, con = 0, alpha = 0.05, dif = T,pool=T) ] then it pools the data for you and then calls the function rmmcp, with the following output:

#$POOLED.RESULTS$test
#     Group Group       test   p.value     p.crit       se
#[1,]     1     2  1.3241169 0.1984734 0.01666667 1.827659
#[2,]     1     3  0.3905554 0.6997198 0.02500000 2.238924
#[3,]     2     3 -0.1246142 0.9019117 0.05000000 2.017569

#$POOLED.RESULTS$psihat
#     Group Group     psihat  ci.lower ci.upper
#[1,]     1     2  2.4200337 -2.299012 7.139080
#[2,]     1     3  0.8744241 -4.906517 6.655365
#[3,]     2     3 -0.2514177 -5.460815 4.957980 

Contrasts: Data pooled from b1 & b2

1. w1 vs w2
2. w1 vs w3
3. w2 vs w3

• bwimcp(2,3, data_wide, tr = 0.2)

• For the interaction, this compares trimmed means using a nonbootstrap method (the bootstrap version is spmcpi), with the following output:

#     A A B B     psihat    p.value     p.crit
#[1,] 1 2 1 2 -6.2937063 0.09074896 0.01666667
#[2,] 1 2 1 3 -4.2437626 0.33564988 0.02500000
#[3,] 1 2 2 3  0.2734374 0.93824609 0.05000000

Performs 3 contrasts:

1. b1w1 vs b2w1
2. b1w2 vs b2w2
3. b1w3 vs b2w3 So it compares the between subjects groups in each of the within subjects conditions

--

From what I can tell, it isn't possible to do planned contrasts or simple effects analysis, which is a shame for anyone with specific hypotheses to test, given the severity of the p-value corrections needed to account for multiple comparisons. Perhaps I'm incorrect here?

--

If you spot any mistakes in this response (or have anything to add to it) then please do point them out.

--

The sources for this are as follows:

• WRS2 stuff: Mair, P., & Wilcox, R. (2020). Robust statistical methods in R using the WRS2 package. Behavior research methods, 52(2), 464-488.

• Rallfun-v40.txt stuff: Wilcox, R. R. (2021). Introduction to robust estimation and hypothesis testing (Fifth Edition). Academic press.

The former is aimed at researchers and so is a bit easier to follow. It also contains some handy supplementary materials, relying upon datasets already present in the WRS2 package.