Can I use multiple T tests as an alternative to post-hoc testing after AN(C)OVA? If so, how can I do this in R?

Question

I have a dataset where I have measured a phenotype across multiple years and multiple varieties of a crop. The first way I have searched for a relationship is a two-way ANOVA, where I ask "Is there an effect of variety on the phenotype, having accounted for the sampling structure?" i.e.

phenotype ~ sampling event + variety

The question I am really interested in, is if any varieties are significantly better or worse than the overall population mean.

I have been struggling to find the appropriate post-hoc test for this, given I get a significant relationship between phenotype and variety. TukeyHSD hasn't been working for me, as I am not interested in comparisons between the different varieties, but rather comparison to the mean.

As a work-around I have been using T tests with a p value adjusted for multiple testing. However, I am not sure whether this is a) valid, and b) possible to do in a more streamlined way in R. At the moment I have to manually create a new column where the variety I am interested in is "a" and everything else is "b", so the T test compares sample "a" to the mean. However, given I have 30 different varieties, and 8 different phenotypes, this will take a long time to do manually!

Any help is more than welcome! I have attached some example data. I have only included one phenotype, for simplicity

variety <- c(A,A,A,A,B,B,B,B,C,C,C,C,D,D,D,D,E,E,E,E,F,F,F,F,G,G,G,G,H,H,H,H)
sampling_event <- c(A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2)
phenotype1 <- c(13.86,14.48,15.5,16.22,14.72,14.72,16.6,16.98,16.98,12.34,15.6,17.82,17.6,9.26,13.46,12.24,13.1,16.22,15.94,10.86,12.44,10.58,17.3,13.38,15.2,13.66,18.2,14.9,15.68,18.8,15.94,13.38)

dummydata <-data.frame(variety, sampling_event, phenotype1)

I followed the instructions on https://www.datanovia.com/en/blog/how-to-perform-multiple-t-test-in-r-for-different-variables/. From this, I managed to generate T tests for each variety across all of the phenotypes, which is good, but I would rather have T tests for all the varieties for each phenotype.

Answer 1

I'm not sure if this is the best stat for your particular question but have you considered using estimated marginal means to look for differences between groups?

library(emmeans)

dummydata <- data.frame(
  variety = c('A','A','A','A',
              'B','B','B','B',
              'C','C','C','C',
              'D','D','D','D',
              'E','E','E','E',
              'F','F','F','F',
              'G','G','G','G',
              'H','H','H','H'),
  sampling_event = c('A1','A1','A2','A2','A1','A1','A2','A2','A1','A1','A2','A2',
                     'A1','A1','A2','A2','A1','A1','A2','A2','A1','A1','A2','A2','A1',
                     'A1','A2','A2','A1','A1','A2','A2'),
  phenotype1 <- c(13.86,14.48,15.5,16.22,14.72,14.72,16.6,16.98,16.98,12.34,15.6,17.82,17.6,9.26,13.46,12.24,13.1,16.22,15.94,10.86,12.44,10.58,17.3,13.38,15.2,13.66,18.2,14.9,15.68,18.8,15.94,13.38)
)

# Create linear model (setup based on what was asked in post)
lm_mod <- lm(phenotype1 ~ sampling_event + variety, data = dummydata)

# Estimated marginal means
emm_lm_mod <- emmeans(lm_mod, pairwise ~ sampling_event | variety)
emm_lm_mod
#> $emmeans
#> variety = A:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               14.6 1.23 23     12.0     17.1
#>  A2               15.5 1.23 23     12.9     18.0
#> 
#> variety = B:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               15.3 1.23 23     12.8     17.8
#>  A2               16.2 1.23 23     13.7     18.8
#> 
#> variety = C:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               15.2 1.23 23     12.7     17.8
#>  A2               16.1 1.23 23     13.6     18.7
#> 
#> variety = D:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               12.7 1.23 23     10.1     15.2
#>  A2               13.6 1.23 23     11.1     16.1
#> 
#> variety = E:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               13.6 1.23 23     11.0     16.1
#>  A2               14.5 1.23 23     11.9     17.0
#> 
#> variety = F:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               13.0 1.23 23     10.4     15.5
#>  A2               13.9 1.23 23     11.3     16.4
#> 
#> variety = G:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               15.0 1.23 23     12.5     17.6
#>  A2               15.9 1.23 23     13.4     18.5
#> 
#> variety = H:
#>  sampling_event emmean   SE df lower.CL upper.CL
#>  A1               15.5 1.23 23     12.9     18.0
#>  A2               16.4 1.23 23     13.9     19.0
#> 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#> variety = A:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = B:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = C:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = D:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = E:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = F:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = G:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745
#> 
#> variety = H:
#>  contrast estimate   SE df t.ratio p.value
#>  A1 - A2    -0.917 0.82 23  -1.120  0.2745

^{Created on 2024-05-27 with reprex v2.1.0}