Ngày 3: Hồi qui tuyến tính

Việc 1. Phân tích phương sai

1.1 Giả sử bạn có dữ liệu về cân nặng của 4 nhóm A, B, C và D như sau:

NhomA <- c(8, 9, 11, 4, 7, 8, 5)
NhomB <- c(7, 17, 10, 14, 12, 24, 11, 22)
NhomC <- c(28, 21, 26, 11, 24, 19)
NhomD <- c(26, 16, 13, 12, 9, 10, 11, 17, 15)

data <- data.frame(can_nang = c(NhomA, NhomB, NhomC, NhomD), nhom = rep(c("A", "B", "C", "D"), times = c(length(NhomA), length(NhomB), length(NhomC), length(NhomD))))

print(data)

##    can_nang nhom
## 1         8    A
## 2         9    A
## 3        11    A
## 4         4    A
## 5         7    A
## 6         8    A
## 7         5    A
## 8         7    B
## 9        17    B
## 10       10    B
## 11       14    B
## 12       12    B
## 13       24    B
## 14       11    B
## 15       22    B
## 16       28    C
## 17       21    C
## 18       26    C
## 19       11    C
## 20       24    C
## 21       19    C
## 22       26    D
## 23       16    D
## 24       13    D
## 25       12    D
## 26        9    D
## 27       10    D
## 28       11    D
## 29       17    D
## 30       15    D

1.2 Mô tả cân nặng giữa 4 nhóm

Chúng ta sẽ xem thống kê mô tả (trung bình, trung vị…) theo từng nhóm và vẽ biểu đồ boxplot.

tapply(data$can_nang, data$nhom, summary)

## $A
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.000   6.000   8.000   7.429   8.500  11.000 
## 
## $B
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00   10.75   13.00   14.62   18.25   24.00 
## 
## $C
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    11.0    19.5    22.5    21.5    25.5    28.0 
## 
## $D
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.00   11.00   13.00   14.33   16.00   26.00

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data %>% group_by(nhom) %>% summarize(so_luong = n(), trung_binh = mean(can_nang),trung_vi = median(can_nang), do_lech_chuan = sd(can_nang))

## # A tibble: 4 × 5
##   nhom  so_luong trung_binh trung_vi do_lech_chuan
##   <chr>    <int>      <dbl>    <dbl>         <dbl>
## 1 A            7       7.43      8            2.37
## 2 B            8      14.6      13            5.95
## 3 C            6      21.5      22.5          6.09
## 4 D            9      14.3      13            5.15

Vẽ Boxplot để xem phân phối

boxplot(can_nang ~ nhom, data = data,
        main = "Phân phối cân nặng theo 4 nhóm",
        xlab = "Nhóm",
        ylab = "Cân nặng")

### 1.3 Phân tích sự khác biệt về cân nặng giữa 4 nhóm. Diễn giải kết quả. Chúng ta dùng phép kiểm ANOVA (Analysis of Variance) để xem có ít nhất 2 nhóm có cân nặng trung bình khác nhau hay không.

model_anova <- aov(can_nang ~ nhom, data = data)

summary(model_anova)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## nhom         3  642.3  214.09   8.197 0.000528 ***
## Residuals   26  679.1   26.12                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diễn giải kết quả:

Nhìn vào hàng nhom trong bảng kết quả.
Tìm cột Pr(>F) (đây chính là p-value).
Nếu p < 0.05: Chúng ta kết luận: “Có sự khác biệt đáng kể về mặt thống kê về cân nặng trung bình giữa ít nhất hai trong bốn nhóm.” -Nếu p >= 0.05: Chúng ta kết luận: “Không có đủ bằng chứng để nói rằng có sự khác biệt về cân nặng trung bình giữa các nhóm.”

1.4 Thực hiện phân tích hậu định (post-hoc analysis) để xác định cụ thể nhóm có khác biệt về cân nặng. Diễn giải kết quả.

TukeyHSD(model_anova)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = can_nang ~ nhom, data = data)
## 
## $nhom
##           diff          lwr        upr     p adj
## B-A  7.1964286  -0.05969765 14.4525548 0.0525014
## C-A 14.0714286   6.27132726 21.8715299 0.0002134
## D-A  6.9047619  -0.16073856 13.9702624 0.0571911
## C-B  6.8750000  -0.69675602 14.4467560 0.0850381
## D-B -0.2916667  -7.10424368  6.5209103 0.9994049
## D-C -7.1666667 -14.55594392  0.2226106 0.0597131

Diễn giải kết quả:

