Mô hình hồi quy tuyến tính dự đoán tuổi bào ngư

1 Giới thiệu bộ dữ liệu

Bộ dữ liệu là kết quả của việc dự đoán tuổi của bào ngư. Tuổi của bào ngư được xác định bằng cách cắt vỏ qua hình nón, nhuộm màu và đếm số vòng qua kính hiển vi.Ở đây ta sẽ sử dụng cách đo khác - sử dụng các phép đo vật lý để xác định tuổi của bào ngư.

2 Load thư viện

library(dplyr)
library(car)
library(psych)
library(ggplot2)
library(gridExtra)
library(MASS)
library(Hmisc)

3 Load dữ liệu

bộ dữ liệu ablone.csv được lưu vào biến df, bao gồm 4177 hàng tương ứng với số cá thể được nghiên cứu và 9 cột tương ứng với các biến sau : sex,length, diameter, heipght, whole_weight(tổng trọng lượng), shucked_weight, viscera_weight, shell_weight, rings.

df = read.csv('https://raw.githubusercontent.com/pnhuy/bioinfo/master/datasets/abalone/abalone.csv', header = TRUE)
attach(df)
head(df)

sex	length	diameter	height	whole_weight	shucked_weight	viscera_weight	shell_weight	rings
M	0.455	0.365	0.095	0.5140	0.2245	0.1010	0.150	15
M	0.350	0.265	0.090	0.2255	0.0995	0.0485	0.070	7
F	0.530	0.420	0.135	0.6770	0.2565	0.1415	0.210	9
M	0.440	0.365	0.125	0.5160	0.2155	0.1140	0.155	10
I	0.330	0.255	0.080	0.2050	0.0895	0.0395	0.055	7
I	0.425	0.300	0.095	0.3515	0.1410	0.0775	0.120	8

str(df )

## 'data.frame':    4177 obs. of  9 variables:
##  $ sex           : Factor w/ 3 levels "F","I","M": 3 3 1 3 2 2 1 1 3 1 ...
##  $ length        : num  0.455 0.35 0.53 0.44 0.33 0.425 0.53 0.545 0.475 0.55 ...
##  $ diameter      : num  0.365 0.265 0.42 0.365 0.255 0.3 0.415 0.425 0.37 0.44 ...
##  $ height        : num  0.095 0.09 0.135 0.125 0.08 0.095 0.15 0.125 0.125 0.15 ...
##  $ whole_weight  : num  0.514 0.226 0.677 0.516 0.205 ...
##  $ shucked_weight: num  0.2245 0.0995 0.2565 0.2155 0.0895 ...
##  $ viscera_weight: num  0.101 0.0485 0.1415 0.114 0.0395 ...
##  $ shell_weight  : num  0.15 0.07 0.21 0.155 0.055 0.12 0.33 0.26 0.165 0.32 ...
##  $ rings         : int  15 7 9 10 7 8 20 16 9 19 ...

Dưới đây là một vài chỉ số thống kê cơ bản:

df = na.omit(df)
summary(df)

##  sex          length         diameter          height        whole_weight   
##  F:1307   Min.   :0.075   Min.   :0.0550   Min.   :0.0000   Min.   :0.0020  
##  I:1342   1st Qu.:0.450   1st Qu.:0.3500   1st Qu.:0.1150   1st Qu.:0.4415  
##  M:1528   Median :0.545   Median :0.4250   Median :0.1400   Median :0.7995  
##           Mean   :0.524   Mean   :0.4079   Mean   :0.1395   Mean   :0.8287  
##           3rd Qu.:0.615   3rd Qu.:0.4800   3rd Qu.:0.1650   3rd Qu.:1.1530  
##           Max.   :0.815   Max.   :0.6500   Max.   :1.1300   Max.   :2.8255  
##  shucked_weight   viscera_weight    shell_weight        rings       
##  Min.   :0.0010   Min.   :0.0005   Min.   :0.0015   Min.   : 1.000  
##  1st Qu.:0.1860   1st Qu.:0.0935   1st Qu.:0.1300   1st Qu.: 8.000  
##  Median :0.3360   Median :0.1710   Median :0.2340   Median : 9.000  
##  Mean   :0.3594   Mean   :0.1806   Mean   :0.2388   Mean   : 9.934  
##  3rd Qu.:0.5020   3rd Qu.:0.2530   3rd Qu.:0.3290   3rd Qu.:11.000  
##  Max.   :1.4880   Max.   :0.7600   Max.   :1.0050   Max.   :29.000

Để xây dựng một mô hình hồi quy tuyến tính, trước hết ta phải xem mức độ tương quan giữa các biến thông qua biểu đồ và hệ số tương quan.

pairs.panels(df)

p1 = ggplot(data = df , aes(x = rings)) + 
  geom_histogram(aes(color = sex , fill = sex), alpha = 0.7 )

p2 = ggplot(data = df, aes(x = rings, fill = sex)) + 
  geom_density(position = 'stack', alpha = 0.7)

p3 = ggplot(data = df, aes(x = rings, fill = sex)) + geom_histogram(position = 'dodge')

p4 = ggplot(data = df , aes(x = rings)) + 
  geom_histogram(aes(color = sex , fill = sex), alpha = 0.7 ) +
  facet_wrap(~sex, ncol = 3)

grid.arrange(p1,p2,p3,p4)

