Tiểu luận

File Pdf

https://drive.google.com/file/d/1LU3Z7TgUWQYyDcP3BC7nfjZcBdPVzFuY/view?usp=sharing

File dữ liệu excel

https://docs.google.com/spreadsheets/d/1f_ZV6ymniLs4GNraO2HCuUA_B9W-ufss/edit?usp=sharing&ouid=111976047388761728103&rtpof=true&sd=true

1 PHẦN MỞ ĐẦU

1.1 Lý do chọn đề tài

Sữa là một trong những thực phẩm quan trọng và phổ biến trong chế độ dinh dưỡng của con người, đặc biệt là trẻ em. Sữa cung cấp nhiều chất dinh dưỡng cần thiết cho sự phát triển và bảo vệ sức khỏe của cơ thể. Tuy nhiên, chất lượng sữa có thể bị ảnh hưởng bởi nhiều yếu tố khác nhau, từ quá trình sản xuất, chế biến, bảo quản đến tiêu thụ. Việc đánh giá và kiểm soát chất lượng sữa là rất cần thiết để đảm bảo an toàn và hiệu quả cho người tiêu dùng.

Có nhiều yếu tố có thể ảnh hưởng đến chất lượng sữa, từ giai đoạn nuôi bò sữa, thu hoạch, vận chuyển, chế biến, bảo quản đến tiêu thụ. Nếu không được kiểm soát và quản lý tốt, các yếu tố này có thể làm giảm chất lượng sữa hoặc gây ra các nguy cơ an toàn thực phẩm cho người tiêu dùng. Ví dụ, nếu bò sữa bị nhiễm khuẩn hoặc sử dụng thuốc kháng sinh quá liều, sữa có thể chứa các vi sinh vật gây bệnh hoặc các dư lượng thuốc kháng sinh. Nếu sữa không được vận chuyển và bảo quản ở nhiệt độ thích hợp, sữa có thể bị ôi thiu hoặc chua. Nếu sữa không được chế biến theo quy trình tiêu chuẩn, sữa có thể bị pha trộn hoặc giảm hàm lượng dinh dưỡng.

Với mong muốn góp phần nâng cao chất lượng cuộc sống của người. Tôi đã nghĩ đến đè tài nghiên cứu Phân tích các yếu tố ảnh hưởng đến chất lượng sữa để góp phần tạo nên nguồn năng lượng từ sữa cho con người và phần nào giải quyết các vấn đề thực tiễn trong ngành công nghiệp sữa. Tôi hy vọng rằng qua đề tài này, tôi có thể đóng góp vào việc nâng cao nhận thức và kiến thức của người tiêu dùng về chất lượng sữa. Tôi cũng mong muốn rằng đề tài này sẽ là một nguồn thông tin hữu ích cho các nhà sản xuất, chế biến và quản lý sữa trong việc cải thiện và duy trì chất lượng sữa.

1.2 Mục tiêu nghiên cứu

Đánh giá các yếu tố có thể ảnh hưởng đến chất lượng sữa. Từ đó, giúp người nông dân chăn nuôi đơn lẻ hoặc các trang trại bò sữa có thể hiểu thêm về các yếu tố làm ảnh hưởng đến chất lượng nguồn sữa. Giúp họ có thể nhận ra các điểm yếu làm ảnh hưởng đến chất lượng sữa cũng như việc làm như thế nào để tạo ra nguồn sữa nhiều và chất lượng cung cấp ra thị trường.

1.3 Đối tượng và phạm vi nghiên cứu

• Đối tượng:

  • Chất lượng sữa.

• Phạm vi nghiên cứu:

  • Để thực hiện đề tài này, tôi đã thu thập được một bộ dữ liệu từ https://www.kaggle.com.
  • Bộ dữ liệu gồm có 1059 mẫu quan sát, có 7 biến là các yếu tố ảnh hưởng đến chất lượng sữa, bao gồm: pH, Temprature, Odor, Fat, Turbidity, Colour và Taste. Biến Grade chất lượng của sữa chính là biến phụ thuộc mà chúng ta cần phân tích .

1.4 Phương pháp nghiên cứu

Phương pháp phân tích dữ liệu định tính.

Tôi đã sử dụng các phương pháp thống kê toán học để mô tả và suy luận để khám phá và phân tích các mối quan hệ giữa các biến này. Tôi cũng đã xây dựng các mô hình hồi quy và phân loại để dự báo và phân loại chất lượng sữa dựa trên các biến đầu vào.

Phần mềm R & RStudio là một công cụ mới nhưng khá được ưa dùng trong thống kê và xử lý dữ liệu - sẽ vì các tính năng hỗ trợ phân tích cũng như vẽ biểu đồ nhanh chóng. Trong bài phân tích này R & RStudio sẽ hỗ trợ chính cho việc nghiên cứu, đánh giá.

1.5 Kết cấu đề tài

• Bài nghiên cứu gồm có 4 chương:

  • Chương 1. Tổng quan lý thuyết.
  • Chương 2. Phương pháp nghiên cứu
  • Chương 3. Kết quả mô hình nghiên cứu.
  • Chương 4. Kết luận và hạn chế.

2 CHƯƠNG 1. TỔNG QUAN LÝ THUYẾT

2.1 Ngành công nghiệp sữa trong nước

Ngành công nghiệp sữa Việt Nam là một trong những ngành có vai trò quan trọng trong nền kinh tế và đời sống của người dân. Ngành này không chỉ cung cấp các sản phẩm dinh dưỡng, chăm sóc sức khỏe cho hàng triệu người tiêu dùng, mà còn góp phần phát triển nông nghiệp, tạo việc làm, tăng thu nhập cho các hộ chăn nuôi bò sữa, và đóng góp vào ngân sách nhà nước. Ngành sữa Việt Nam cũng đã có những bước tiến vượt bậc trong việc nâng cao chất lượng sản phẩm, đa dạng hóa mẫu mã và chủng loại, áp dụng công nghệ hiện đại, và mở rộng thị trường xuất khẩu.

Ngành sữa Việt Nam bắt đầu phát triển từ những năm 1920, khi các trang trại bò sữa nhỏ lẻ được thành lập để cung cấp sữa tươi cho thị trường địa phương. Sau đó, vào những năm 1970, các doanh nghiệp sản xuất sữa được hình thành và phát triển, trong đó có Công ty Cổ phần Sữa Việt Nam (Vinamilk) là doanh nghiệp tiên phong và lớn nhất trong ngành. Đến nay, ngành sữa Việt Nam đã có hơn 200 doanh nghiệp hoạt động với quy mô lớn nhỏ khác nhau, trong đó có các doanh nghiệp trong nước như Vinamilk, TH True Milk, Nutifood, Mộc Châu Milk, và các doanh nghiệp nước ngoài như FrieslandCampina (Dutch Lady), Nestlé, Abbott, Mead Johnson.

Theo Euromonitor, thị trường sữa Việt Nam đạt giá trị 135.000 tỷ đồng vào năm 2020, tăng hơn 8% so với năm 2019. Các sản phẩm chủ lực của ngành bao gồm sữa bột, sữa thanh trùng, sữa chua, pho mát, bơ và các sản phẩm từ sữa khác. Trong đó, sữa thanh trùng là phân khúc có giá trị lớn nhất và tăng trưởng cao nhất trong ngành.

2.2 Thuận lợi và khó khăn

2.2.1 Thuận lợi

Sự tăng trưởng của thị trường được thúc đẩy bởi các yếu tố sau:

• Nhu cầu tiêu dùng sữa ngày càng tăng do dân số tăng, thu nhập người dân tăng và nhận thức về chăm sóc sức khỏe của người tiêu dùng cải thiện.

• Chính sách hỗ trợ của nhà nước cho ngành sữa thông qua các chương trình như “Sữa học đường”, “Sản xuất bò sữa”, “Hỗ trợ giá sữa cho trẻ em” và các ưu đãi thuế cho các doanh nghiệp sản xuất sữa.

• Sự đầu tư của các doanh nghiệp trong và ngoài nước vào ngành sữa thông qua việc xây dựng các nhà máy sản xuất hiện đại, mở rộng chuỗi cung ứng bò sữa, phát triển các sản phẩm mới theo xu hướng tiêu dùng, và tham gia các hiệp định thương mại tự do để mở rộng thị trường xuất khẩu.

