Gọi các packages liên quan
library(psych)
library(beeswarm)
library(readr)
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
Load dữ liệu
Litho <- "https://raw.githubusercontent.com/ngocdlu/Litho/master/Litho.csv"
df <- read_csv(Litho)
## Rows: 57 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): species
## dbl (5): len, wid, rat, cir, pet
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
attach(df)
Kiểm tra và mô tả dữ liệu
Sử dụng hàm View và hàm Summary () để show dữ liệu và summary từng biến trong tập dữ liệu:
summary(df)
## len wid rat cir
## Min. : 6.29 Min. :2.520 Min. :2.080 Min. :0.4700
## 1st Qu.:11.19 1st Qu.:4.000 1st Qu.:2.330 1st Qu.:0.6000
## Median :13.04 Median :4.800 Median :2.610 Median :0.6600
## Mean :14.78 Mean :5.495 Mean :2.712 Mean :0.6547
## 3rd Qu.:19.40 3rd Qu.:7.490 3rd Qu.:2.990 3rd Qu.:0.7200
## Max. :26.46 Max. :9.880 Max. :4.270 Max. :0.7800
## pet species
## Min. :0.840 Length:57
## 1st Qu.:1.100 Class :character
## Median :1.230 Mode :character
## Mean :1.546
## 3rd Qu.:2.070
## Max. :3.400
Mô tả dữ liệu theo biến phân nhóm
describe.by(len, species)
## Warning: describe.by is deprecated. Please use the describeBy function
##
## Descriptive statistics by group
## group: L_cadamensis
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 19 21.38 2.92 21.14 21.39 2.79 16.12 26.46 10.34 0.19 -0.89 0.67
## ------------------------------------------------------------
## group: L_campylolepis
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 19 12.74 2.9 12.93 12.68 2.49 7.05 19.4 12.35 0.17 -0.06 0.66
## ------------------------------------------------------------
## group: L_gougerotae
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 19 10.23 1.72 10.93 10.29 2.05 6.29 13.18 6.89 -0.4 -0.65 0.4
Thực hiện so sánh 2 nhóm độc lập và kiểm định t
Giả định trong kiểm định t
- Dữ liệu tuân theo luật phân phối chuẩn (normoal distribution)
- Hai nhóm độc lập
- Cỡ mẫu đủ lớn
- Mẫu được thu thập ngẫu nhiễn (simple random samples)
Tạo dataset mới có tên là “df1” gồm 2 nhóm
- Sau khi mô tả dữ liệu chúng ta thấy dataset df biến phân nhóm “species” có 3 nhóm, vì vậy chúng ta sẽ tách thành 1 dataset mới bao gồm 2 nhóm “L_cadamensis” và “L_campylolepis”.
- 2 nhóm “L_cadamensis” và “L_campylolepis” có thứ tự từ 1 đến 38 nên chúng ta thực hiện lệnh tách dataframe như sau.
df1 <- df[(1:38),]
Kiểm tra dữ liệu trong df1
attach(df1)
## The following objects are masked from df:
##
## cir, len, pet, rat, species, wid
describeBy(len, species)
##
## Descriptive statistics by group
## group: L_cadamensis
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 19 21.38 2.92 21.14 21.39 2.79 16.12 26.46 10.34 0.19 -0.89 0.67
## ------------------------------------------------------------
## group: L_campylolepis
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 19 12.74 2.9 12.93 12.68 2.49 7.05 19.4 12.35 0.17 -0.06 0.66
- Nhận dạng dữ liệu có tuân theo luật phân bố chuẩn hay không
# Nhận dạng bằng biểu đồ histogram hoặc density plot
hist(len)
plot(density(len), col=2)
# Kiểm tra phân bố chuẩn bằng qqnorm
qqnorm(len)
# Thêm đường kỳ giá trị vọng
qqline(len, col=3)
# Kiểm tra bằng qplot trong ggplot2 đối với từng nhóm
qplot(len, data = df1, geom = "density",
color = species, linetype = species)
- 3 biểu đồ trên cho thấy dữ liệu biến “len” ko tuân theo luật phân bố chuẩn.
Thử vẽ biểu đồ và thực hiện t.test xem thế nào.
- So sánh chiều dài phiến lá của 2 loài “L_cadamensis” và “L_campylolepis” bằng biểu đồ hộp
# Biểu diễn giá trị biến len bằng beeswarm
beeswarm(len~species, data =df1, col=c("Red", "Blue"))
# Vẽ biểu đồ hộp
boxplot(len~species, add=T, data =df1)
- Thực hiện t.test cho 2 nhóm độc lập
t.test(len~species)
##
## Welch Two Sample t-test
##
## data: len by species
## t = 9.1514, df = 35.997, p-value = 6.29e-11
## alternative hypothesis: true difference in means between group L_cadamensis and group L_campylolepis is not equal to 0
## 95 percent confidence interval:
## 6.723995 10.552847
## sample estimates:
## mean in group L_cadamensis mean in group L_campylolepis
## 21.37737 12.73895
- Có giá trị p = 6.29e-11 < 0.05 suy ra sự khác biệt chiều dài phiến lá giữa 2 loài có ý nghĩa thống kê.
Phân tích phương sai trong so sánh nhiều hơn 2 nhóm
So sánh chiều dài cuống lá của 3 loaì thực vật trong df
attach(df)
## The following objects are masked from df1:
##
## cir, len, pet, rat, species, wid
## The following objects are masked from df (pos = 4):
##
## cir, len, pet, rat, species, wid
ggplot(data=df,aes(x=species, y=pet, fill=species)) + geom_boxplot()
Phân tích phuong sai và tích hậu định
av <-aov(pet~species)
summary(av)
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 17.055 8.528 114.5 <2e-16 ***
## Residuals 54 4.022 0.074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- Giá trị của p < 0.05 chứng tỏ có sự khác biệt có ý nghĩa thống kê giữa các nhóm. Tuy nhiên chưa xác định chính xác nhóm nào khac biệt có ý nghĩa kê nhóm nào, vì vậy chúng ta thực hiện phân tích hậu định bằng phương pháp TukeyHSD
Phân tích hậu định bằng hàm TukeyHSD
t<-TukeyHSD(av)
t
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pet ~ species)
##
## $species
## diff lwr upr p adj
## L_campylolepis-L_cadamensis -1.1068421 -1.3202456 -0.8934386 0.000000
## L_gougerotae-L_cadamensis -1.2073684 -1.4207719 -0.9939649 0.000000
## L_gougerotae-L_campylolepis -0.1005263 -0.3139298 0.1128772 0.496839
Biểu diễn kết quả phân tích hậh định bằng biểu đồ
plot(t, xlim=c(-5, 2))
- Nhìn vào biểu đồ chúng ta thấy, phép so sánh nào có khoảng tin cậy 95% nằm lệch về một phía (trái hoặc phải) của giá trị “0” thì sự khác biệt đó có ý nghĩa kê.