-Kết quả sẽ hiển thị một bảng so sánh từng cặp (ví dụ: B-A, C-A, D-A, v.v.). - Nhìn vào cột p adj (p-value đã điều chỉnh). - Nếu p adj < 0.05 cho một cặp nào đó (ví dụ: C-A): Ta kết luận “Có sự khác biệt đáng kể về cân nặng trung bình giữa nhóm C và nhóm A.” - Nếu p adj >= 0.05 (ví dụ: B-A): Ta kết luận “Không có sự khác biệt đáng kể về cân nặng trung bình giữa nhóm B và nhóm A.” ## Việc 2. Phân tích tương quan ### 2.1 Đọc dữ liệu “Demo data.csv” vào R và gọi dữ liệu là “df”

df <- read.csv("C:\\Users\\Admin\\Desktop\\TH phan tich du lieu\\Demo.csv")

head(df)

##   X id age gender weight height pcfat
## 1 1  1  53      F     49    150  37.3
## 2 2  2  65      M     52    165  16.8
## 3 3  3  64      F     57    157  34.0
## 4 4  4  56      F     53    156  33.8
## 5 5  5  54      M     51    160  14.8
## 6 6  6  52      F     47    153  32.2

2.2 Mô tả đặc điểm cân nặng (weight) và chiều cao (height).

summary(df[, c("weight", "height")])

##      weight          height     
##  Min.   :34.00   Min.   :136.0  
##  1st Qu.:49.00   1st Qu.:151.0  
##  Median :54.00   Median :155.0  
##  Mean   :55.14   Mean   :156.7  
##  3rd Qu.:61.00   3rd Qu.:162.0  
##  Max.   :95.00   Max.   :185.0

2.3 Vẽ biểu đồ tán xạ đánh giá mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.

head(df)

##   X id age gender weight height pcfat
## 1 1  1  53      F     49    150  37.3
## 2 2  2  65      M     52    165  16.8
## 3 3  3  64      F     57    157  34.0
## 4 4  4  56      F     53    156  33.8
## 5 5  5  54      M     51    160  14.8
## 6 6  6  52      F     47    153  32.2

plot(x = df$height, y = df$weight,
     main = "Biểu đồ tán xạ giữa Chiều cao và Cân nặng",
     xlab = "Chiều cao (height)",
     ylab = "Cân nặng (weight)",
     pch = 19,
     col = "blue")

Nhận xét kết quả: Nhìn vào xu hướng của các dấu chấm.

Tương quan dương: Các dấu chấm đi lên từ trái sang phải (chiều cao tăng thì cân nặng có xu hướng tăng).
Tương quan âm: Các dấu chấm đi xuống từ trái sang phải.
Không tương quan: Các dấu chấm phân tán ngẫu nhiên, không có xu hướng rõ rệt. ### 2.4 Tiến hành phân tích tương quan định lượng mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.

cor.test(df$weight, df$height)

## 
##  Pearson's product-moment correlation
## 
## data:  df$weight and df$height
## t = 25.984, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5602911 0.6326135
## sample estimates:
##       cor 
## 0.5976667

Nhận xét kết quả: p-value:

Nếu p < 0.05, mối tương quan này có ý nghĩa thống kê. cor (hệ số tương quan r):
Giá trị này chạy từ -1 đến 1.
Gần +1: Tương quan dương mạnh.
Gần -1: Tương quan âm mạnh.
Gần 0: Tương quan yếu hoặc không có tương quan tuyến tính. ### 2.5 Tiến hành phân tích tương quan định lượng mối liên quan giữa chiều cao (height) và tỉ trọng mỡ (pcfat). Nhận xét kết quả.

cor.test(df$height, df$pcfat)

## 
##  Pearson's product-moment correlation
## 
## data:  df$height and df$pcfat
## t = -19.063, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5219407 -0.4353664
## sample estimates:
##        cor 
## -0.4798206

Nhận xét kết quả: (Tương tự như mục 2.4)

Xem p-value để biết có ý nghĩa thống kê hay không.
Xem cor (hệ số r) để biết độ mạnh và chiều hướng của mối tương quan. ## Việc 3. Hồi qui tuyến tính ### 3.1 Đọc dữ liệu “gapminder” vào R từ gói lệnh gapminder. Tạo bộ dữ liệu “vn” gồm dữ liệu của Việt Nam.

library(gapminder)

vn <- subset(gapminder, country == "Vietnam")

head(vn)

