# Show a random sample
set.seed(1234)
dplyr::sample_n(my_data, 10)
# Show the levels
levels(my_data$group)
[1] "ctrl" "trt1" "trt2"
my_data$group <- ordered(my_data$group,levels = c("ctrl", "trt1", "trt2"))
group_by(my_data, group) %>%
  summarise(
    count = n(),
    mean = mean(weight, na.rm = TRUE),
    sd = sd(weight, na.rm = TRUE)
  )
# Box plots
# ++++++++++++++++++++
# Plot weight by group and color by group
ggboxplot(my_data, x = "group", y = "weight", 
          color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment")

# Mean plots
# ++++++++++++++++++++
# Plot weight by group
# Add error bars: mean_se
# (other values include: mean_sd, mean_ci, median_iqr, ....)
#library("ggpubr")
ggline(my_data, x = "group", y = "weight", 
       add = c("mean_se", "jitter"), 
       order = c("ctrl", "trt1", "trt2"),
       ylab = "Weight", xlab = "Treatment")

# Box plot
boxplot(weight ~ group, data = my_data,
        xlab = "Treatment", ylab = "Weight",
        frame = FALSE, col = c("#00AFBB", "#E7B800", "#FC4E07"))

# plotmeans
#library("gplots")
plotmeans(weight ~ group, data = my_data, #frame = FALSE,
          xlab = "Treatment", ylab = "Weight",
          main="Mean Plot with 95% CI") 

# Compute the analysis of variance
res.aov <- aov(weight ~ group, data = my_data)
# Summary of the analysis
summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(res.aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ group, data = my_data)

$`group`
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064
#diff: difference between means of the two groups
#lwr, upr: the lower and the upper end point of the confidence interval at 95% (default)
#p adj: p-value after adjustment for the multiple comparisons.

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = weight ~ group, data = my_data)

Linear Hypotheses:
                 Estimate Std. Error t value Pr(>|t|)  
trt1 - ctrl == 0  -0.3710     0.2788  -1.331   0.3909  
trt2 - ctrl == 0   0.4940     0.2788   1.772   0.1979  
trt2 - trt1 == 0   0.8650     0.2788   3.103   0.0121 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
pairwise.t.test(my_data$weight, my_data$group,
                 p.adjust.method = "BH")

    Pairwise comparisons using t tests with pooled SD 

data:  my_data$weight and my_data$group 

     ctrl  trt1 
trt1 0.194 -    
trt2 0.132 0.013

P value adjustment method: BH 
# 1. Homogeneity of variances
plot(res.aov, 1)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.1192 0.3412
      27               
oneway.test(weight ~ group, data = my_data)

    One-way analysis of means (not assuming equal variances)

data:  weight and group
F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739
pairwise.t.test(my_data$weight, my_data$group,
                 p.adjust.method = "BH", pool.sd = FALSE)

    Pairwise comparisons using t tests with non-pooled SD 

data:  my_data$weight and my_data$group 

     ctrl  trt1 
trt1 0.250 -    
trt2 0.072 0.028

P value adjustment method: BH 
# 2. Normality
plot(res.aov, 2)

# Extract the residuals
aov_residuals <- residuals(object = res.aov )
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )

    Shapiro-Wilk normality test

data:  aov_residuals
W = 0.96607, p-value = 0.4379
kruskal.test(weight ~ group, data = my_data)

    Kruskal-Wallis rank sum test

data:  weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842
LS0tDQp0aXRsZTogIk9uZS1XYXkgQU5PVkEgVGVzdCBpbiBSIChBTk9WQSBOTyBSKSBwYXJhIG8gbWV1IGFtaWdvIE5JTE8gZGEgVUVSSiAtIDIwMTkiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQoNCmBgYHtyIGVjaG89RkFMU0UgLHdhcm1pbmc9Rn0NCmxpYnJhcnkoImRwbHlyIikNCmxpYnJhcnkoImdncHViciIpDQpsaWJyYXJ5KCJncGxvdHMiKQ0KIyBJZiAudHh0IHRhYiBmaWxlLCB1c2UgdGhpcw0KI215X2RhdGEgPC0gcmVhZC5kZWxpbShmaWxlLmNob29zZSgpKQ0KIyBPciwgaWYgLmNzdiBmaWxlLCB1c2UgdGhpcw0KI215X2RhdGEgPC0gcmVhZC5jc3YoZmlsZS5jaG9vc2UoKSkNCg0KbXlfZGF0YSA8LSBQbGFudEdyb3d0aA0KYGBgDQoNCg0KYGBge3J9DQojIFNob3cgYSByYW5kb20gc2FtcGxlDQpzZXQuc2VlZCgxMjM0KQ0KZHBseXI6OnNhbXBsZV9uKG15X2RhdGEsIDEwKQ0KDQojIFNob3cgdGhlIGxldmVscw0KbGV2ZWxzKG15X2RhdGEkZ3JvdXApDQoNCm15X2RhdGEkZ3JvdXAgPC0gb3JkZXJlZChteV9kYXRhJGdyb3VwLGxldmVscyA9IGMoImN0cmwiLCAidHJ0MSIsICJ0cnQyIikpDQpgYGANCg0KDQpgYGB7cn0NCmdyb3VwX2J5KG15X2RhdGEsIGdyb3VwKSAlPiUNCiAgc3VtbWFyaXNlKA0KICAgIGNvdW50ID0gbigpLA0KICAgIG1lYW4gPSBtZWFuKHdlaWdodCwgbmEucm0gPSBUUlVFKSwNCiAgICBzZCA9IHNkKHdlaWdodCwgbmEucm0gPSBUUlVFKQ0KICApDQpgYGANCg0KDQoNCmBgYHtyfQ0KIyBCb3ggcGxvdHMNCiMgKysrKysrKysrKysrKysrKysrKysNCiMgUGxvdCB3ZWlnaHQgYnkgZ3JvdXAgYW5kIGNvbG9yIGJ5IGdyb3VwDQpnZ2JveHBsb3QobXlfZGF0YSwgeCA9ICJncm91cCIsIHkgPSAid2VpZ2h0IiwgDQogICAgICAgICAgY29sb3IgPSAiZ3JvdXAiLCBwYWxldHRlID0gYygiIzAwQUZCQiIsICIjRTdCODAwIiwgIiNGQzRFMDciKSwNCiAgICAgICAgICBvcmRlciA9IGMoImN0cmwiLCAidHJ0MSIsICJ0cnQyIiksDQogICAgICAgICAgeWxhYiA9ICJXZWlnaHQiLCB4bGFiID0gIlRyZWF0bWVudCIpDQpgYGANCg0KDQpgYGB7cn0NCiMgTWVhbiBwbG90cw0KIyArKysrKysrKysrKysrKysrKysrKw0KIyBQbG90IHdlaWdodCBieSBncm91cA0KIyBBZGQgZXJyb3IgYmFyczogbWVhbl9zZQ0KIyAob3RoZXIgdmFsdWVzIGluY2x1ZGU6IG1lYW5fc2QsIG1lYW5fY2ksIG1lZGlhbl9pcXIsIC4uLi4pDQojbGlicmFyeSgiZ2dwdWJyIikNCmdnbGluZShteV9kYXRhLCB4ID0gImdyb3VwIiwgeSA9ICJ3ZWlnaHQiLCANCiAgICAgICBhZGQgPSBjKCJtZWFuX3NlIiwgImppdHRlciIpLCANCiAgICAgICBvcmRlciA9IGMoImN0cmwiLCAidHJ0MSIsICJ0cnQyIiksDQogICAgICAgeWxhYiA9ICJXZWlnaHQiLCB4bGFiID0gIlRyZWF0bWVudCIpDQpgYGANCg0KDQpgYGB7cn0NCiMgQm94IHBsb3QNCmJveHBsb3Qod2VpZ2h0IH4gZ3JvdXAsIGRhdGEgPSBteV9kYXRhLA0KICAgICAgICB4bGFiID0gIlRyZWF0bWVudCIsIHlsYWIgPSAiV2VpZ2h0IiwNCiAgICAgICAgZnJhbWUgPSBGQUxTRSwgY29sID0gYygiIzAwQUZCQiIsICIjRTdCODAwIiwgIiNGQzRFMDciKSkNCiMgcGxvdG1lYW5zDQojbGlicmFyeSgiZ3Bsb3RzIikNCnBsb3RtZWFucyh3ZWlnaHQgfiBncm91cCwgZGF0YSA9IG15X2RhdGEsICNmcmFtZSA9IEZBTFNFLA0KICAgICAgICAgIHhsYWIgPSAiVHJlYXRtZW50IiwgeWxhYiA9ICJXZWlnaHQiLA0KICAgICAgICAgIG1haW49Ik1lYW4gUGxvdCB3aXRoIDk1JSBDSSIpIA0KYGBgDQoNCg0KYGBge3J9DQojIENvbXB1dGUgdGhlIGFuYWx5c2lzIG9mIHZhcmlhbmNlDQpyZXMuYW92IDwtIGFvdih3ZWlnaHQgfiBncm91cCwgZGF0YSA9IG15X2RhdGEpDQojIFN1bW1hcnkgb2YgdGhlIGFuYWx5c2lzDQpzdW1tYXJ5KHJlcy5hb3YpDQpgYGANCg0KDQpgYGB7cn0NClR1a2V5SFNEKHJlcy5hb3YpDQoNCiNkaWZmOiBkaWZmZXJlbmNlIGJldHdlZW4gbWVhbnMgb2YgdGhlIHR3byBncm91cHMNCiNsd3IsIHVwcjogdGhlIGxvd2VyIGFuZCB0aGUgdXBwZXIgZW5kIHBvaW50IG9mIHRoZSBjb25maWRlbmNlIGludGVydmFsIGF0IDk1JSAoZGVmYXVsdCkNCiNwIGFkajogcC12YWx1ZSBhZnRlciBhZGp1c3RtZW50IGZvciB0aGUgbXVsdGlwbGUgY29tcGFyaXNvbnMuDQoNCmBgYA0KYGBge3IgZWNobz1GQUxTRX0NCmxpYnJhcnkobXVsdGNvbXApDQpzdW1tYXJ5KGdsaHQocmVzLmFvdiwgbGluZmN0ID0gbWNwKGdyb3VwID0gIlR1a2V5IikpKQ0KYGBgDQoNCmBgYHtyfQ0KcGFpcndpc2UudC50ZXN0KG15X2RhdGEkd2VpZ2h0LCBteV9kYXRhJGdyb3VwLA0KICAgICAgICAgICAgICAgICBwLmFkanVzdC5tZXRob2QgPSAiQkgiKQ0KYGBgDQoNCg0KDQpgYGB7cn0NCiMgMS4gSG9tb2dlbmVpdHkgb2YgdmFyaWFuY2VzDQpwbG90KHJlcy5hb3YsIDEpDQpgYGANCg0KDQpgYGB7ciBlY2hvPUZBTFNFfQ0KbGlicmFyeShjYXIpDQpsZXZlbmVUZXN0KHdlaWdodCB+IGdyb3VwLCBkYXRhID0gbXlfZGF0YSkNCmBgYA0KYGBge3J9DQpvbmV3YXkudGVzdCh3ZWlnaHQgfiBncm91cCwgZGF0YSA9IG15X2RhdGEpDQpgYGANCg0KDQpgYGB7cn0NCnBhaXJ3aXNlLnQudGVzdChteV9kYXRhJHdlaWdodCwgbXlfZGF0YSRncm91cCwNCiAgICAgICAgICAgICAgICAgcC5hZGp1c3QubWV0aG9kID0gIkJIIiwgcG9vbC5zZCA9IEZBTFNFKQ0KYGBgDQoNCmBgYHtyfQ0KIyAyLiBOb3JtYWxpdHkNCnBsb3QocmVzLmFvdiwgMikNCmBgYA0KDQoNCmBgYHtyfQ0KIyBFeHRyYWN0IHRoZSByZXNpZHVhbHMNCmFvdl9yZXNpZHVhbHMgPC0gcmVzaWR1YWxzKG9iamVjdCA9IHJlcy5hb3YgKQ0KIyBSdW4gU2hhcGlyby1XaWxrIHRlc3QNCnNoYXBpcm8udGVzdCh4ID0gYW92X3Jlc2lkdWFscyApDQpgYGANCg0KYGBge3J9DQprcnVza2FsLnRlc3Qod2VpZ2h0IH4gZ3JvdXAsIGRhdGEgPSBteV9kYXRhKQ0KYGBgDQoNCg==