TukeyHSD() function in R

Question

I am trying to run TukeyHSD for a data set (40 observations). The outcome is 'mbvs' and the two factors are sex (0 or 1) and treatment (0 or 1) ?

I am trying to find out what is the pairwise difference in the mean outcome between treatment and control groups in men (sex=0)?

Could someone help me with the code for this ? Thank you.

Sure. Can you post a small dataframe from your original data? — 12666727b9, Dec 08 '21 at 20:00
It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. — MrFlick, Dec 08 '21 at 20:00
So for example I tried, m1=aov(mbvs~sex+tmt==0,data=.....) however that only gives 13 observations but in reality it must be 20 observations: sex=0 with tmt = 0 or 1 — user16841302, Dec 08 '21 at 20:09
That means that, if I am not wrong, that sex should be equal == 0 not tmt — 12666727b9, Dec 08 '21 at 20:16
I tried it that way as well (only sex==0) but it still gives 13 observations — user16841302, Dec 08 '21 at 20:20
Don't post a picture of your data. Use `dput(yourdataframe)` and paste the results into your question. — dcarlson, Dec 08 '21 at 20:37

12666727b9 · Accepted Answer · 2021-12-08T20:46:14.600

Look, I do not whether this answer can help you, but I've tried to semplify the dataset you provided with a smallest one, just to make you understand the logic presiding over it.

library(multcomp)
library(dplyr) 

mbvs <- c(6.4, 6.2, 6.9, 6.9, 5.4, 8.4)
sex <- c(0,0,0,0,0,1)
tmt <- c(0,1,0,0,0,1) 


   
#let's make this as your dataset 

df1 <- data.frame(mbvs, sex, tmt, stringsAsFactors = FALSE)

postHocs = df1 %>% 
  mutate(sex = factor(sex, levels = c('0', '1')), 
         tmt = factor(tmt, levels = c('0', '1'))) %>% 
  filter(sex == '0') %>%
  aov(mbvs ~ tmt, data = .) %>% 
  glht(., linfct = mcp(tmt = 'Tukey'))

The final result I've obtained is

     General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Linear Hypotheses:
           Estimate
1 - 0 == 0     -0.2

Alternatively, with the TukeyHSD function you can run this code

postHocs = df1 %>% 
  mutate(sex = factor(sex, levels = c('0', '1')), 
         tmt = factor(tmt, levels = c('0', '1'))) %>% 
  filter(sex == '0') %>%
  aov(mbvs ~ tmt, data = .) %>% 
  TukeyHSD(., conf.level=.95)

  Tukey multiple comparisons of means
    95% family-wise confidence level
Fit: aov(formula = mbvs ~ tmt, data = .)

$tmt
    diff       lwr      upr     p adj
1-0 -0.2 -2.715945 2.315945 0.8166266

score 1 · Answer 2 · answered Dec 10 '21 at 09:24

here is a chapter I wrote where Tukey tests are conducted using yet another function compared to those in @mały_statystyczny's answer. As you can see, the example data PlantGrowth has weight ~ group instead of your mbvs ~ tmt.

library(emmeans)
library(multcomp)
library(multcompView)

# set up model
model <- lm(weight ~ group, data = PlantGrowth)

# get (adjusted) weight means per group
model_means <- emmeans(object = model,
                       specs = "group")

# show differences
pairs(model_means, adjust = "Tukey", alpha = 0.05)
#>  contrast    estimate    SE df t.ratio p.value
#>  ctrl - trt1    0.371 0.279 27   1.331  0.3909
#>  ctrl - trt2   -0.494 0.279 27  -1.772  0.1980
#>  trt1 - trt2   -0.865 0.279 27  -3.103  0.0120
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates

# add letters to each mean
cld(object = model_means,
    adjust = "Tukey",
    Letters = letters,
    alpha = 0.05)
#>  group emmean    SE df lower.CL upper.CL .group
#>  trt1    4.66 0.197 27     4.16     5.16  a    
#>  ctrl    5.03 0.197 27     4.53     5.53  ab   
#>  trt2    5.53 0.197 27     5.02     6.03   b   
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 3 estimates 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
#> significance level used: alpha = 0.05 
#> NOTE: Compact letter displays can be misleading
#>       because they show NON-findings rather than findings.
#>       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

^{Created on 2021-12-10 by the reprex package (v2.0.1)}

TukeyHSD() function in R

2 Answers2