SO SÁNH T-TEST THÍ NGHIỆM CHO TÔM ĂN

LOAD DỮ LIỆU

library(readxl)
shrimp_data <- read_excel("compare-hai-trung-binh_V2222.xlsx")
shrimp_data$treatment <- as.factor(shrimp_data$treatment)
print(shrimp_data, n = 40)
## # A tibble: 40 × 2
##    treatment    FC
##    <fct>     <dbl>
##  1 diet_1    0.539
##  2 diet_1    0.682
##  3 diet_1    0.548
##  4 diet_1    0.771
##  5 diet_1    0.771
##  6 diet_1    0.441
##  7 diet_1    0.521
##  8 diet_1    0.244
##  9 diet_1    0.369
## 10 diet_1    0.45 
## 11 diet_1    0.441
## 12 diet_1    0.664
## 13 diet_1    0.253
## 14 diet_1    0.477
## 15 diet_1    0.459
## 16 diet_1    0.387
## 17 diet_1    0.441
## 18 diet_1    0.45 
## 19 diet_1    0.539
## 20 diet_1    0.414
## 21 diet_6    0.271
## 22 diet_6    0.578
## 23 diet_6    0.754
## 24 diet_6    0.64 
## 25 diet_6    0.815
## 26 diet_6    0.254
## 27 diet_6    0.227
## 28 diet_6    0.438
## 29 diet_6    0.499
## 30 diet_6    0.368
## 31 diet_6    0.245
## 32 diet_6    0.227
## 33 diet_6    0.684
## 34 diet_6    0.517
## 35 diet_6    0.403
## 36 diet_6    0.324
## 37 diet_6    0.429
## 38 diet_6    0.517
## 39 diet_6    0.218
## 40 diet_6    0.543

XỬ LÝ THỐNG KÊ

CÁCH 1: SỬ DỤNG PACKAGE Hmisc

# library("Hmisc")
# Hmisc::describe(shrimp_data)

CÁCH 2: SỬ DỤNG PACKAGE psych

options(width = 200)
library("psych")
psych::describe(shrimp_data ~ treatment)
## 
##  Descriptive statistics by group 
## treatment: diet_1
##            vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## treatment*    1 20 1.00 0.00   1.00    1.00 0.00 1.00 1.00  0.00  NaN      NaN 0.00
## FC            2 20 0.49 0.14   0.45    0.49 0.11 0.24 0.77  0.53 0.36    -0.47 0.03
## ------------------------------------------------------------------------------------------------------------------------------------------------------ 
## treatment: diet_6
##            vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## treatment*    1 20 2.00 0.00   2.00    2.00 0.00 2.00 2.00   0.0  NaN      NaN 0.00
## FC            2 20 0.45 0.18   0.43    0.43 0.23 0.22 0.81   0.6 0.36    -1.09 0.04

CÁCH 3: SỬ DỤNG PACKAGE pastecs

# library("pastecs")
# pastecs::stat.desc(shrimp_data)

CÁCH 4: SỬ DỤNG PACKAGE doBy

# library("doBy")# 
# doBy::summaryBy(FC ~ treatment, data = shrimp_data,
#                 FUN = function(x) { c("Mean" = mean(x), 
#                                       "Std Dev" = sd(x)) })

CÁCH 5: SỬ DỤNG PACKAGE parameters

# library("parameters")
# if (require("lme4", quietly = TRUE)) {
#   model <- lmer(Petal.Length ~ Sepal.Length + (1 | Species), data = iris)
#   p_value_satterthwaite(model)
# }

CÁCH 6: SỬ DỤNG PACKAGE easystats

library("tidyverse")
library("easystats")
shrimp_data %>%
  group_by(treatment) %>%
  report() %>%
  summary()
## The data contains 40 observations, grouped by treatment, of the following 2 variables:
## 
## - diet_1 (n = 20):
##   - FC: Mean = 0.49, SD = 0.14, range: [0.24, 0.77]
## 
## - diet_6 (n = 20):
##   - FC: Mean = 0.45, SD = 0.18, range: [0.22, 0.81]
report(t.test(shrimp_data$FC ~ shrimp_data$treatment))
## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of shrimp_data$FC by shrimp_data$treatment (mean in group diet_1 = 0.49, mean in group diet_6 = 0.45) suggests that the effect is positive,
## statistically not significant, and small (difference = 0.05, 95% CI [-0.06, 0.15], t(35.90) = 0.87, p = 0.390; Cohen's d = 0.28, 95% CI [-0.35, 0.90])

CÁCH 7: SỬ DỤNG PACKAGE packHV

library("packHV")
packHV::desc(as.data.frame(shrimp_data),
     vars = c("FC", "treatment"),
     group = "treatment")
## mean (sd) and t-test/anova for numeric variables
## N (%) and khi-2 or Fisher exact test for categorical variables
##                                                                             
##  Covariate            Whole sample treatment diet_1 treatment diet_6 p-value
##                             (N=40)           (N=20)           (N=20)        
##         FC             0.47 (0.16)      0.49 (0.14)      0.45 (0.18) 0.39   
##  treatment diet_1         20 (50%)        20 (100%)           0 (0%) <0.001 
##            diet_6         20 (50%)           0 (0%)        20 (100%)

VẼ ĐỒ THỊ HISTOGRAM BY GROUP WITH BOXPLOT

library("packHV")
subset(shrimp_data, treatment == "diet_1") -> diet_1
subset(shrimp_data, treatment == "diet_6") -> diet_6

par(mfrow = c(1, 2))
hist_boxplot(diet_1$FC, 
             col = "lightblue",
             freq = FALSE,
             density = TRUE)
hist_boxplot(diet_6$FC, 
             col = "lightgreen",
             freq = FALSE,
             density = TRUE)

CÁCH 8: SỬ DỤNG PACKAGE ggplot2

VẼ ĐỒ THỊ QQ PLOT

library("ggplot2")
shrimp_data %>%
  ggplot(aes(sample = FC)) +
  geom_qq() + geom_qq_line() +
  facet_wrap(~treatment, scales = "free_y")

CÁCH 9: SỬ DỤNG PACKAGE ggpubr

library("ggpubr")
library("ggplot2")

# compare_means(residual ~ treatment, data = shrimp_data, method = "wilcox.test")

compare_means(FC ~ treatment, data = shrimp_data, method = "t.test")
## # A tibble: 1 × 8
##   .y.   group1 group2     p p.adj p.format p.signif method
##   <chr> <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr> 
## 1 FC    diet_1 diet_6 0.390  0.39 0.39     ns       T-test
p <- ggboxplot(data = shrimp_data, x = "treatment", y = "FC",
               color = "treatment", palette = "jco",
               add = "jitter")
#  Add p-value
# p + stat_compare_means(method = "wilcox.test")
# Change method
p + stat_compare_means(method = "t.test")