p5 = ggplot(data = df, aes(y = rings, x = length, color = sex)) + geom_point()
p6 = ggplot(data = df, aes(y = rings, x = height, color = sex)) + geom_point()
p7 = ggplot(data = df, aes(y = rings, x = diameter, color = sex)) + geom_point()
p8 = ggplot(data = df, aes(y = rings, x = whole_weight, color = sex)) + geom_point()
p9 = ggplot(data = df, aes(y = rings, x = shucked_weight, color = sex)) + geom_point()
p10 = ggplot(data = df, aes(y = rings, x = viscera_weight, color = sex)) + geom_point()
p11 = ggplot(data = df, aes(y = rings, x = shell_weight, color = sex)) + geom_point()

grid.arrange(p5,p6,p7,p8,p9,p10,p11)

Từ hệ số tương quan và các đồ thị ở phía trên cho ta biết được mối tương quan giữa biến ring và các biến còn lại.

4 Kiểm định giả thuyết thống kê

anova = aov(rings ~ sex)
summary(anova)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## sex            2   8381    4191   499.3 <2e-16 ***
## Residuals   4174  35030       8                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rings ~ sex)
## 
## $sex
##           diff        lwr        upr     p adj
## I-F -3.2388418 -3.5027939 -2.9748896 0.0000000
## M-F -0.4238064 -0.6797094 -0.1679034 0.0003087
## M-I  2.8150354  2.5609373  3.0691334 0.0000000

Từ kết quả ở trên, ta thấy rằng giữa các nhóm giới tính sex có sự khác biệt về số vòng rings

5 Xây dựng mô hình hồi quy tuyến tính

Ở đây, ta xây dựng mô hình hồi quy tuyến tính với biến mục tiêu ring với các biến còn lại

m1 = lm(rings ~ ., data = df)
model_final = stepAIC(m1, direction = 'both',trace = FALSE)
summary(model_final)

## 
## Call:
## lm(formula = rings ~ sex + diameter + height + whole_weight + 
##     shucked_weight + viscera_weight + shell_weight, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4644  -1.3041  -0.3427   0.8633  13.9571 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.87038    0.27536  14.056  < 2e-16 ***
## sexI            -0.82644    0.10220  -8.087 7.96e-16 ***
## sexM             0.05755    0.08333   0.691     0.49    
## diameter        10.56951    0.98887  10.688  < 2e-16 ***
## height          10.74911    1.53525   7.002 2.94e-12 ***
## whole_weight     8.97751    0.72528  12.378  < 2e-16 ***
## shucked_weight -19.80258    0.81490 -24.301  < 2e-16 ***
## viscera_weight -10.61279    1.28782  -8.241 2.27e-16 ***
## shell_weight     8.75078    1.12405   7.785 8.73e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.194 on 4168 degrees of freedom
## Multiple R-squared:  0.5379, Adjusted R-squared:  0.537 
## F-statistic: 606.4 on 8 and 4168 DF,  p-value: < 2.2e-16

6 Diễn giải mô hình

Gọi các biến :

• sexI, sexM, diameter, height, whole_weight, shucked_weight, viscera_weight, shell_weight lần lượt là \(X_1\), \(X_2\), \(X_3\), \(X_4\), \(X_5\), \(X_6\) ,\(X_7\), \(X_8\)

• rings là \(Y\)

Mô hình hồi quy tuyến tính với biến rings có dạng như sau : \[Y = 3.87 - 0.82{X_1} + 0.05{X_2} + 10.56{X_3} + 10.75{X_4} + 8.97{X_5} - 19.80{X_6} - 10.61{X_7} + 8.75{X_8}\]

Dựa vào kết quả phân tích từ hàm summary() :

• Phần dư (residuals) : cho ta biết một vài thông số như Min, Max, Median

• Phần hệ số (Cofficients) : cho ta biết các ước số của mô hình, một vài chỉ số về sai số chuẩn ,giá trị P value,..

• Phần còn lại cho ta biết kết quả của các chỉ số như : R-square = 0.5379, Adjusted R squared = 0.537, F statistic = 606.4

Tiếp theo ta sẽ tính các giá trị Y dự đoán \(\hat{Y}\)

head(predict(model_final), 20)

##         1         2         3         4         5         6         7         8 
##  9.216457  7.848181 11.095069  9.640896  6.729240  7.826991 13.541942 11.467419 
##         9        10        11        12        13        14        15        16 
##  9.719014 13.133518 11.267105  9.459321 10.604171 10.437677 10.413905 11.098297 
##        17        18        19        20 
##  8.230932  9.079185  8.697501  8.639560

Cuối cùng ta thực hiện một vài phép kiểm với phần dư :

op = par(mfrow = c(2,2))
p12 = plot(model_final, c(1,2))
p13 = spreadLevelPlot(model_final)

## Warning in spreadLevelPlot.lm(model_final): 
## 2 negative fitted values removed

Từ các hình trên, có thể thấy :

• Trung bình phần dư không bằng 0
• Khồng tuân theo quy luật phân phối chuẩn
• Phương sai bất biến

7 Kết luận

Mô hình hồi quy tuyến tính xây dựng cho biến rings của bộ dữ liệu tương đối tốt khi đánh giá thông qua các chỉ số đã tính toán trong mô hình, tuy nhiên các kiểm định về các giả thiết của phần dư thì không đạt đối với một vài giả thiết chính.