2.2.2 Khó khăn

Tuy nhiên, ngành sữa Việt Nam cũng đang đối mặt với một số thách thức và cơ hội trong tương lai, như:

• Sự cạnh tranh khốc liệt từ các doanh nghiệp nước ngoài có uy tín và chất lượng cao, đặc biệt là trong phân khúc sữa bột cho trẻ em.

• Sự thiếu hụt nguồn cung sữa tươi trong nước do năng suất bò sữa thấp, chi phí chăn nuôi cao và thiếu hụt đất canh tác. Hiện nay, ngành sữa Việt Nam phải nhập khẩu khoảng 70% nguyên liệu sữa từ các nước như New Zealand, Mỹ, Australia và châu Âu.

• Sự thay đổi của thị hiếu và nhu cầu của người tiêu dùng, đòi hỏi các doanh nghiệp phải liên tục cải tiến và đổi mới sản phẩm để phù hợp với xu hướng tiêu dùng. Một số xu hướng tiêu dùng hiện nay bao gồm các sản phẩm sữa hữu cơ, sữa không lactose, sữa chức năng và sữa thực vật.

Ngành công nghiệp sữa Việt Nam là một ngành có tiềm năng lớn để phát triển trong bối cảnh kinh tế toàn cầu hội nhập và hợp tác. Để duy trì và nâng cao vị thế của mình trên thị trường trong và ngoài nước, các doanh nghiệp ngành sữa cần phải tận dụng các cơ hội, vượt qua các thách thức, và không ngừng nâng cao năng lực sản xuất, chất lượng sản phẩm và hiệu quả kinh doanh.

3 CHƯƠNG 2. PHƯƠNG PHÁP NGHIÊN CỨU

3.1 Các package sử dụng trong bài

library(ggplot2)

library(readxl)

library(epitools)

library(caret)

## Loading required package: lattice

library(DescTools)

## 
## Attaching package: 'DescTools'

## The following objects are masked from 'package:caret':
## 
##     MAE, RMSE

3.2 Dữ liệu nghiên cứu

3.2.1 Mô tả bộ dữ liệu

Bộ dữ liệu này được sử dụng để phân tích và khám phá các mối quan hệ giữa các yếu tố ảnh hưởng đến chất lượng sữa, cũng như để xây dựng các mô hình dự báo và phân loại chất lượng sữa. Bộ dữ liệu này cũng có thể giúp cho các nhà sản xuất, chế biến và tiêu thụ sữa nâng cao hiểu biết và kiểm soát chất lượng sữa một cách tốt hơn.

Bộ dữ liệu này gồm có 1059 mẫu quan sát, có 7 biến là: pH, Temprature, Taste, Odor, Fat, Turbidity, Colour và 1 biến Grade là biến định tính phụ thuộc vào các biến này.

3.2.2 Mô tả các biến

• Biến định tính
  • Taste: Đây là biến định tính có nghĩa là hương vị của sữa. Sữa có hương vị khi nếm vào không ngon sẽ được đánh số (0), sữa có hương vị khi nếm vào ngon sẽ được đánh số (1).
  • Odor: Đây là biến định tính có nghĩa là mùi của sữa. Sữa có mùi khi ngửi không thơm sẽ được đánh số (0), sữa có mùi khi ngửi thơm sẽ được đánh số (1).
  • Fat: Đây là biến định tính có nghĩa là độ béo của sữa. Sữa có độ béo thấp sẽ được đánh số (0), sữa có độ béo cao sẽ được đánh số (1).
  • Turbidity: Đây là biến định tính có nghĩa là độ đục của sữa. Sữa có độ đục thấp sẽ được đánh số (0), sữa có độ đục cao sẽ được đánh số (1).
    
  • Grade: Đây là biến định tính có nghĩa là chất lượng của sữa. Sữa có chất lượng thấp sẽ đánh số (0), sữa có chất lượng cao sẽ đánh số (1).
• Biến định lượng:
  • pH: Đây là biến định lượng có nghĩa là độ pH của sữa. Độ pH của sữa dao động từ 3 đến 9,5, độ pH phù hợp nhất nằm ở khoảng 6,25 đến 6,90.
  • Temprature: Đây là biến định lượng có nghĩa là nhiệt độ của sữa. Nhiệt độ của sữa dao động từ 34 đến 90, nhiệt độ phù hợp nhất nằm ở khoảng 34 đến 45,20.
  • Colour: Đây là biến định lượng có nghĩa là màu của sữa. Màu của sữa dao động từ 240 đến 255, phân màu phù hợp nhất nằm ở khoảng 245 đến 255.

3.2.3 Dữ liệu

milknew <- read_excel("C:/Users/ASUS/Downloads/milknew.xlsx")
View(milknew)
library(DT)
milknew %>% DT::datatable(milknew)

3.3 Xác định biến phụ thuộc

3.3.1 Biến định tính

Biến định tính Grade sẽ được chọn là biến định tính phụ thuộc bởi vì:

• Mục đích của bài là phân tích các yếu tố ảnh hưởng đến chất lượng sữa. Chính vì lẽ đó nên ta sẽ quan sát xem các biến pH, Temprature, Taste, Odor, Fat, Turbidity và Colour sẽ ảnh hưởng đến biến Grade thế nào.
• Xét trên phương diện thực tế, chất lượng sữa cũng sẽ được đánh giá dựa trên các tiêu chí nồng độ, hương vị, độ béo, màu sắc,…

3.3.2 Biến định lượng

Biến định lượng Colour sẽ được chọn làm biến phụ thuộc bởi vì:

• Màu sữa phần nào sẽ thể hiện được chất lượng của sữa và nó có thể bị ảnh hưởng bởi các yếu tố được nêu trong bộ dữ liệu.

3.4 Mô hình nghiên cứu

3.4.1 Mô hình tuyến tính

GLM (General Linear Model) là một mô hình tuyến tính tổng quát giải thích mối quan hệ giữa một hàm kỳ vọng có điều kiện với các biến giải thích thông qua một phương trình dự đoán có dạng tuyến tính.

GLM phổ biến vì nó có thể xử lý nhiều loại dữ liệu với các kiểu biến phản ứng khác nhua như nhị thức, Poisson hay đa thức.

Để mô hình hóa ảnh hưởng của X là sử dụng dạng hồi quy tuyến tính cổ điển, theo đó giá trị kỳ vọng có điều kiện của Y là một hàm tuyến tính của X.

\(\pi(x) = \beta_0 + \beta_1.x\)

Đây là GLM với thành phần ngẫu nhiên nhị thức và hàm liên kết đồng nhất. Mô hình này được gọi là mô hình xác suất tuyến tính hay LPM (Linear Probability Model). \(\beta_1\) biểu thị cho sự thay đổi xác suất cho mỗi đơn vị thay đổi của x.

3.4.2 Mô hình hồi quy Logistic

Để khắc phục một số hạn chế của LPM, người ta đề xuất một hàm liên kết mà nó đồng biến với xác suất \(\pi\) và giá trị biến thiên của nó là tùy ý đó là \(log(\frac{\pi(x)}{1 - \pi(x)})\), được gọi là hàm logistic hay logit(\(\pi\)).

Mô hình hồi quy logit:

\(log(\frac{\pi(x)}{1 - \pi(x)}) = \beta_0 + \beta_1.x_1 + \beta_2.x_2 + ... +\beta_m.x_m\)

Trong đó \(\beta_0, \beta_1,..., \beta_m\) là hệ số hồi quy.

3.4.3 Mô hình Probit

Mô hình Probit được sử dung để ước lượng mô hình có biến phụ thuộc dạng nhị phân. Mô hình Probit là một GLM với thành phần ngẫu nhiên nhị phân và liên kết Probit. Hàm hôi quy Probit có dạng:

\(\pi(x) = \phi(\alpha + \beta x)\)

Trong đó, \(\phi(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{1}{2}t^2}dt\) là hàm phân phối chuẩn chính tắc. Với \(Probit(\pi(x)) = \phi^{-1}(\pi(x))\), ta có mô hình Probit:

\(Probit(\pi(x)) = \alpha + \beta x\)