## # A tibble: 6 × 6
##   country continent  year lifeExp      pop gdpPercap
##   <fct>   <fct>     <int>   <dbl>    <int>     <dbl>
## 1 Vietnam Asia       1952    40.4 26246839      605.
## 2 Vietnam Asia       1957    42.9 28998543      676.
## 3 Vietnam Asia       1962    45.4 33796140      772.
## 4 Vietnam Asia       1967    47.8 39463910      637.
## 5 Vietnam Asia       1972    50.3 44655014      700.
## 6 Vietnam Asia       1977    55.8 50533506      714.

3.2 Thực hiện phân tích hồi qui tuyến tính để xác định xem mỗi năm, trong thời gian từ 1952-2007, tuổi thọ của người Việt Nam gia tăng như thế nào?

Chúng ta muốn xem year (năm) dự đoán lifeExp (tuổi thọ) như thế nào.

model_vn <- lm(lifeExp ~ year, data = vn)

summary(model_vn)

## 
## Call:
## lm(formula = lifeExp ~ year, data = vn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1884 -0.5840  0.1335  0.7396  1.7873 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.272e+03  4.349e+01  -29.25 5.10e-11 ***
## year         6.716e-01  2.197e-02   30.57 3.29e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared:  0.9894, Adjusted R-squared:  0.9884 
## F-statistic: 934.5 on 1 and 10 DF,  p-value: 3.289e-11

summary(model_vn)

## 
## Call:
## lm(formula = lifeExp ~ year, data = vn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1884 -0.5840  0.1335  0.7396  1.7873 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.272e+03  4.349e+01  -29.25 5.10e-11 ***
## year         6.716e-01  2.197e-02   30.57 3.29e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared:  0.9894, Adjusted R-squared:  0.9884 
## F-statistic: 934.5 on 1 and 10 DF,  p-value: 3.289e-11

3.3 Kiểm tra các giả định của mô hình hồi qui tuyến tính

Một mô hình hồi qui tốt cần tuân thủ một số giả định.

par(mfrow = c(2, 2))

plot(model_vn)

shapiro.test(residuals(model_vn))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_vn)
## W = 0.95663, p-value = 0.7349

par(mfrow = c(1, 1))

3.4 Viết phương trình đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007. Nhận xét kết quả.

Từ kết quả summary(model_vn) ở mục 3.2:Tìm giá trị (Intercept) trong cột Estimate. Đây là \(B_0\).Tìm giá trị year trong cột Estimate. Đây là \(B_1\).Phương trình (ví dụ):Giả sử (Intercept) là -1000 và year là 0.5.Phương trình sẽ là:\[Tuổi Thọ = -1000 + (0.5 \times Năm)\]

3.5 Sử dụng ChatGPT viết code để đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007 bằng dữ liệu sau:

year: 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, 2002, 2007 lifeExp: 40.4, 42.9, 45.4, 47.8, 50.3, 55.8, 58.8, 62.8, 67.7, 70.7, 73.0, 74.2

year <- c(1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, 2002, 2007)
lifeExp <- c(40.4, 42.9, 45.4, 47.8, 50.3, 55.8, 58.8, 62.8, 67.7, 70.7, 73.0, 74.2)
data_chat <- data.frame(year, lifeExp)

model_chat <- lm(lifeExp ~ year, data = data_chat)

summary(model_chat)

## 
## Call:
## lm(formula = lifeExp ~ year, data = data_chat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1494 -0.5944  0.1387  0.7324  1.8268 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.271e+03  4.377e+01  -29.04 5.47e-11 ***
## year         6.712e-01  2.211e-02   30.35 3.53e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.322 on 10 degrees of freedom
## Multiple R-squared:  0.9893, Adjusted R-squared:  0.9882 
## F-statistic: 921.4 on 1 and 10 DF,  p-value: 3.527e-11

Việc 4. Hồi qui tuyến tính đa biến

4.1 Giả sử bạn có dữ liệu về ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1) và thời gian lưu (X2) như sau:

Y: 12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6 X1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 X2: 7, 4, 4, 6, 4, 2, 1, 1, 1, 0 Nhập dữ liệu vào R và đặt tên tập dữ liệu là df.

Y <- c(12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6)
X1 <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
X2 <- c(7, 4, 4, 6, 4, 2, 1, 1, 1, 0)

df <- data.frame(Y, X1, X2)
print(df)

##       Y X1 X2
## 1  12.1  0  7
## 2  11.9  1  4
## 3  10.2  2  4
## 4   8.0  3  6
## 5   7.7  4  4
## 6   5.3  5  2
## 7   7.9  6  1
## 8   7.8  7  1
## 9   5.5  8  1
## 10  2.6  9  0

4.2 Đánh giá mối liên quan giữa ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1).

