CODEWITHSUNDEEP
r
user16841302
user16841302

Reputation: 11

TukeyHSD() function in R

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.

Upvotes: 0

Views: 763

Answers (2)

Paul Schmidt
Paul Schmidt

Reputation: 1229

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)

Upvotes: 1

12666727b9
12666727b9

Reputation: 1139

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

Upvotes: 1

Related Questions