3.4.4 Mô hình Cloglog

Trong trường hợp tỷ lệ cược phụ thuộc vào các biến giải thích dưới dạng hàm lũy thừa:

\(Odds(x)=\ \frac{\pi(x)}{1-\pi(x)}=\lambda.x_1^{\beta_1}.x_2^{\beta_2}....x_m^{\beta_m}\)

Ta nhận được mô hình logistic có các biến dự báo ở dạng log, gọi là mô hình logistic dạng log-log.

\(\log{\left(\frac{\pi(x)}{1-\pi(x)}\right)=\beta_0+\beta_1.logx_1+\beta_2.logx_2+...+\beta_m.logx_m}\)

3.4.5 Phân phối nhị thức

Phân phối nhị thức là một phân phối xác suất rời rạc với số lần xảy ra của một sự kiện thành công trong một số lần thử cố định, các tham số là n và p (n là số tự nhiên > 0, 0 < p < 1).

Nếu gọi Y là số lần thành công trong n lần lặp lại phép thử, thì theo công thức Becnoulli, xác suất để có k lần thành công là:

\(p_k = P(Y = k) = C\binom{n}{k}.p^k.(1 - p)^(n - k), k = 1,2,...,n\)

Như vậy mô hình lặp lại \(n\) lần một phép thử đã cho ra một biến ngẫu nhiên có phân phối Nhị thức \(B(n,p)\), đó là Y: Số lần thành công trong n lần thử.

Đối với biến \(Y ~ B(n,p)\), thì giá trị trung bình của nó là \(EY = n.p\) và phương sai là \(VarY = n.p.(1 - p)\).

Nếu \(Y_1, Y_2,..., Y_m\) là các biến ngẫu nhiên độc lập có cùng phân phối nhị thức \(B(n,p)\) thì \(Y = Y_1 + Y_2 + ... + Y_m\) là biến ngẫu nhiên có phân phối nhị thức \(B(mn,p)\).

3.4.6 Bảng tần số - Bảng tần suất

Bảng tần số - Bảng tần suất còn được gọi là bảng ngẫu nhiên. Khi lập bảng ngẫu nhiên cho 2 biến thì bảng đó được gọi là bảng ngẫu nhiên 2 chiều, nếu lập cho 3 biến thì gọi là bảng ngẫu nhiên 3 chiều và cứ thế tăng lên.

Đối với bảng tần số chúng ta quy ước biến phụ thuộc (dependent/outcome/response variable) được xắp xếp theo cột, biến độc lập (independent/explanatory/predictor variable) được xắp xếp theo hàng.

3.4.7 Phân phối Poisson