model_X1 <- lm(Y ~ X1, data = df)
summary(model_X1)

## 
## Call:
## lm(formula = Y ~ X1, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1606 -1.0735  0.1742  0.8621  2.0970 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.8545     0.8283  14.312 5.54e-07 ***
## X1           -0.8788     0.1552  -5.664 0.000474 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.409 on 8 degrees of freedom
## Multiple R-squared:  0.8004, Adjusted R-squared:  0.7755 
## F-statistic: 32.08 on 1 and 8 DF,  p-value: 0.0004737

Nhận xét: Nhìn vào hệ số Estimate và p-value của X1. ### 4.3 Đánh giá mối liên quan giữa ethylene oxide (Y) và thời gian lưu (X2).

model_X2 <- lm(Y ~ X2, data = df)
summary(model_X2)

## 
## Call:
## lm(formula = Y ~ X2, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.702 -1.533 -0.034  1.667  3.066 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   5.0980     1.1222   4.543  0.00189 **
## X2            0.9340     0.2999   3.114  0.01436 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.121 on 8 degrees of freedom
## Multiple R-squared:  0.548,  Adjusted R-squared:  0.4915 
## F-statistic: 9.698 on 1 and 8 DF,  p-value: 0.01436

4.4 Bạn muốn đánh giá mối liên quan độc lập giữa thời gian lưu (X2) và ethylene oxide (Y). Nhận xét kết quả.

model_multi <- lm(Y ~ X1 + X2, data = df)

summary(model_multi)

## 
## Call:
## lm(formula = Y ~ X1 + X2, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.46078 -0.33384  0.00026  0.81856  1.98476 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  14.7076     2.9785   4.938  0.00168 **
## X1           -1.2042     0.3614  -3.332  0.01255 * 
## X2           -0.4629     0.4642  -0.997  0.35187   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.41 on 7 degrees of freedom
## Multiple R-squared:  0.8252, Adjusted R-squared:  0.7753 
## F-statistic: 16.53 on 2 and 7 DF,  p-value: 0.002232

Nhận xét kết quả:

Tập trung vào hàng X2 trong bảng Coefficients của summary(model_multi).
Hệ số Estimate của X2: Cho biết: “Sau khi kiểm soát ảnh hưởng của X1, cứ mỗi 1 đơn vị X2 (thời gian lưu) tăng lên, thì Y (ethylene oxide) thay đổi [tăng/giảm] \(B_2\) đơn vị”.p-value của X2: Cho biết mối liên quan độc lập này có ý nghĩa thống kê hay không (nếu p < 0.05).
So sánh: Hãy so sánh hệ số Estimate và p-value của X2 trong mô hình này (model_multi) với mô hình đơn biến ở mục 4.3 (model_X2). Chúng có thể sẽ khác nhau!

Phân tích dữ liệu - day3

Lê Thủy Tiên - DAU

2025-10-23

Ngày 3: Hồi qui tuyến tính

Việc 1. Phân tích phương sai

1.1 Giả sử bạn có dữ liệu về cân nặng của 4 nhóm A, B, C và D như sau:

1.2 Mô tả cân nặng giữa 4 nhóm

Vẽ Boxplot để xem phân phối

Diễn giải kết quả:

1.4 Thực hiện phân tích hậu định (post-hoc analysis) để xác định cụ thể nhóm có khác biệt về cân nặng. Diễn giải kết quả.

Diễn giải kết quả:

2.2 Mô tả đặc điểm cân nặng (weight) và chiều cao (height).

2.3 Vẽ biểu đồ tán xạ đánh giá mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.

Nhận xét kết quả: Nhìn vào xu hướng của các dấu chấm.

Nhận xét kết quả: p-value:

Nhận xét kết quả: (Tương tự như mục 2.4)

3.2 Thực hiện phân tích hồi qui tuyến tính để xác định xem mỗi năm, trong thời gian từ 1952-2007, tuổi thọ của người Việt Nam gia tăng như thế nào?

3.3 Kiểm tra các giả định của mô hình hồi qui tuyến tính

3.4 Viết phương trình đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007. Nhận xét kết quả.

3.5 Sử dụng ChatGPT viết code để đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007 bằng dữ liệu sau:

Việc 4. Hồi qui tuyến tính đa biến

4.1 Giả sử bạn có dữ liệu về ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1) và thời gian lưu (X2) như sau:

4.2 Đánh giá mối liên quan giữa ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1).

4.4 Bạn muốn đánh giá mối liên quan độc lập giữa thời gian lưu (X2) và ethylene oxide (Y). Nhận xét kết quả.

Nhận xét kết quả: