Phân Tích Dữ Liệu Sản Xuất Thực Phẩm và Khí Thải
Sinh viên Thống Kê
2025-10-05
1 Giới Thiệu
Bài phân tích này nghiên cứu dữ liệu về sản xuất thực phẩm và tác động môi trường của các loại thực phẩm khác nhau. Dữ liệu bao gồm thông tin về khí thải nhà kính, sử dụng đất, sử dụng nước và các chỉ số môi trường khác.
Mục tiêu: Áp dụng các phương pháp thống kê để phân tích và so sánh tác động môi trường của các nhóm thực phẩm.
# Cài đặt và load các package cần thiết
if(!require(tidyverse)) install.packages("tidyverse")
if(!require(knitr)) install.packages("knitr")
if(!require(ggplot2)) install.packages("ggplot2")
library(tidyverse)
library(knitr)
library(ggplot2)# Đọc dữ liệu
data <- read.csv("Food_Production.csv", stringsAsFactors = FALSE)
# Xem cấu trúc dữ liệu
str(data)## 'data.frame': 43 obs. of 23 variables:
## $ Food.product : chr "Wheat & Rye (Bread)" "Maize (Meal)" "Barley (Beer)" "Oatmeal" ...
## $ Land.use.change : num 0.1 0.3 0 0 0 0 0.6 1.2 0 0 ...
## $ Animal.Feed : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Farm : num 0.8 0.5 0.2 1.4 3.6 0.2 0.2 0.5 0.5 1.1 ...
## $ Processing : num 0.2 0.1 0.1 0 0.1 0 0 0 0.2 0 ...
## $ Transport : num 0.1 0.1 0 0.1 0.1 0.1 0.1 0.8 0.6 0.1 ...
## $ Packging : num 0.1 0.1 0.5 0.1 0.1 0 0 0.1 0.1 0.4 ...
## $ Retail : num 0.1 0 0.3 0 0.1 0 0 0 0 0 ...
## $ Total_emissions : num 1.4 1.1 1.1 1.6 4 0.3 0.9 2.6 1.4 1.6 ...
## $ Eutrophying.emissions.per.1000kcal..gPO.eq.per.1000kcal. : num NA NA NA 4.28 9.51 ...
## $ Eutrophying.emissions.per.kilogram..gPO.eq.per.kilogram. : num NA NA NA 11.2 35.1 ...
## $ Eutrophying.emissions.per.100g.protein..gPO.eq.per.100.grams.protein. : num NA NA NA 8.64 49.39 ...
## $ Freshwater.withdrawals.per.1000kcal..liters.per.1000kcal. : num NA NA NA 184 610 ...
## $ Freshwater.withdrawals.per.100g.protein..liters.per.100g.protein. : num NA NA NA 371 3167 ...
## $ Freshwater.withdrawals.per.kilogram..liters.per.kilogram. : num NA NA NA 482 2248 ...
## $ Greenhouse.gas.emissions.per.1000kcal..kgCO.eq.per.1000kcal. : num NA NA NA 0.945 1.207 ...
## $ Greenhouse.gas.emissions.per.100g.protein..kgCO.eq.per.100g.protein. : num NA NA NA 1.91 6.27 ...
## $ Land.use.per.1000kcal..m..per.1000kcal. : num NA NA NA 2.9 0.76 ...
## $ Land.use.per.kilogram..m..per.kilogram. : num NA NA NA 7.6 2.8 ...
## $ Land.use.per.100g.protein..m..per.100g.protein. : num NA NA NA 5.85 3.94 ...
## $ Scarcity.weighted.water.use.per.kilogram..liters.per.kilogram. : num NA NA NA 18786 49576 ...
## $ Scarcity.weighted.water.use.per.100g.protein..liters.per.100g.protein. : num NA NA NA 14451 69826 ...
## $ Scarcity.weighted.water.use.per.1000kcal..liters.per.1000.kilocalories.: num NA NA NA 7162 13450 ...# Làm sạch dữ liệu: chỉ giữ các cột số có đủ dữ liệu
# Chọn biến Total_emissions làm biến phân tích chính
data_clean <- data %>%
select(Food.product, Total_emissions, Land.use.change, Farm,
Processing, Transport, Packging, Retail) %>%
filter(!is.na(Total_emissions) & Total_emissions > 0)
# Tạo nhóm thực phẩm
data_clean <- data_clean %>%
mutate(Food_Group = case_when(
grepl("Beef|Lamb|Pig|Poultry|Meat", Food.product) ~ "Thịt",
grepl("Milk|Cheese|Eggs", Food.product) ~ "Sữa & Trứng",
grepl("Rice|Wheat|Maize|Barley|Oatmeal", Food.product) ~ "Ngũ cốc",
grepl("Oil|Palm|Sunflower|Rapeseed|Olive|Soybean", Food.product) ~ "Dầu thực vật",
grepl("Nuts|Groundnuts", Food.product) ~ "Hạt",
grepl("Fish|Shrimps", Food.product) ~ "Hải sản",
TRUE ~ "Khác"
))
head(data_clean)## Food.product Total_emissions Land.use.change Farm Processing Transport
## 1 Wheat & Rye (Bread) 1.4 0.1 0.8 0.2 0.1
## 2 Maize (Meal) 1.1 0.3 0.5 0.1 0.1
## 3 Barley (Beer) 1.1 0.0 0.2 0.1 0.0
## 4 Oatmeal 1.6 0.0 1.4 0.0 0.1
## 5 Rice 4.0 0.0 3.6 0.1 0.1
## 6 Potatoes 0.3 0.0 0.2 0.0 0.1
## Packging Retail Food_Group
## 1 0.1 0.1 Ngũ cốc
## 2 0.1 0.0 Ngũ cốc
## 3 0.5 0.3 Ngũ cốc
## 4 0.1 0.0 Ngũ cốc
## 5 0.1 0.1 Ngũ cốc
## 6 0.0 0.0 Khác2 Phần 1: Thống Kê Mô Tả
2.1 Bảng Tần Số
2.1.1 Bảng tần số nhóm thực phẩm
# Tạo bảng tần số cho nhóm thực phẩm
freq_table <- data_clean %>%
group_by(Food_Group) %>%
summarise(
Tan_so = n(),
Ti_le = round(n()/nrow(data_clean) * 100, 2)
) %>%
arrange(desc(Tan_so))
kable(freq_table, caption = "Bảng tần số các nhóm thực phẩm")| Food_Group | Tan_so | Ti_le |
|---|---|---|
| Khác | 21 | 48.84 |
| Dầu thực vật | 5 | 11.63 |
| Ngũ cốc | 5 | 11.63 |
| Thịt | 5 | 11.63 |
| Sữa & Trứng | 3 | 6.98 |
| Hạt | 2 | 4.65 |
| Hải sản | 2 | 4.65 |
Nhận xét: - Nhóm “Khác” chiếm tỉ lệ cao nhất, bao gồm các loại rau củ quả - Các nhóm thịt, ngũ cốc, và dầu thực vật cũng được đại diện tốt trong dữ liệu
2.1.2 Bảng tần số phân loại mức khí thải
# Phân loại mức khí thải
data_clean <- data_clean %>%
mutate(Emission_Level = cut(Total_emissions,
breaks = c(0, 2, 5, 10, Inf),
labels = c("Thấp", "Trung bình", "Cao", "Rất cao"),
include.lowest = TRUE))
# Bảng tần số
emission_freq <- data_clean %>%
group_by(Emission_Level) %>%
summarise(
Tan_so = n(),
Ti_le_phan_tram = round(n()/nrow(data_clean) * 100, 2)
)
kable(emission_freq, caption = "Phân loại mức khí thải")| Emission_Level | Tan_so | Ti_le_phan_tram |
|---|---|---|
| Thấp | 22 | 51.16 |
| Trung bình | 8 | 18.60 |
| Cao | 6 | 13.95 |
| Rất cao | 7 | 16.28 |
2.2 Biểu Đồ Mô Tả
2.2.1 Biểu đồ cột - Khí thải theo nhóm thực phẩm
# Tính trung bình khí thải theo nhóm
avg_emissions <- data_clean %>%
group_by(Food_Group) %>%
summarise(TB_Khi_thai = mean(Total_emissions, na.rm = TRUE)) %>%
arrange(desc(TB_Khi_thai))
ggplot(avg_emissions, aes(x = reorder(Food_Group, TB_Khi_thai),
y = TB_Khi_thai, fill = Food_Group)) +
geom_col() +
coord_flip() +
labs(title = "Trung bình khí thải theo nhóm thực phẩm",
x = "Nhóm thực phẩm",
y = "Khí thải trung bình (kgCO2eq)") +
theme_minimal() +
theme(legend.position = "none")
Nhận xét: - Nhóm thịt có mức khí thải cao nhất - Ngũ cốc và rau củ có mức khí thải thấp hơn nhiều
2.2.2 Biểu đồ tròn - Tỉ lệ các nhóm thực phẩm
ggplot(freq_table, aes(x = "", y = Tan_so, fill = Food_Group)) +
geom_bar(stat = "identity", width = 1) +
coord_polar("y", start = 0) +
labs(title = "Phân bố các nhóm thực phẩm") +
theme_void() +
theme(legend.position = "right")
2.2.3 Histogram - Phân bố khí thải
ggplot(data_clean, aes(x = Total_emissions)) +
geom_histogram(bins = 20, fill = "steelblue", color = "black", alpha = 0.7) +
labs(title = "Phân bố tổng khí thải",
x = "Tổng khí thải (kgCO2eq)",
y = "Tần số") +
theme_minimal()
Nhận xét: Phân bố lệch phải, đa số thực phẩm có mức khí thải thấp, một số ít có mức rất cao.
2.2.4 Boxplot - So sánh khí thải giữa các nhóm
ggplot(data_clean, aes(x = Food_Group, y = Total_emissions, fill = Food_Group)) +
geom_boxplot() +
coord_flip() +
labs(title = "So sánh khí thải giữa các nhóm thực phẩm",
x = "Nhóm thực phẩm",
y = "Tổng khí thải (kgCO2eq)") +
theme_minimal() +
theme(legend.position = "none")
Nhận xét: - Nhóm thịt có độ biến động lớn và giá trị ngoại lai cao - Ngũ cốc có phân bố đồng đều hơn
3 Phần 2: Các Đại Lượng Thống Kê Cơ Bản
# Tính các đại lượng thống kê cho Total_emissions
stats_summary <- data_clean %>%
summarise(
So_quan_sat = n(),
Trung_binh = round(mean(Total_emissions), 3),
Trung_vi = round(median(Total_emissions), 3),
Mode = names(sort(table(round(Total_emissions, 1)), decreasing = TRUE))[1],
Do_lech_chuan = round(sd(Total_emissions), 3),
Phuong_sai = round(var(Total_emissions), 3),
Min = round(min(Total_emissions), 3),
Q1 = round(quantile(Total_emissions, 0.25), 3),
Q3 = round(quantile(Total_emissions, 0.75), 3),
Max = round(max(Total_emissions), 3),
IQR = round(IQR(Total_emissions), 3),
He_so_bien_dong = round(sd(Total_emissions)/mean(Total_emissions) * 100, 2)
)
# Hiển thị bảng
kable(t(stats_summary), col.names = c("Giá trị"),
caption = "Thống kê mô tả cho tổng khí thải")| Giá trị | |
|---|---|
| So_quan_sat | 43 |
| Trung_binh | 5.972 |
| Trung_vi | 1.6 |
| Mode | 0.3 |
| Do_lech_chuan | 10.502 |
| Phuong_sai | 110.287 |
| Min | 0.2 |
| Q1 | 0.85 |
| Q3 | 6 |
| Max | 59.6 |
| IQR | 5.15 |
| He_so_bien_dong | 175.85 |
Giải thích: - Trung bình: Giá trị trung bình của khí thải - Trung vị: Giá trị ở giữa khi sắp xếp dữ liệu - Độ lệch chuẩn: Đo mức độ phân tán của dữ liệu - Phương sai: Bình phương của độ lệch chuẩn - IQR: Khoảng tứ phân vị (Q3 - Q1) - Hệ số biến động: (SD/Mean) × 100%, cho biết mức độ biến thiên tương đối
Nhận xét: Dữ liệu có độ biến động lớn do sự khác biệt giữa các loại thực phẩm.
4 Phần 3: Ước Lượng Điểm và Khoảng Tin Cậy
4.1 Ước lượng Trung bình
# Ước lượng khoảng tin cậy 95% cho trung bình
n <- nrow(data_clean)
mean_val <- mean(data_clean$Total_emissions)
sd_val <- sd(data_clean$Total_emissions)
se <- sd_val / sqrt(n)
# Khoảng tin cậy
ci_95 <- t.test(data_clean$Total_emissions, conf.level = 0.95)$conf.int
cat("Ước lượng điểm (Trung bình):", round(mean_val, 3), "\n")## Ước lượng điểm (Trung bình): 5.972## Độ lệch chuẩn mẫu: 10.502## Sai số chuẩn: 1.602## Khoảng tin cậy 95%: [ 2.74 , 9.204 ]Giải thích: - Với độ tin cậy 95%, trung bình tổng khí thải của tất cả các loại thực phẩm nằm trong khoảng trên - Sai số chuẩn đo độ chính xác của ước lượng
4.2 Ước lượng Tỉ lệ
# Tính tỉ lệ thực phẩm có khí thải cao (> 5 kgCO2eq)
high_emission <- sum(data_clean$Total_emissions > 5)
n_total <- nrow(data_clean)
p_hat <- high_emission / n_total
# Khoảng tin cậy cho tỉ lệ (phương pháp Wilson)
prop_test <- prop.test(high_emission, n_total, conf.level = 0.95)
cat("Số thực phẩm có khí thải cao:", high_emission, "\n")## Số thực phẩm có khí thải cao: 13## Tổng số thực phẩm: 43## Ước lượng tỉ lệ: 0.302cat("Khoảng tin cậy 95%: [", round(prop_test$conf.int[1], 3), ",",
round(prop_test$conf.int[2], 3), "]\n")## Khoảng tin cậy 95%: [ 0.177 , 0.463 ]Nhận xét: Tỉ lệ thực phẩm có khí thải cao trong tổng thể được ước lượng với độ tin cậy 95%.
4.3 Ước lượng Phương sai
# Khoảng tin cậy cho phương sai (dùng phân phối Chi-square)
alpha <- 0.05
df <- n - 1
var_val <- var(data_clean$Total_emissions)
chi_lower <- qchisq(1 - alpha/2, df)
chi_upper <- qchisq(alpha/2, df)
ci_var_lower <- (df * var_val) / chi_lower
ci_var_upper <- (df * var_val) / chi_upper
cat("Ước lượng điểm (Phương sai):", round(var_val, 3), "\n")## Ước lượng điểm (Phương sai): 110.287cat("Khoảng tin cậy 95% cho phương sai: [", round(ci_var_upper, 3), ",",
round(ci_var_lower, 3), "]\n")## Khoảng tin cậy 95% cho phương sai: [ 178.165 , 74.98 ]cat("Khoảng tin cậy 95% cho độ lệch chuẩn: [", round(sqrt(ci_var_upper), 3), ",",
round(sqrt(ci_var_lower), 3), "]\n")## Khoảng tin cậy 95% cho độ lệch chuẩn: [ 13.348 , 8.659 ]Giải thích: Phương sai đo mức độ phân tán của dữ liệu, khoảng tin cậy giúp ước lượng phương sai trong tổng thể.
5 Phần 4: Kiểm Định 1 Tổng Thể
5.1 Kiểm định Trung bình (t-test)
Giả thuyết: Kiểm tra xem trung bình khí thải có lớn hơn 5 kgCO2eq không?
- H₀: μ = 5 (trung bình bằng 5)
- H₁: μ > 5 (trung bình lớn hơn 5)
- Mức ý nghĩa: α = 0.05
# Kiểm định t một mẫu
t_test_result <- t.test(data_clean$Total_emissions, mu = 5, alternative = "greater")
print(t_test_result)##
## One Sample t-test
##
## data: data_clean$Total_emissions
## t = 0.60699, df = 42, p-value = 0.2736
## alternative hypothesis: true mean is greater than 5
## 95 percent confidence interval:
## 3.278442 Inf
## sample estimates:
## mean of x
## 5.972093##
## --- KẾT LUẬN ---## Giá trị p: 0.2736if(t_test_result$p.value < 0.05) {
cat("Kết luận: Bác bỏ H0. Có bằng chứng thống kê rằng trung bình khí thải > 5.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng để bác bỏ H0.\n")
}## Kết luận: Chưa đủ bằng chứng để bác bỏ H0.5.2 Kiểm định Tỉ lệ (Proportion test)
Giả thuyết: Kiểm tra xem tỉ lệ thực phẩm có khí thải cao (>5) có khác 0.5 không?
- H₀: p = 0.5
- H₁: p ≠ 0.5
- Mức ý nghĩa: α = 0.05
# Kiểm định tỉ lệ
prop_result <- prop.test(high_emission, n_total, p = 0.5, alternative = "two.sided")
print(prop_result)##
## 1-sample proportions test with continuity correction
##
## data: high_emission out of n_total, null probability 0.5
## X-squared = 5.9535, df = 1, p-value = 0.01469
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.1766968 0.4629894
## sample estimates:
## p
## 0.3023256##
## --- KẾT LUẬN ---## Giá trị p: 0.0147if(prop_result$p.value < 0.05) {
cat("Kết luận: Bác bỏ H0. Tỉ lệ thực phẩm có khí thải cao khác 0.5.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng để bác bỏ H0.\n")
}## Kết luận: Bác bỏ H0. Tỉ lệ thực phẩm có khí thải cao khác 0.5.6 Phần 5: Kiểm Định 2 Tổng Thể
6.1 So sánh Trung bình 2 nhóm (t-test độc lập)
Bài toán: So sánh khí thải trung bình giữa nhóm Thịt và Ngũ cốc
- H₀: μ_Thịt = μ_Ngũ_cốc
- H₁: μ_Thịt ≠ μ_Ngũ_cốc
- Mức ý nghĩa: α = 0.05
# Tách dữ liệu 2 nhóm
meat_data <- data_clean %>% filter(Food_Group == "Thịt")
grain_data <- data_clean %>% filter(Food_Group == "Ngũ cốc")
# Kiểm định t 2 mẫu độc lập
t_test_2sample <- t.test(meat_data$Total_emissions,
grain_data$Total_emissions,
alternative = "two.sided",
var.equal = FALSE)
print(t_test_2sample)##
## Welch Two Sample t-test
##
## data: meat_data$Total_emissions and grain_data$Total_emissions
## t = 2.2522, df = 4.0256, p-value = 0.08701
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.021444 48.741444
## sample estimates:
## mean of x mean of y
## 23.70 1.84##
## --- KẾT LUẬN ---## Trung bình Thịt: 23.7## Trung bình Ngũ cốc: 1.84## Giá trị p: 0.087if(t_test_2sample$p.value < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa thống kê về khí thải giữa Thịt và Ngũ cốc.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt.\n")
}## Kết luận: Chưa đủ bằng chứng về sự khác biệt.6.2 So sánh Tỉ lệ 2 nhóm
Bài toán: So sánh tỉ lệ thực phẩm có khí thải cao (>5) giữa nhóm Thịt và Ngũ cốc
# Đếm số lượng
meat_high <- sum(meat_data$Total_emissions > 5)
grain_high <- sum(grain_data$Total_emissions > 5)
meat_n <- nrow(meat_data)
grain_n <- nrow(grain_data)
# Kiểm định 2 tỉ lệ
prop_2sample <- prop.test(c(meat_high, grain_high), c(meat_n, grain_n))
print(prop_2sample)##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(meat_high, grain_high) out of c(meat_n, grain_n)
## X-squared = 6.4, df = 1, p-value = 0.01141
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.8 1.0
## sample estimates:
## prop 1 prop 2
## 1 0##
## --- KẾT LUẬN ---## Tỉ lệ Thịt: 1## Tỉ lệ Ngũ cốc: 0## Giá trị p: 0.0114if(prop_2sample$p.value < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa về tỉ lệ khí thải cao giữa 2 nhóm.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt.\n")
}## Kết luận: Có sự khác biệt có ý nghĩa về tỉ lệ khí thải cao giữa 2 nhóm.6.3 So sánh Phương sai 2 nhóm (F-test)
Bài toán: Kiểm tra xem phương sai của Thịt và Ngũ cốc có bằng nhau không?
- H₀: σ²_Thịt = σ²_Ngũ_cốc
- H₁: σ²_Thịt ≠ σ²_Ngũ_cốc
# F-test cho phương sai
var_test <- var.test(meat_data$Total_emissions, grain_data$Total_emissions)
print(var_test)##
## F test to compare two variances
##
## data: meat_data$Total_emissions and grain_data$Total_emissions
## F = 312.41, num df = 4, denom df = 4, p-value = 6.095e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 32.52755 3000.56888
## sample estimates:
## ratio of variances
## 312.4118##
## --- KẾT LUẬN ---## Phương sai Thịt: 469.555## Phương sai Ngũ cốc: 1.503## Giá trị p: 1e-04if(var_test$p.value < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa về phương sai giữa 2 nhóm.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt phương sai.\n")
}## Kết luận: Có sự khác biệt có ý nghĩa về phương sai giữa 2 nhóm.7 Phần 6: Phân Tích ANOVA và TukeyHSD
7.1 ANOVA một chiều
Bài toán: So sánh khí thải trung bình giữa TẤT CẢ các nhóm thực phẩm
- H₀: μ₁ = μ₂ = μ₃ = … (tất cả các nhóm có trung bình bằng nhau)
- H₁: Có ít nhất một cặp khác nhau
- Mức ý nghĩa: α = 0.05
# Phân tích ANOVA
anova_model <- aov(Total_emissions ~ Food_Group, data = data_clean)
anova_summary <- summary(anova_model)
print(anova_summary)## Df Sum Sq Mean Sq F value Pr(>F)
## Food_Group 6 1993 332.1 4.53 0.00162 **
## Residuals 36 2639 73.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## --- KẾT LUẬN ANOVA ---## Giá trị p: 0.001619if(p_value_anova < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa về khí thải giữa các nhóm thực phẩm.\n")
cat("=> Cần thực hiện kiểm định hậu nghiệm (Post-hoc) để tìm nhóm nào khác biệt.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt.\n")
}## Kết luận: Có sự khác biệt có ý nghĩa về khí thải giữa các nhóm thực phẩm.
## => Cần thực hiện kiểm định hậu nghiệm (Post-hoc) để tìm nhóm nào khác biệt.7.1.1 Biểu đồ ANOVA
ggplot(data_clean, aes(x = Food_Group, y = Total_emissions, fill = Food_Group)) +
geom_boxplot() +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "red") +
labs(title = "So sánh khí thải giữa các nhóm thực phẩm (ANOVA)",
subtitle = "Hình thoi đỏ = Trung bình",
x = "Nhóm thực phẩm",
y = "Tổng khí thải (kgCO2eq)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
7.2 Kiểm định TukeyHSD (Post-hoc)
Mục đích: Tìm cặp nhóm nào có sự khác biệt có ý nghĩa
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Total_emissions ~ Food_Group, data = data_clean)
##
## $Food_Group
## diff lwr upr p adj
## Hải sản-Dầu thực vật 3.0900000 -19.26712691 25.44713 0.9994386
## Hạt-Dầu thực vật -4.0600000 -26.41712691 18.29713 0.9973821
## Khác-Dầu thực vật -2.7742857 -16.07144884 10.52288 0.9944076
## Ngũ cốc-Dầu thực vật -3.5200000 -20.42039938 13.38040 0.9944592
## Sữa & Trứng-Dầu thực vật 4.1400000 -15.37490026 23.65490 0.9938822
## Thịt-Dầu thực vật 18.3400000 1.43960062 35.24040 0.0260969
## Hạt-Hải sản -7.1500000 -33.87187770 19.57188 0.9794674
## Khác-Hải sản -5.8642857 -25.63881912 13.91025 0.9658783
## Ngũ cốc-Hải sản -6.6100000 -28.96712691 15.74713 0.9663778
## Sữa & Trứng-Hải sản 1.0500000 -23.34362533 25.44363 0.9999994
## Thịt-Hải sản 15.2500000 -7.10712691 37.60713 0.3584015
## Khác-Hạt 1.2857143 -18.48881912 21.06025 0.9999932
## Ngũ cốc-Hạt 0.5400000 -21.81712691 22.89713 1.0000000
## Sữa & Trứng-Hạt 8.2000000 -16.19362533 32.59363 0.9385449
## Thịt-Hạt 22.4000000 0.04287309 44.75713 0.0492959
## Ngũ cốc-Khác -0.7457143 -14.04287741 12.55145 0.9999972
## Sữa & Trứng-Khác 6.9142857 -9.57881527 23.40739 0.8435934
## Thịt-Khác 21.1142857 7.81712259 34.41145 0.0003209
## Sữa & Trứng-Ngũ cốc 7.6600000 -11.85490026 27.17490 0.8797418
## Thịt-Ngũ cốc 21.8600000 4.95960062 38.76040 0.0046164
## Thịt-Sữa & Trứng 14.2000000 -5.31490026 33.71490 0.2854710# Trích xuất các cặp có p < 0.05
tukey_df <- as.data.frame(tukey_result$Food_Group)
tukey_df$Comparison <- rownames(tukey_df)
significant_pairs <- tukey_df %>%
filter(`p adj` < 0.05) %>%
arrange(`p adj`)
cat("\n--- CÁC CẶP CÓ SỰ KHÁC BIỆT CÓ Ý NGHĨA (p < 0.05) ---\n")##
## --- CÁC CẶP CÓ SỰ KHÁC BIỆT CÓ Ý NGHĨA (p < 0.05) ---| diff | lwr | upr | p adj | Comparison | |
|---|---|---|---|---|---|
| Thịt-Khác | 21.1143 | 7.8171 | 34.4114 | 0.0003 | Thịt-Khác |
| Thịt-Ngũ cốc | 21.8600 | 4.9596 | 38.7604 | 0.0046 | Thịt-Ngũ cốc |
| Thịt-Dầu thực vật | 18.3400 | 1.4396 | 35.2404 | 0.0261 | Thịt-Dầu thực vật |
| Thịt-Hạt | 22.4000 | 0.0429 | 44.7571 | 0.0493 | Thịt-Hạt |
8 Phần 7: Kiểm Định Phi Tham Số
8.1 Kiểm định Chi-square (χ²) - Độc lập
Bài toán: Kiểm tra mối liên hệ giữa nhóm thực phẩm và mức khí thải
- H₀: Nhóm thực phẩm và mức khí thải độc lập (không liên quan)
- H₁: Có mối liên hệ giữa nhóm thực phẩm và mức khí thải
# Tạo bảng chéo
contingency_table <- table(data_clean$Food_Group, data_clean$Emission_Level)
cat("Bảng chéo:\n")## Bảng chéo:##
## Thấp Trung bình Cao Rất cao
## Dầu thực vật 0 2 3 0
## Hải sản 0 0 1 1
## Hạt 1 1 0 0
## Khác 17 2 0 2
## Ngũ cốc 4 1 0 0
## Sữa & Trứng 0 2 0 1
## Thịt 0 0 2 3##
## Pearson's Chi-squared test
##
## data: contingency_table
## X-squared = 46.946, df = 18, p-value = 0.0002156##
## --- KẾT LUẬN ---## Giá trị χ²: 46.946## Bậc tự do: 18## Giá trị p: 0.000216if(chi_test$p.value < 0.05) {
cat("Kết luận: Có mối liên hệ có ý nghĩa giữa nhóm thực phẩm và mức khí thải.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về mối liên hệ.\n")
}## Kết luận: Có mối liên hệ có ý nghĩa giữa nhóm thực phẩm và mức khí thải.8.1.1 Biểu đồ Chi-square
# Biểu đồ cột chồng
contingency_df <- as.data.frame(contingency_table)
names(contingency_df) <- c("Food_Group", "Emission_Level", "Count")
ggplot(contingency_df, aes(x = Food_Group, y = Count, fill = Emission_Level)) +
geom_bar(stat = "identity", position = "fill") +
labs(title = "Phân bố mức khí thải theo nhóm thực phẩm",
x = "Nhóm thực phẩm",
y = "Tỉ lệ",
fill = "Mức khí thải") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(labels = scales::percent)
8.2 Kiểm định Goodness-of-Fit (χ² GOF)
Bài toán: Kiểm tra xem các nhóm thực phẩm có phân bố đều không?
- H₀: Các nhóm thực phẩm phân bố đều
- H₁: Các nhóm thực phẩm không phân bố đều
# Đếm số lượng từng nhóm
observed <- table(data_clean$Food_Group)
n_groups <- length(observed)
expected_prob <- rep(1/n_groups, n_groups)
# Kiểm định GOF
gof_test <- chisq.test(observed, p = expected_prob)
print(gof_test)##
## Chi-squared test for given probabilities
##
## data: observed
## X-squared = 43.767, df = 6, p-value = 8.22e-08##
## --- KẾT LUẬN ---## Giá trị p: 0if(gof_test$p.value < 0.05) {
cat("Kết luận: Bác bỏ H0. Các nhóm thực phẩm không phân bố đều.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng để bác bỏ H0.\n")
}## Kết luận: Bác bỏ H0. Các nhóm thực phẩm không phân bố đều.# So sánh quan sát vs kỳ vọng
comparison <- data.frame(
Nhom = names(observed),
Quan_sat = as.numeric(observed),
Ky_vong = round(gof_test$expected, 2)
)
kable(comparison, caption = "So sánh tần số quan sát và kỳ vọng")| Nhom | Quan_sat | Ky_vong | |
|---|---|---|---|
| Dầu thực vật | Dầu thực vật | 5 | 6.14 |
| Hải sản | Hải sản | 2 | 6.14 |
| Hạt | Hạt | 2 | 6.14 |
| Khác | Khác | 21 | 6.14 |
| Ngũ cốc | Ngũ cốc | 5 | 6.14 |
| Sữa & Trứng | Sữa & Trứng | 3 | 6.14 |
| Thịt | Thịt | 5 | 6.14 |
8.3 Kiểm định Wilcoxon (Mann-Whitney U)
Bài toán: So sánh phân bố khí thải giữa Thịt và Ngũ cốc (không giả định phân phối chuẩn)
- H₀: Phân bố của 2 nhóm giống nhau
- H₁: Phân bố của 2 nhóm khác nhau
# Kiểm định Wilcoxon (Mann-Whitney)
wilcox_test <- wilcox.test(meat_data$Total_emissions,
grain_data$Total_emissions,
alternative = "two.sided")
print(wilcox_test)##
## Wilcoxon rank sum test with continuity correction
##
## data: meat_data$Total_emissions and grain_data$Total_emissions
## W = 25, p-value = 0.01193
## alternative hypothesis: true location shift is not equal to 0##
## --- KẾT LUẬN ---## Giá trị p: 0.011925if(wilcox_test$p.value < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa về phân bố khí thải giữa Thịt và Ngũ cốc.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt.\n")
}## Kết luận: Có sự khác biệt có ý nghĩa về phân bố khí thải giữa Thịt và Ngũ cốc.##
## Median Thịt: 21.1## Median Ngũ cốc: 1.4Lưu ý: Kiểm định Wilcoxon phù hợp khi dữ liệu không phân phối chuẩn hoặc có ngoại lai.
8.3.1 Kiểm định Kruskal-Wallis (ANOVA phi tham số)
Bài toán: So sánh nhiều nhóm khi không đáp ứng giả định ANOVA
# Kruskal-Wallis test
kruskal_test <- kruskal.test(Total_emissions ~ Food_Group, data = data_clean)
print(kruskal_test)##
## Kruskal-Wallis rank sum test
##
## data: Total_emissions by Food_Group
## Kruskal-Wallis chi-squared = 24.095, df = 6, p-value = 0.0005018##
## --- KẾT LUẬN ---## Giá trị p: 0.000502if(kruskal_test$p.value < 0.05) {
cat("Kết luận: Có sự khác biệt có ý nghĩa về khí thải giữa các nhóm thực phẩm.\n")
} else {
cat("Kết luận: Chưa đủ bằng chứng về sự khác biệt.\n")
}## Kết luận: Có sự khác biệt có ý nghĩa về khí thải giữa các nhóm thực phẩm.9 Tổng Kết
9.1 Kết quả chính
- Thống kê mô tả:
- Nhóm Thịt có khí thải cao nhất
- Phân bố khí thải lệch phải với nhiều ngoại lai
- Ước lượng:
- Trung bình khí thải được ước lượng với khoảng tin cậy 95%
- Tỉ lệ thực phẩm có khí thải cao được xác định
- Kiểm định:
- Có sự khác biệt có ý nghĩa về khí thải giữa các nhóm thực phẩm
- Nhóm Thịt và Ngũ cốc khác biệt rõ rệt cả về trung bình và phân bố
- ANOVA và TukeyHSD chỉ ra các cặp nhóm có sự khác biệt
- Kiểm định phi tham số:
- Kết quả nhất quán với các kiểm định tham số
- Có mối liên hệ giữa loại thực phẩm và mức khí thải
9.2 Khuyến nghị
Dựa trên phân tích, để giảm khí thải: - Ưu tiên tiêu thụ thực phẩm từ nhóm Ngũ cốc, Rau củ - Hạn chế thực phẩm từ nhóm Thịt, đặc biệt là thịt bò - Cân nhắc các nguồn protein thay thế (đậu, hạt)
10 Phụ Lục
10.1 Kiểm tra giả định phân phối chuẩn
# Shapiro-Wilk test (chỉ với n < 5000)
if(nrow(data_clean) <= 5000) {
shapiro_test <- shapiro.test(data_clean$Total_emissions)
cat("Shapiro-Wilk test:\n")
print(shapiro_test)
cat("\n")
}## Shapiro-Wilk test:
##
## Shapiro-Wilk normality test
##
## data: data_clean$Total_emissions
## W = 0.55601, p-value = 3.335e-10# Q-Q plot
qqnorm(data_clean$Total_emissions, main = "Q-Q Plot - Kiểm tra phân phối chuẩn")
qqline(data_clean$Total_emissions, col = "red")
Nhận xét: Nếu p < 0.05 trong Shapiro test hoặc Q-Q plot không thẳng hàng => dữ liệu không phân phối chuẩn, nên dùng kiểm định phi tham số.
10.2 Thông tin phiên làm việc
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Asia/Saigon
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.50 lubridate_1.9.4 forcats_1.0.1 stringr_1.5.2
## [5] dplyr_1.1.4 purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
## [9] tibble_3.3.0 ggplot2_4.0.0 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.5.1 tidyselect_1.2.1
## [5] jquerylib_0.1.4 scales_1.4.0 yaml_2.3.10 fastmap_1.2.0
## [9] R6_2.6.1 labeling_0.4.3 generics_0.1.4 bslib_0.9.0
## [13] pillar_1.11.1 RColorBrewer_1.1-3 tzdb_0.5.0 rlang_1.1.6
## [17] stringi_1.8.7 cachem_1.1.0 xfun_0.53 sass_0.4.10
## [21] S7_0.2.0 timechange_0.3.0 cli_3.6.5 withr_3.0.2
## [25] magrittr_2.0.4 digest_0.6.37 grid_4.5.1 rstudioapi_0.17.1
## [29] hms_1.1.3 lifecycle_1.0.4 vctrs_0.6.5 evaluate_1.0.5
## [33] glue_1.8.0 farver_2.1.2 rmarkdown_2.30 tools_4.5.1
## [37] pkgconfig_2.0.3 htmltools_0.5.8.1HẾT
Bài phân tích này được thực hiện với mục đích học tập môn Xác suất Thống kê.