Phân phối Poisson với tham số \(\lambda (\lambda > 0\) có hàm phân phối xác suất như sau:

\(p_k = P(Y = k) = \frac{\lambda^k.e^{-\lambda}}{k!}, k = 0,1,2,...\)

Khi đó ta viết: \(Y\) ~ \(P(\lambda)\) để chỉ rằng: \(Y\) có phân phối xác suất Poisson với tham số \(\lambda\). Giá trị trung bình khi đó của nó là \(EY = \lambda\) và phương sai là \(VarY = \lambda\).

Nếu \(Y_1, Y_2,..., Y_m\) là các biến ngẫu nhiên độc lập có cùng phân phối Poisson \(P(\lambda)\) thì \(Y = Y_1 + Y_2 + ... + Y_m\) là biến ngẫu nhiên có phân phối Poisson \(P(m\lambda)\).

3.4.8 Rủi ro tương đối (Relative Risk)

Sự khác biệt giữa hai tỷ lệ \(\pi_1\) và \(\pi_2\) được đánh giá qua sai số tuyệt đối (\(\pi_1 - \pi_2\)) chỉ dựa vào khoảng cách giữa hai tỷ lệ này và bỏ qua độ lớn và ý nghĩa của chúng, nghĩa là khoảng cách như nhau thì đánh giá sai lệch như nhau.

Giả sử \(\pi_1\) và \(\pi_2\) là xác suất thành công của hai nhóm thì rủi ro tương đối giữa hai nhím này là tỷ lệ xác suất thành công của hai nhóm: \(\frac{\pi_1}{\pi_2}\)

Các tỷ lệ tổng thể \(\pi_1\) và \(\pi_2\) là chưa biết, ta nhận được các tỷ lệ mẫu tương ứng là \(f_1\) và \(f_2\). Khi đó, ta dùng rủi ro tương đối mẫu \(\frac{f_1}{f_2}\) để xấp xỉ cho rủi ro tương đối \(\frac{\pi_1}{\pi_2}\) trên tổng thể.

3.4.9 Tỷ lệ chênh (Odds Ratios)

Tỷ lệ cược (Odds) là tỷ lệ giữa khả năng thành công và khả năng không thành công trong cùng một điều kiện tác động của biến giải thích, ký hiệu \(odds_i\):

\(odds_i = \frac{\pi_i}{1 - \pi_i}\)

Tỷ lệ chênh (Odds Ratio) của hàng 1 và hàng 2, được ký hiệu là \(\theta\), là tỷ số giữa tỷ lệ cược của thành công ở hàng 1 và tỷ lệ cược của thành công ở hàng 2:

\(\theta = \frac{odds_1}{odds_2} = \frac{\pi_1(1 - \pi_2)}{\pi_2(1 - \pi_1)}\)

3.4.10 Chỉ số đánh giá mô hình AIC (Akaike Information Criterian)

Công thức được sử dụng cho AIC trong R là:

\(AIC = Deviance + 2.k\) (với k là tham số của mô hình)

AIC càng bé thì mô hình càng chính xác.

3.4.11 Chỉ số đánh giá mô hình Bier

Chỉ số Brier là trung bình của tổng bình phương các độ chênh lệch giữa giá trị của biến phụ thuộc Y và giá trị ước lượng của nó.

\(B = \frac{1}{n}\sum_{i = 1}^n(Y_i - \hat\pi_i)^2 = \frac{sum square resid}{n}\)

3.4.12 Chỉ số đánh giá mô hình Deviance

Deviance cũng là một tiêu chí rất phổ biến được sử dụng để đánh giá một mô hình hồi quy được ước lượng bởi phương pháp Hợp lý cực đại (ML). Một cách tổng quá, cũng giống như chỉ tiêu AIC, giá trị của Deviance càng nhỏ thì mô hình càng tốt.

3.4.13 Ma trận nhầm lẫn

Ma trận nhầm lẫn là yếu tố giúp so sánh giá trị thực tế và giá trị dự báo. Trong ma trận nhầm lẫn người ta xét các chỉ số: Độ chính xác toàn thể, độ nhạy, độ đặc hiệu,…

Trong ma trận nhầm lẫn, các chỉ số quan trọng được chú ý là độ nhạy và độ đặc hiệu. Các chỉ số này càng lớn càng tốt.

3.4.14 Độ nhạy và độ đặc hiệu

Độ nhạy (sensitivity) của một thí nghiệm: Là tỷ lệ (%) của số ca bị bệnh thực sự khi xét nghiệm và cho kết quả dương tính với tổng số ca bị bệnh. Công thức để tính độ nhạy:

Độ nhạy = số dương tính thật/(số đương tính thật + số âm tính giả)

Độ đặc hiệu (specificity) của một thí nghiệm: Là tỷ lệ (%) của số ca không bị bệnh và kết quả xét nghiệm không bị bệnh với tổng số người không bị bệnh. Công thức tính độ đặc hiệu:

Độ đặc hiệu = Số trường hợp âm tính thật/ (số trường hợp âm tính thật + số trường hợp dương tính giả)

4 CHƯƠNG 3. KẾT QUẢ MÔ HÌNH NGHIÊN CỨU

4.1 Thống kê mô tả từng biến

4.1.1 Biến định tính Taste

table(milknew$Taste)

## 
##   0   1 
## 480 579

table(milknew$Taste)/sum(table(milknew$Taste))*100

## 
##        0        1 
## 45.32578 54.67422

milknew |> ggplot( aes( x = Taste, y= after_stat(count))) + geom_bar(color='green', fill='blue') + geom_text(aes(label= scales :: percent(after_stat(count/sum(count)),accuracy=.01)), stat = 'count', color= 'black', vjust= -.5) + theme_classic() + xlab('Hương vị của sữa') + ylab('Số mẫu')

Nhìn vào các bảng, ta thấy:

• Có 480 mẫu sữa có hương vị không ngon, chiếm 45.33%.
  
• Có 579 mẫu sữa có hương vị ngon, chiếm 54.67%.

4.1.2 Biến định tính Odor

table(milknew$Odor)

## 
##   0   1 
## 601 458

table(milknew$Odor)/sum(table(milknew$Odor))*100

## 
##        0        1 
## 56.75165 43.24835

milknew |> ggplot( aes( x = Odor, y= after_stat(count))) + geom_bar(color='green', fill='blue') + geom_text(aes(label= scales :: percent(after_stat(count/sum(count)),accuracy=.01)), stat = 'count', color= 'black', vjust= -.5) + theme_classic() + xlab('Mùi của sữa') + ylab('Số mẫu')

Nhìn vào các bảng, ta thấy:

• Có 601 mẫu sữa có mùi không thơm, chiếm 56.75%.
  
• Có 458 mẫu sữa có mùi thơm, chiếm 43.25%.

4.1.3 Biến định tính Fat

table(milknew$Fat)

## 
##   0   1 
## 348 711

table(milknew$Fat)/sum(table(milknew$Fat))*100

## 
##        0        1 
## 32.86119 67.13881

milknew |> ggplot( aes( x = Fat, y= after_stat(count))) + geom_bar(color='green', fill='blue') + geom_text(aes(label= scales :: percent(after_stat(count/sum(count)),accuracy=.01)), stat = 'count', color= 'black', vjust= -.5) + theme_classic() + xlab('Độ béo của sữa') + ylab('Số mẫu')

Nhìn vào các bảng, ta thấy:

• Có 348 mẫu sữa có độ béo thấp, chiếm 32.86%.
  
• Có 711 mẫu sữa có độ béo cao, chiếm 67.14%.

4.1.4 Biến định tính Turbidity

table(milknew$Turbidity)

## 
##   0   1 
## 539 520

table(milknew$Turbidity)/sum(table(milknew$Turbidity))*100

## 
##        0        1 
## 50.89707 49.10293

milknew |> ggplot( aes( x = Turbidity, y= after_stat(count))) + geom_bar(color='green', fill='blue') + geom_text(aes(label= scales :: percent(after_stat(count/sum(count)),accuracy=.01)), stat = 'count', color= 'black', vjust= -.5) + theme_classic() + xlab('Độ đục của sữa') + ylab('Số mẫu')

Nhìn vào các bảng, ta thấy:

• Có 539 mẫu sữa có độ đục thấp, chiếm 50.90%.
  
• Có 520 mẫu sữa có độ đục cao, chiếm 49.1%.

4.1.5 Biến định tính Grade

table(milknew$Grade)

## 
##   0   1 
## 548 511

table(milknew$Grade)/sum(table(milknew$Grade))*100

## 
##        0        1 
## 51.74693 48.25307

milknew |> ggplot( aes( x = Grade, y= after_stat(count))) + geom_bar(color='green', fill='blue') + geom_text(aes(label= scales :: percent(after_stat(count/sum(count)),accuracy=.01)), stat = 'count', color= 'black', vjust= -.5) + theme_classic() + xlab('Chất lượng của sữa') + ylab('Số lượng')

Nhìn vào các bảng, ta thấy:

• Có 548 mẫu sữa có chất lượng kém, chiếm 51.75%.
  
• Có 511 mẫu sữa có chất lượng cao, chiếm 48.25%.

4.1.6 Biến định lượng pH

summary(milknew$pH)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3.00    6.50    6.70    6.63    6.80    9.50

sd(milknew$pH)

## [1] 1.399679

milknew |> ggplot(aes(pH))+geom_bar(aes(y=(..count..)),color="orange",fill="blue")

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Nhìn vào bảng, ta thấy:

• Độ pH của sữa dao động ở mức 3 đến 9.5, mức độ trung bình là 6.63.

• Có 25% số mẫu trong dữ liệu có độ pH thấp hơn 6.5.

• Có 50% số mẫu trong dữ liệu có độ pH thấp hơn 6.7.

• Có 75% số mẫu trong dữ liệu có độ pH thấp hơn 6.8.

• Độ pH có độ lệch chuẩn là 1.399679.

4.1.7 Biến định lượng Temprature

summary(milknew$Temprature)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   34.00   38.00   41.00   44.23   45.00   90.00

sd(milknew$Temprature)

## [1] 10.09836

milknew |> ggplot(aes(Temprature))+geom_bar(aes(y=(..count..)),color="orange",fill="blue")

Nhìn vào bảng, ta thấy:

• Nhiệt độ của sữa dao động ở mức 34 đến 90, mức độ trung bình là 44.23.

• Có 25% số mẫu trong dữ liệu có nhiệt độ thấp hơn 38.

• Có 50% số mẫu trong dữ liệu có nhiệt độ thấp hơn 41.

• Có 75% số mẫu trong dữ liệu có nhiệt độ thấp hơn 45.

• Độ pH có độ lệch chuẩn là 10,09836.

4.1.8 Biến định lượng Colour

summary(milknew$Colour)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   240.0   250.0   255.0   251.8   255.0   255.0

sd(milknew$Colour)

## [1] 4.307424

milknew |> ggplot(aes(Colour))+geom_bar(aes(y=(..count..)),color="orange",fill="blue")

Nhìn vào bảng, ta thấy:

• Màu của sữa dao động ở mức 240 đến 255, mức độ trung bình là 251.8.

• Có 25% số mẫu trong dữ liệu có màu thấp hơn 250.

• Có 50% số mẫu trong dữ liệu có màu thấp hơn 255.

• Có 75% số mẫu trong dữ liệu có màu thấp hơn 255.

• Độ pH có độ lệch chuẩn là 10.09836.

4.2 Thống kê mô tả và thống kê suy diễn hai biến

4.2.1 Mối quan hệ giữa biến Grade và biến Taste

1. Bảng tần số, tần suất

ts1 <- table(milknew$Grade, milknew$Taste)
ts1

##    
##       0   1
##   0 251 297
##   1 229 282

ts11 <- prop.table(ts1)
addmargins(ts11)

##      
##               0         1       Sum
##   0   0.2370161 0.2804533 0.5174693
##   1   0.2162417 0.2662890 0.4825307
##   Sum 0.4532578 0.5467422 1.0000000

Nhìn vào bảng, ta thấy:

• Có 282 mẫu sữa có hương vị ngon được đánh giá chất lượng cao, chiếm 26.63%.
• Có 297 mẫu sữa có hương vị ngon được đánh giá chất lượng thấp, chiếm 28.05%.
• Có 229 mẫu sữa có hương vị không ngon được đánh giá chất lượng cao, chiếm 21.62%.
• Có 251 mẫu sữa có hương vị không ngon được đánh giá chất lượng thấp, chiếm 23.70%.

2. Rủi ro tương đối và tỷ lệ chênh lệch

epitab(ts1, method = "riskratio")

## $tab
##    
##       0        p0   1        p1 riskratio    lower    upper   p.value
##   0 251 0.4580292 297 0.5419708  1.000000       NA       NA        NA
##   1 229 0.4481409 282 0.5518591  1.018245 0.912474 1.136277 0.7577342
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, rủi ro tương đối của chất lượng sữa theo hương vị là 1.0182. Tức là, tỉ lệ chất lượng sữa theo hương vị ngon cao hơn tỉ lệ chất lượng sữa theo hương vị không ngon. Tuy nhiên, giá trị p-value = 0.7577 cho thấy sự khác biệt này không đáng kể thống kê. Tức là, không có đủ bằng chứng thống kê để kết luận rằng tỷ lệ rủi ro giữa 2 nhóm khác nhau.

epitab(ts1, method = "oddsratio")

## $tab
##    
##       0        p0   1        p1 oddsratio     lower    upper   p.value
##   0 251 0.5229167 297 0.5129534  1.000000        NA       NA        NA
##   1 229 0.4770833 282 0.4870466  1.040713 0.8168966 1.325851 0.7577342
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, tỷ lệ chênh lệch của chất lượng sữa theo hương vị là 1,0407. Tức là, tỉ lệ chất lượng sữa theo hương vị ngon cao hơn tỉ lệ chất lượng sữa theo hương vị không ngon. Tuy nhiên, giá trị p-value = 0.7577 cho thấy sự khác biệt này không đáng kể thống kê. Tức là, không có đủ bằng chứng thống kê để kết luận rằng tỷ lệ chênh lệch giữa 2 nhóm khác nhau.

3. Biểu đồ

ggplot(milknew, aes(Grade, fill = Taste)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

4. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Taste độc lập với nhau

\(H_1\): Biến Grade và biến Taste không độc lập với nhau

gil <- table (milknew$Grade, milknew$Taste)
gil

##    
##       0   1
##   0 251 297
##   1 229 282

chisq.test(gil)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gil
## X-squared = 0.068246, df = 1, p-value = 0.7939

Vì p_value = 0.7939 > 0.05 nên ta chấp nhận giả thuyết \(H_0\). Nghĩa là biến Grade và biến Taste độc lập với nhau.

4.2.2 Mối quan hệ giữa biến Grade và biến Odor

1. Bảng tần số, tần suất

ts2 <- table(milknew$Grade, milknew$Odor)
ts2

##    
##       0   1
##   0 338 210
##   1 263 248

ts22 <- prop.table(ts2)
addmargins(ts22)

##      
##               0         1       Sum
##   0   0.3191690 0.1983003 0.5174693
##   1   0.2483475 0.2341832 0.4825307
##   Sum 0.5675165 0.4324835 1.0000000

Nhìn vào bảng, ta thấy:

• Có 248 mẫu sữa có mùi thơm được đánh giá chất lượng cao, chiếm 23.42%.
• Có 210 mẫu sữa có mùi thơm được đánh giá chất lượng thấp, chiếm 19.83%.
• Có 263 mẫu sữa có mùi không thơm được đánh giá chất lượng cao, chiếm 24.83%.
• Có 338 mẫu sữa có mùi không thơm được đánh giá chất lượng thấp, chiếm 31.92%.

2. Rủi ro tương đối và tỷ lệ chênh lệch

epitab(ts2, method = "riskratio")

## $tab
##    
##       0        p0   1        p1 riskratio    lower    upper      p.value
##   0 338 0.6167883 210 0.3832117  1.000000       NA       NA           NA
##   1 263 0.5146771 248 0.4853229  1.266462 1.102373 1.454975 0.0009941705
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, rủi ro tương đối của chất lượng sữa theo mùi thơm là 1,2665. Tức là, tỉ lệ chất lượng sữa theo mùi thơm cao hơn tỉ lệ chất lượng sữa theo mùi không thơm. Giá trị p-value = 0.001 rất nhỏ (gần bằng 0) cho thấy sự khác biệt này có ý nghĩa thống kê đáng kể.

epitab(ts2, method = "oddsratio")

## $tab
##    
##       0       p0   1        p1 oddsratio    lower    upper      p.value
##   0 338 0.562396 210 0.4585153  1.000000       NA       NA           NA
##   1 263 0.437604 248 0.5414847  1.517726 1.188589 1.938006 0.0009941705
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, tỷ lệ chênh lệch của chất lượng sữa theo hương vị là 1.5177. Tức là, tỉ lệ chất lượng sữa theo mùi thơm cao hơn tỉ lệ chất lượng sữa theo mùi không thơm. Giá trị p-value = 0.001 rất nhỏ (gần bằng 0) cho thấy sự khác biệt này có ý nghĩa thống kê đáng kể.

3. Biểu đồ

ggplot(milknew, aes(Grade, fill = Odor)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

4. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Odor độc lập với nhau

\(H_1\): Biến Grade và biến Odor không độc lập với nhau

gil <- table (milknew$Grade, milknew$Odor)
gil

##    
##       0   1
##   0 338 210
##   1 263 248

chisq.test(gil)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gil
## X-squared = 10.821, df = 1, p-value = 0.001004

Vì p_value = 0.001 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến Odor không độc lập với nhau.

4.2.3 Mối quan hệ giữa biến Grade và biến Fat

1. Bảng tần số, tần suất

ts3 <- table(milknew$Grade, milknew$Fat)
ts3

##    
##       0   1
##   0 193 355
##   1 155 356

ts33 <- prop.table(ts3)
addmargins(ts33)

##      
##               0         1       Sum
##   0   0.1822474 0.3352219 0.5174693
##   1   0.1463645 0.3361662 0.4825307
##   Sum 0.3286119 0.6713881 1.0000000

Nhìn vào bảng, ta thấy:

• Có 356 mẫu sữa có độ béo cao được đánh giá chất lượng cao, chiếm 33.62%.
• Có 355 mẫu sữa có độ béo cao được đánh giá chất lượng thấp, chiếm 33.52%.
• Có 155 mẫu sữa có độ béo thấp được đánh giá chất lượng cao, chiếm 14.64%.
• Có 193 mẫu sữa có độ béo thấp được đánh giá chất lượng thấp, chiếm 18.22%.

2. Rủi ro tương đối và tỷ lệ chênh lệch

epitab(ts3, method = "riskratio")

## $tab
##    
##       0        p0   1        p1 riskratio     lower    upper   p.value
##   0 193 0.3521898 355 0.6478102  1.000000        NA       NA        NA
##   1 155 0.3033268 356 0.6966732  1.075428 0.9886166 1.169862 0.1016743
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, rủi ro tương đối của chất lượng sữa theo độ béo là 1.0754. Tức là, tỉ lệ chất lượng sữa theo độ béo cao cao hơn tỉ lệ chất lượng sữa theo độ béo thấp. Tuy nhiên, giá trị p-value = 0.1017 cho thấy sự khác biệt này không đáng kể thống kê. Tức là không có đủ bằng chứng thống kể để kết luận rằng tỷ lệ rủi ro giữa 2 nhóm khác nhau.

epitab(ts3, method = "oddsratio")

## $tab
##    
##       0        p0   1        p1 oddsratio     lower   upper   p.value
##   0 193 0.5545977 355 0.4992968  1.000000        NA      NA        NA
##   1 155 0.4454023 356 0.5007032  1.248669 0.9652113 1.61537 0.1016743
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, tỷ lệ chênh lệch của chất lượng sữa theo độ béo là 1.2487. Tức là, tỉ lệ chất lượng sữa theo độ béo cao cao hơn tỉ lệ chất lượng sữa theo độ béo thấp. Tuy nhiên, giá trị p-value = 0.1017 cho thấy sự khác biệt này không đáng kể thống kê. Tức là không có đủ bằng chứng thống kể để kết luận rằng tỷ lệ rủi ro giữa 2 nhóm khác nhau.

3. Biểu đồ

ggplot(milknew, aes(Grade, fill = Fat)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

4. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Fat độc lập với nhau

\(H_1\): Biến Grade và biến Fat không độc lập với nhau

gil <- table (milknew$Grade, milknew$Fat)
gil

##    
##       0   1
##   0 193 355
##   1 155 356

chisq.test(gil)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gil
## X-squared = 2.6444, df = 1, p-value = 0.1039

Vì p_value = 2.2e-16 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến không độc lập với nhau.

4.2.4 Mối quan hệ giữa biến Grade và biến Turbidity

1. Bảng tần số, tần suất

ts4 <- table(milknew$Grade, milknew$Turbidity)
ts4

##    
##       0   1
##   0 222 326
##   1 317 194

ts44 <- prop.table(ts4)
addmargins(ts44)

##      
##               0         1       Sum
##   0   0.2096317 0.3078376 0.5174693
##   1   0.2993390 0.1831917 0.4825307
##   Sum 0.5089707 0.4910293 1.0000000

Nhìn vào bảng, ta thấy:

• Có 194 mẫu sữa có độ đục cao được đánh giá chất lượng cao, chiếm 18.32%.
• Có 326 mẫu sữa có độ đục cao được đánh giá chất lượng thấp, chiếm 30.78%.
• Có 317 mẫu sữa có độ đục thấp được đánh giá chất lượng cao, chiếm 29.93%.
• Có 222 mẫu sữa có độ đục thấp được đánh giá chất lượng thấp, chiếm 20.96%.

2. Rủi ro tương đối và tỷ lệ chênh lệch

epitab(ts4, method = "riskratio")

## $tab
##    
##       0        p0   1        p1 riskratio     lower     upper      p.value
##   0 222 0.4051095 326 0.5948905 1.0000000        NA        NA           NA
##   1 317 0.6203523 194 0.3796477 0.6381809 0.5600452 0.7272178 3.030392e-12
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, rủi ro tương đối của chất lượng sữa theo độ béo là 0.6382. Tức là, tỉ lệ chất lượng sữa theo độ đục cao thấp hơn tỉ lệ chất lượng sữa theo độ đục thấp. Giá trị p-value = 3.030392e-12 rất nhỏ (gần bằng 0) cho thấy sự khác biệt này là có ý nghĩa thống kê đáng kể.

epitab(ts4, method = "oddsratio")

## $tab
##    
##       0        p0   1        p1 oddsratio     lower     upper      p.value
##   0 222 0.4118738 326 0.6269231 1.0000000        NA        NA           NA
##   1 317 0.5881262 194 0.3730769 0.4167521 0.3255434 0.5335152 3.030392e-12
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Nhìn vào kết quả ta thấy, tỷ lệ chênh lệch của chất lượng sữa theo độ béo là 0.4168. Tức là, tỉ lệ chất lượng sữa theo độ đục cao thấp hơn tỉ lệ chất lượng sữa theo độ đục thấp. Giá trị p-value = 3.030392e-12 rất nhỏ (gần bằng 0) cho thấy sự khác biệt này là có ý nghĩa thống kê đáng kể.

3. Biểu đồ

ggplot(milknew, aes(Grade, fill = Turbidity)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

4. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Turbidity độc lập với nhau

\(H_1\): Biến Grade và biến Turbidity không độc lập với nhau

gil <- table (milknew$Grade, milknew$Turbidity)
gil

##    
##       0   1
##   0 222 326
##   1 317 194

chisq.test(gil)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gil
## X-squared = 48.161, df = 1, p-value = 3.926e-12

Vì p_value = 3.926e-12 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến Turbidity không độc lập với nhau.

4.2.5 Mối quan hệ giữa biến Grade và biến pH

Mã hóa biến định lượng pH thành định tính nhằm để thống kê mô tả biến

pH1 <- cut(milknew$pH, breaks = c(2,6.25,10), labels = c("không", "một"))

1. Bảng tần số, tần suất

ts5 <- table(milknew$Grade, pH1)
ts5

##    pH1
##     không một
##   0   169 379
##   1     0 511

ts55 <- prop.table(ts5)
addmargins(ts55)

##      pH1
##           không       một       Sum
##   0   0.1595845 0.3578848 0.5174693
##   1   0.0000000 0.4825307 0.4825307
##   Sum 0.1595845 0.8404155 1.0000000

Nhìn vào bảng, ta thấy:

• Có 511 mẫu sữa có độ pH phù hợp được đánh giá chất lượng cao, chiếm 48.25%.
• Có 379 mẫu sữa có độ pH phù hợp được đánh giá chất lượng thấp, chiếm 35.79%.
• Có 0 mẫu sữa có độ pH phù hợp được đánh giá chất lượng cao, chiếm 0%.
• Có 169 mẫu sữa có độ pH không phù hợp được đánh giá chất lượng thấp, chiếm 16%.

2. Biểu đồ

ggplot(milknew, aes(Grade, fill = pH1)) + geom_bar(position = 'dodge')

3. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến pH độc lập với nhau

\(H_1\): Biến Grade và biến pH không độc lập với nhau

gil <- table (milknew$Grade, milknew$pH)
gil

##    
##       3 4.5 4.7 5.5 5.6 6.4 6.5 6.6 6.7 6.8 7.4 8.1 8.5 8.6   9 9.5
##   0  70  37  20  23  19   0  65  40   8  58  37  24  22  40  61  24
##   1   0   0   0   0   0   1 124 119  74 191   2   0   0   0   0   0

chisq.test(gil)

## Warning in chisq.test(gil): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  gil
## X-squared = 553.63, df = 15, p-value < 2.2e-16

Vì p_value = 2.2e-16 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến pH không độc lập với nhau.

4.2.6 Mối quan hệ giữa biến Grade và biến Temprature

Mã hóa biến định lượng Temprature thành định tính nhằm để thống kê mô tả biến

Temprature1 <- cut(milknew$Temprature, breaks = c(33,45.5,91), labels = c("một", "không"))

1. Bảng tần số, tần suất

ts6 <- table(milknew$Grade, Temprature1)
ts6

##    Temprature1
##     một không
##   0 342   206
##   1 508     3

ts66 <- prop.table(ts6)
addmargins(ts66)

##      Temprature1
##               một       không         Sum
##   0   0.322946176 0.194523135 0.517469311
##   1   0.479697828 0.002832861 0.482530689
##   Sum 0.802644004 0.197355996 1.000000000

Nhìn vào bảng, ta thấy:

• Có 508 mẫu sữa có nhiệt độ phù hợp được đánh giá chất lượng cao, chiếm 47.97%.
• Có 342 mẫu sữa có nhiệt độ phù hợp được đánh giá chất lượng thấp, chiếm 32.29%.
• Có 3 mẫu sữa có nhiệt độ không phù hợp được đánh giá chất lượng cao, chiếm 0.003%.
• Có 206 mẫu sữa có nhiệt độ không phù hợp được đánh giá chất lượng thấp, chiếm 19.45%.

2. Biểu đồ

ggplot(milknew, aes(Grade, fill = Turbidity)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

3. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Turbidity độc lập với nhau

\(H_1\): Biến Grade và biến Turbidity không độc lập với nhau

gil <- table (milknew$Grade, milknew$Turbidity)
gil

##    
##       0   1
##   0 222 326
##   1 317 194

chisq.test(gil)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gil
## X-squared = 48.161, df = 1, p-value = 3.926e-12

Vì p_value = 3.926e-12 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến Turbidity không độc lập với nhau.

4.2.7 Mối quan hệ giữa biến Grade và biến Colour

Mã hóa biến định lượng Colour thành định tính nhằm để thống kê mô tả biến

Colour1 <- cut(milknew$Colour, breaks = c(239,243,256), labels = c("không", "một"))

1. Bảng tần số, tần suất

ts7 <- table(milknew$Grade, Colour1)
ts7

##    Colour1
##     không một
##   0    13 535
##   1    19 492

ts77 <- prop.table(ts7)
addmargins(ts77)

##      Colour1
##            không        một        Sum
##   0   0.01227573 0.50519358 0.51746931
##   1   0.01794145 0.46458924 0.48253069
##   Sum 0.03021719 0.96978281 1.00000000

Nhìn vào bảng, ta thấy:

• Có 492 mẫu sữa có màu phù hợp được đánh giá chất lượng cao, chiếm 46.46%.
• Có 535 mẫu sữa có màu phù hợp được đánh giá chất lượng thấp, chiếm 50.52%.
• Có 19 mẫu sữa có màu không phù hợp được đánh giá chất lượng cao, chiếm 1.79%.
• Có 13 mẫu sữa có màu không phù hợp được đánh giá chất lượng thấp, chiếm 1.23%.

2. Biểu đồ

ggplot(milknew, aes(Grade, fill = Colour)) + geom_bar(position = 'dodge')

## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

3. Kiểm định tính độc lập

Giả thuyết:

\(H_0\): Biến Grade và biến Turbidity độc lập với nhau

\(H_1\): Biến Grade và biến Turbidity không độc lập với nhau

gil <- table (milknew$Grade, milknew$Colour)
gil

##    
##     240 245 246 247 248 250 253 254 255
##   0  13  26  35   8  23 100   0   0 343
##   1  19  89   9  40   0  46  22   1 285

chisq.test(gil)

## Warning in chisq.test(gil): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  gil
## X-squared = 142.55, df = 8, p-value < 2.2e-16

Vì p_value = 2.2e-16 < 0.05 nên ta chấp nhận giả thuyết \(H_1\). Nghĩa là biến Grade và biến Colour không độc lập với nhau.

4.3 Khoảng ước lượng tỷ lệ

4.3.1 Ước lượng tỷ lệ chất lượng sữa

a <- milknew[milknew$Grade == "1",]
prop.test(length(a$Grade), length(milknew$Grade))

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(a$Grade) out of length(milknew$Grade), null probability 0.5
## X-squared = 1.2238, df = 1, p-value = 0.2686
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4520830 0.5131066
## sample estimates:
##         p 
## 0.4825307

Với khoảng tin cậy 95%, tỷ lệ chất lượng sữa cao so với tổng thể nằm trong khoảng 45.21% đến 51.31%. Hay nói cách khác, tỷ lệ chất lượng sữa thấp so với tổng thể nằm trong khoảng 48.69% đến 54.79%.

4.3.2 Ước lượng tỷ lệ mùi của sữa

b <- milknew[milknew$Odor == "1",]
prop.test(length(b$Odor), length(milknew$Odor))

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(b$Odor) out of length(milknew$Odor), null probability 0.5
## X-squared = 19.041, df = 1, p-value = 1.28e-05
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4024756 0.4629871
## sample estimates:
##         p 
## 0.4324835

Với khoảng tin cậy 95%, tỷ lệ mùi sữa thơm so với tổng thể nằm trong khoảng 40.25% đến 46.3%. Hay nói cách khác, tỷ lệ mùi sữa không thơm so với tổng thể nằm trong khoảng 53.7% đến 59.75%.

4.3.3 Ước lượng tỷ lệ độ béo của sữa

c <- milknew[milknew$Fat == "1",]
prop.test(length(c$Fat), length(milknew$Fat))

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(c$Fat) out of length(milknew$Fat), null probability 0.5
## X-squared = 123.74, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6420426 0.6994741
## sample estimates:
##         p 
## 0.6713881

Với khoảng tin cậy 95%, tỷ lệ độ béo sữa cao so với tổng thể nằm trong khoảng 64.20% đến 69.95%. Hay nói cách khác, tỷ lệ độ béo sữa thấp so với tổng thể nằm trong khoảng 30.05% đến 35.80%.

4.3.4 Ước lượng tỷ lệ độ đục của sữa

d <- milknew[milknew$Turbidity == "1",]
prop.test(length(d$Grade), length(milknew$Turbidity))

## 
##  1-sample proportions test with continuity correction
## 
## data:  length(d$Grade) out of length(milknew$Turbidity), null probability 0.5
## X-squared = 0.30595, df = 1, p-value = 0.5802
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4605370 0.5215875
## sample estimates:
##         p 
## 0.4910293

Với khoảng tin cậy 95%, tỷ lệ độ đục sữa cao so với tổng thể nằm trong khoảng 46.05% đến 52.16%. Hay nói cách khác, tỷ lệ độ đục sữa thấp so với tổng thể nằm trong khoảng 47.84% đến 53.95%.

4.4 Ước lượng mô hình hồi quy

4.4.1 Hồi quy với hàm probit

# Hồi quy với hàm probit
fit1 <- glm(Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + milknew$pH + milknew$Temprature + milknew$Colour , family = binomial(link = "probit"), data = milknew)
summary(fit1)

## 
## Call:
## glm(formula = Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + 
##     milknew$pH + milknew$Temprature + milknew$Colour, family = binomial(link = "probit"), 
##     data = milknew)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         7.049541   2.575967   2.737  0.00621 ** 
## milknew$Odor        0.642031   0.103548   6.200 5.63e-10 ***
## milknew$Fat         0.315081   0.097373   3.236  0.00121 ** 
## milknew$Turbidity  -0.818292   0.104527  -7.828 4.94e-15 ***
## milknew$pH          0.139542   0.031035   4.496 6.92e-06 ***
## milknew$Temprature -0.065175   0.006522  -9.993  < 2e-16 ***
## milknew$Colour     -0.021043   0.009951  -2.115  0.03446 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1466.8  on 1058  degrees of freedom
## Residual deviance: 1194.7  on 1052  degrees of freedom
## AIC: 1208.7
## 
## Number of Fisher Scoring iterations: 5

# Kiểm định sự phù hợp của mô hình bằng cách tính giá trị Prob(LR statistic)
lr_test <- anova(fit1, test = "Chisq")
p_value <- lr_test$Pr[2] 
p_value

## [1] 0.0007973808

Kiểm định sự phù hợp của mô hình:

• \(H_0\): Mô hình không phù hợp
• \(H_1\): Mô hình phù hợp

Với P-value = Prob(LR) < 0 chấp nhận giả thuyết \(H_1\) nên mô hình phù hợp với dữ liệu.

# Giá trị BrierScore
BrierScore(fit1)

## [1] 0.1878503

# Ma trận nhầm lẫn
predictions <- predict(fit1, newdata = milknew, type = "response")
predicted_classes <- ifelse(predictions > 0.5, "1", "0")
predictions1<-factor(predicted_classes, levels = c("0","1"))
actual<- factor(milknew$Grade, labels = c("0","1"))
confusionMatrix(table(predictions1, actual))

## Confusion Matrix and Statistics
## 
##             actual
## predictions1   0   1
##            0 374 102
##            1 174 409
##                                           
##                Accuracy : 0.7394          
##                  95% CI : (0.7118, 0.7656)
##     No Information Rate : 0.5175          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4806          
##                                           
##  Mcnemar's Test P-Value : 1.923e-05       
##                                           
##             Sensitivity : 0.6825          
##             Specificity : 0.8004          
##          Pos Pred Value : 0.7857          
##          Neg Pred Value : 0.7015          
##              Prevalence : 0.5175          
##          Detection Rate : 0.3532          
##    Detection Prevalence : 0.4495          
##       Balanced Accuracy : 0.7414          
##                                           
##        'Positive' Class : 0               
## 

4.4.2 Hồi quy với hàm logit

# Hồi quy với hàm logit
fit2 <- glm(Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + milknew$pH + milknew$Temprature + milknew$Colour , family = binomial(link = "logit"), data = milknew)
summary(fit2)

## 
## Call:
## glm(formula = Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + 
##     milknew$pH + milknew$Temprature + milknew$Colour, family = binomial(link = "logit"), 
##     data = milknew)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        13.22758    4.34496   3.044  0.00233 ** 
## milknew$Odor        1.20102    0.18119   6.628 3.39e-11 ***
## milknew$Fat         0.53201    0.16607   3.204  0.00136 ** 
## milknew$Turbidity  -1.41430    0.18200  -7.771 7.79e-15 ***
## milknew$pH          0.24414    0.05255   4.646 3.39e-06 ***
## milknew$Temprature -0.12650    0.01289  -9.815  < 2e-16 ***
## milknew$Colour     -0.03851    0.01672  -2.303  0.02126 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1466.8  on 1058  degrees of freedom
## Residual deviance: 1180.9  on 1052  degrees of freedom
## AIC: 1194.9
## 
## Number of Fisher Scoring iterations: 5

# Kiểm định sự phù hợp của mô hình bằng cách tính giá trị Prob(LR statistic)
lr_test <- anova(fit2, test = "Chisq")
p_value <- lr_test$Pr[2] 
p_value

## [1] 0.0007973808

Kiểm định sự phù hợp của mô hình:

• \(H_0\): Mô hình không phù hợp
• \(H_1\): Mô hình phù hợp

Với P-value = Prob(LR) < 0 chấp nhận giả thuyết \(H_1\) nên mô hình phù hợp với dữ liệu.

# Giá trị BrierScore
BrierScore(fit2)

## [1] 0.1865416

# Ma trận nhầm lẫn
predictions <- predict(fit2, newdata = milknew, type = "response")
predicted_classes <- ifelse(predictions > 0.5, "1", "0")
predictions1<-factor(predicted_classes, levels = c("0","1"))
actual<- factor(milknew$Grade, labels = c("0","1"))
confusionMatrix(table(predictions1, actual))

## Confusion Matrix and Statistics
## 
##             actual
## predictions1   0   1
##            0 357 108
##            1 191 403
##                                           
##                Accuracy : 0.7177          
##                  95% CI : (0.6895, 0.7446)
##     No Information Rate : 0.5175          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4377          
##                                           
##  Mcnemar's Test P-Value : 2.114e-06       
##                                           
##             Sensitivity : 0.6515          
##             Specificity : 0.7886          
##          Pos Pred Value : 0.7677          
##          Neg Pred Value : 0.6785          
##              Prevalence : 0.5175          
##          Detection Rate : 0.3371          
##    Detection Prevalence : 0.4391          
##       Balanced Accuracy : 0.7201          
##                                           
##        'Positive' Class : 0               
## 

4.4.3 Hồi quy với hàm cloglog

# Hồi quy với hàm probit
fit3 <- glm(Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + milknew$pH + milknew$Temprature + milknew$Colour , family = binomial(link = "cloglog"), data = milknew)
summary(fit3)

## 
## Call:
## glm(formula = Grade ~ milknew$Odor + milknew$Fat + milknew$Turbidity + 
##     milknew$pH + milknew$Temprature + milknew$Colour, family = binomial(link = "cloglog"), 
##     data = milknew)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         7.558708   2.925973   2.583 0.009786 ** 
## milknew$Odor        0.933924   0.123004   7.593 3.14e-14 ***
## milknew$Fat         0.476555   0.107909   4.416 1.00e-05 ***
## milknew$Turbidity  -1.147420   0.126655  -9.059  < 2e-16 ***
## milknew$pH          0.140271   0.037860   3.705 0.000211 ***
## milknew$Temprature -0.089932   0.009436  -9.531  < 2e-16 ***
## milknew$Colour     -0.021026   0.011094  -1.895 0.058059 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1466.8  on 1058  degrees of freedom
## Residual deviance: 1177.5  on 1052  degrees of freedom
## AIC: 1191.5
## 
## Number of Fisher Scoring iterations: 16

# Kiểm định sự phù hợp của mô hình bằng cách tính giá trị Prob(LR statistic)
lr_test <- anova(fit3, test = "Chisq")
p_value <- lr_test$Pr[2] 
p_value

## [1] 0.0007973808

Kiểm định sự phù hợp của mô hình:

• \(H_0\): Mô hình không phù hợp
• \(H_1\): Mô hình phù hợp

Với P-value = Prob(LR) < 0 chấp nhận giả thuyết \(H_1\) nên mô hình phù hợp với dữ liệu.

# Giá trị BrierScore
BrierScore(fit3)

## [1] 0.1878635

# Ma trận nhầm lẫn
predictions <- predict(fit3, newdata = milknew, type = "response")
predicted_classes <- ifelse(predictions > 0.5, "1", "0")
predictions1<-factor(predicted_classes, levels = c("0","1"))
actual<- factor(milknew$Grade, labels = c("0","1"))
confusionMatrix(table(predictions1, actual))

## Confusion Matrix and Statistics
## 
##             actual
## predictions1   0   1
##            0 362 138
##            1 186 373
##                                           
##                Accuracy : 0.6941          
##                  95% CI : (0.6653, 0.7217)
##     No Information Rate : 0.5175          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.3893          
##                                           
##  Mcnemar's Test P-Value : 0.009025        
##                                           
##             Sensitivity : 0.6606          
##             Specificity : 0.7299          
##          Pos Pred Value : 0.7240          
##          Neg Pred Value : 0.6673          
##              Prevalence : 0.5175          
##          Detection Rate : 0.3418          
##    Detection Prevalence : 0.4721          
##       Balanced Accuracy : 0.6953          
##                                           
##        'Positive' Class : 0               
## 

4.4.4 Lựa chọn mô hình phù hợp

Mô hình	AIC	Deviance	Brier Score	Độ chính xác	Độ nhạy	Độ đặc hiệu
Probit	1208.7	1194.9	0.1878503	0.7394	0.6825	0.8004
Logit	1194.9	1180.9	0.1865416	0.7177	0.6515	0.7886
Cloglog	1191.5	1177.5	0.1878635	0.6941	0.6606	0.7299
Lựa chọn	MH3	MH3	MH2	MH1	MH1	MH1

Từ kết quả thu được, mô hình Probit là mô hình thỏa mãn các tiêu chi đánh giá nhất trong 3 mô hình đề xuất. Do đó, mô hình Probit là mô hình tối ưu nhất.

4.4.5 Giải thích mô hình

Hệ số chặn là 7.049541, có nghĩa là khi tất cả các biến độc lập bằng không, xác suất chất lượng sữa là \(\phi\)(7.049541), trong đó \(\phi\) là hàm phân phối chuẩn tích lũy.

Hệ số Odor là 0.642031, có nghĩa là khi mùi của sữa tăng lên một đơn vị (từ 0 sang 1), xác suất chất lượng sữa tăng lên \(\phi\)(0.642031) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Hệ số Fat là 0.315081, có nghĩa là khi độ béo của sữa tăng lên một đơn vị (từ 0 sang 1), xác suất chất lượng sữa tăng lên \(\phi\)(0.315081) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Hệ số Turbidity của sữa là -0.818292, có nghĩa là khi độ đục của sữa tăng lên một đơn vị (từ 0 sang 1), xác suất chất lượng sữa giảm đi \(\phi\)(-0.818292) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Hệ số pH là 0.139542, có nghĩa là khi độ pH của sữa tăng lên một đơn vị, xác suất chất lượng sữa tăng lên \(\phi\)(0.139542) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Hệ số Temprature là -0.065175, có nghĩa là khi nhiệt độ của sữa tăng lên một đơn vị, xác suất chất lượng sữa giảm đi \(\phi\)(-0.065175) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Hệ số Colour của sữa là -0.021043, có nghĩa là khi màu của sữa tăng lên một đơn vị, xác suất chất lượng sữa giảm đi \(\phi\)(-0.021043) đơn vị chuẩn hoá, giữ nguyên các biến khác.

Các Pr(>|z|) của các hệ số đều nhỏ hơn 0.05, cho thấy rằng 6 biến đưa vào mô hình hồi quy để phân tích đều có ý nghĩa thống kê.

5 CHƯƠNG 4. KẾT LUẬN VÀ HẠN CHẾ

5.1 Kết luận

Với mục đích chính là Phân tích các yếu tố ảnh hưởng đến chất lượng sữa, tôi đã sử dụng bộ dữ liệu được thu thập từ năm 2022 để tiến hành đánh giá các yếu tố ảnh hưởng đến chất lượng sữa. Mô hình tối ưu giúp việc đánh giá hiệu quả là mô hình Probit.

Kết quả phân tích cho thấy 6 yếu tố trong bộ dữ liệu nhằm phân tích đánh giá chất lượng sữa đều có ý nghĩa. Qua đó, ta có thể phát triển tìm hiểu sâu hơn về các yếu tố này nhằm có thể đánh giá chất lượng sữa chính xác hơn nữa.

Qua bài phân tích, các hộ nông dân chăn nuôi đơn lẻ hoặc trang trại bò sữa có thể dựa theo đó mà tạo ra môi trường thích hợp nhất để đàn bò có thể cho ra các sản phẩm sữa chất lượng, dinh dưỡng nhất đến tay người tiêu dùng.

5.2 Hạn chế

Bộ dữ liệu được thu thập từ nước ngoài và được đăng tải trên internet nên tỷ lệ chính xác sẽ không tuyệt đối nên. Đó cũng chính là nguyên nhân dẫn đến một số bất cập trong quá trình nghiên cứu, làm xuất hiện các biến NA không xác định.

Bài phân tích được thực hiện trên phân mềm R & RStudio. Nhưng do đây cũng là phần mềm mới xuất hiện từ khoảng những năm 1993 nên việc có những nghiên cứu sâu sắc về nó tương đối ít. Cộng thêm sự non kém kinh nghiệm của sinh viên nên khi thực hiện nghiên cứu bằng phần mềm này còn nhiều bất cập và chỗ chưa chính xác.

6 TÀI LIỆU THAM KHẢO

[1] Rajendran, S. (2022) Milk quality prediction, Kaggle. Available at: https://www.kaggle.com/datasets/cpluzshrijayan/milkquality.

[2] Trần Kim Thanh, Trần Mạnh Tường, Vũ Anh Linh Duy: “Giáo trình Phân tích dữ liệu định tính”, trường Đại học Tài chính – Marketing, Thành phố Hồ Chí Minh.