PTDLDT_CT2_TUẦN 1

TASK 1: TÓM TẮT QUYỂN SÁCH “2019_Generalized Linear Models With Examples in R”

CHƯƠNG 1. GIỚI THIỆU MÔ HÌNH THỐNG KÊ

• Trích dẫn câu nói của nhà thống kê người Anh:

“All models are wrong, but some are useful” — George E. P. Box

“Tất cả các mô hình đều sai, nhưng một số mô hình có ích.”

Ý nghĩa:

• Không có mô hình nào phản ánh hoàn toàn thực tế một cách chính xác tuyệt đối. Bất kỳ mô hình thống kê nào cũng là một sự đơn giản hóa của thế giới thực.
• Do đó, tất cả các mô hình đều sai ở một mức độ nào đó, vì chúng luôn bỏ qua hoặc giả định những điều không hoàn toàn đúng.

Vậy tại sao lại “có ích”?

• Mặc dù sai, mô hình vẫn có thể giúp con người hiểu rõ xu hướng, mối quan hệ, hoặc đưa ra dự đoán hợp lý dựa trên dữ liệu.
• Một mô hình được xem là “có ích” khi nó hỗ trợ tốt cho mục đích cụ thể, chẳng hạn như:
  • Dự đoán chính xác
  • Hiểu rõ cơ chế sinh dữ liệu
  • Hỗ trợ ra quyết định

Bài học rút ra:

• Không nên tuyệt đối hóa mô hình.
• Nên dùng mô hình như một công cụ hỗ trợ tư duy và phân tích, thay vì là chân lý cuối cùng.
• Luôn cần kiểm tra giả định mô hình, đánh giá độ phù hợp, và diễn giải kết quả một cách cẩn trọng.

Ví dụ trong R:

model <- lm(mpg ~ wt, data = mtcars)
summary(model)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

1.1 Giới thiệu và tổng quan

Giới thiệu:

• Chương này giới thiệu khái niệm cơ bản về mô hình thống kê như một cách để mô tả cả các đặc điểm ngẫu nhiên và có hệ thống của dữ liệu. Nó nhấn mạnh tầm quan trọng của việc sử dụng các mô hình để phân tích dữ liệu.

• Là công cụ then chốt để hiểu và mô tả dữ liệu. Trong đó, mô hình tuyến tính tổng quát (GLM) là chủ đề chính của cuốn sách. Trước khi đi sâu vào GLM, chương này trình bày các khái niệm cơ bản giúp thiết lập ngôn ngữ, ký hiệu và phương pháp làm việc trong phân tích mô hình.

(“Data analysis: The need for models?” - Reese, 1986).

Mục đích của mô hình thống kê:

• Mục đích của một một một hình ảnh hưởng đến cách nó được phát triển (Mục 1.9). Các mục đích có thể bao gồm mô tả, dự đoán và hiểu mối quan hệ giữa các biến.

1.2 Quy ước mô tả dữ liệu

Giới thiệu cách mô tả dữ liệu dưới dạng toán học:

• Các biến phản hồi (\(y\)) và biến giải thích (\(x\)). Cách sắp xếp dữ liệu thành bảng quan sát và sử dụng ký hiệu chuẩn hóa để trình bày mô hình một cách thống nhất.

library(GLMsData)
data(lungcap)
summary(lungcap)

##       Age              FEV              Ht        Gender      Smoke        
##  Min.   : 3.000   Min.   :0.791   Min.   :46.00   F:318   Min.   :0.00000  
##  1st Qu.: 8.000   1st Qu.:1.981   1st Qu.:57.00   M:336   1st Qu.:0.00000  
##  Median :10.000   Median :2.547   Median :61.50           Median :0.00000  
##  Mean   : 9.931   Mean   :2.637   Mean   :61.14           Mean   :0.09939  
##  3rd Qu.:12.000   3rd Qu.:3.119   3rd Qu.:65.50           3rd Qu.:0.00000  
##  Max.   :19.000   Max.   :5.793   Max.   :74.00           Max.   :1.00000

str(lungcap)

## 'data.frame':    654 obs. of  5 variables:
##  $ Age   : int  3 4 4 4 4 4 4 5 5 5 ...
##  $ FEV   : num  1.072 0.839 1.102 1.389 1.577 ...
##  $ Ht    : num  46 48 48 48 49 49 50 46.5 49 49 ...
##  $ Gender: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Smoke : int  0 0 0 0 0 0 0 0 0 0 ...

1.3 Vẽ dữ liệu

Nhấn mạnh vai trò của việc trực quan hóa dữ liệu thông qua biểu đồ. Việc vẽ đồ thị giúp phát hiện xu hướng, bất thường và hiểu mối liên hệ giữa các biến trước khi mô hình hóa.

plot(FEV ~ Age, data = lungcap, main = "FEV vs Age", las = 1)

boxplot(FEV ~ Smoke + Gender, data = lungcap, las = 2)

interaction.plot(lungcap$Smoke, lungcap$Gender, lungcap$FEV,
                 xlab = "Smoke", ylab = "FEV", trace.label = "Gender", las = 1)

1.4 Mã hóa biến không phải số

Các biến phân loại (như giới tính, nhóm tuổi) cần được chuyển thành số để đưa vào mô hình. Mục này giải thích cách mã hóa chúng bằng biến giả (dummy variables) hoặc các phương pháp mã hóa khác.

lungcap$Gender <- factor(lungcap$Gender)
contrasts(lungcap$Gender)

##   M
## F 0
## M 1

lungcap$Smoke <- factor(lungcap$Smoke, levels = c(0, 1),
                        labels = c("Non-smoker", "Smoker"))
contrasts(lungcap$Smoke)

##            Smoker
## Non-smoker      0
## Smoker          1

1.5 Hai thành phần của mô hình thống kê

Mô hình thống kê gồm:\ - Thành phần hệ thống: mô tả xu hướng chính (predictor).\ - Thành phần ngẫu nhiên: mô tả sai số hoặc sự biến động không giải thích được.

model <- lm(FEV ~ Age + Ht + Gender + Smoke, data = lungcap)
summary(model)

## 
## Call:
## lm(formula = FEV ~ Age + Ht + Gender + Smoke, data = lungcap)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.37656 -0.25033  0.00894  0.25588  1.92047 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.456974   0.222839 -20.001  < 2e-16 ***
## Age          0.065509   0.009489   6.904 1.21e-11 ***
## Ht           0.104199   0.004758  21.901  < 2e-16 ***
## GenderM      0.157103   0.033207   4.731 2.74e-06 ***
## SmokeSmoker -0.087246   0.059254  -1.472    0.141    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4122 on 649 degrees of freedom
## Multiple R-squared:  0.7754, Adjusted R-squared:  0.774 
## F-statistic:   560 on 4 and 649 DF,  p-value: < 2.2e-16

1.6 Mô hình hồi quy

Giới thiệu lớp mô hình hồi quy – nền tảng cho tất cả các mô hình được trình bày trong sách. Mô hình hồi quy mô tả mối quan hệ tuyến tính giữa biến phản hồi và các biến giải thích.

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

Trình bày cách hiểu và giải thích ý nghĩa của các hệ số hồi quy. Cho biết mức độ ảnh hưởng của từng biến giải thích đến biến phản hồi.

1.8 So sánh mô hình vật lý và mô hình thống kê

So sánh hai loại mô hình: vật lý (dựa vào định luật tự nhiên) và thống kê (dựa vào dữ liệu). Mô hình thống kê có tính linh hoạt hơn nhưng cũng phụ thuộc vào dữ liệu nhiều hơn.

1.9 Mục đích của mô hình thống kê

Tùy mục tiêu mà mô hình có thể được dùng để: mô tả, dự đoán, hoặc kiểm định giả thuyết. Mục tiêu ảnh hưởng đến cách xây dựng, lựa chọn và đánh giá mô hình.

1.10 Đánh giá mô hình: Chính xác và đơn giản

Một mô hình tốt cần vừa chính xác (dự đoán hoặc mô tả tốt dữ liệu), vừa đơn giản (ít biến, dễ hiểu). Đây là hai tiêu chí quan trọng trong lựa chọn mô hình.

1.11 Hiểu giới hạn của mô hình

Mô hình thống kê không thể nắm bắt toàn bộ thực tế. Việc phân biệt giữa dữ liệu thực nghiệm (có kiểm soát) và quan sát (không kiểm soát) là rất quan trọng để hiểu được tính hợp lệ của suy luận.

1.12 Khả năng khái quát hóa mô hình

Mô hình tốt cần có khả năng áp dụng vào dữ liệu mới hoặc tổng thể. Mục này thảo luận về sự liên kết giữa cách lấy mẫu dữ liệu và khả năng khái quát hóa của mô hình.

1.13 Giới thiệu sử dụng R cho mô hình hóa

R là công cụ chính được sử dụng xuyên suốt sách. Các ví dụ và phân tích mô hình trong sách đều được trình bày bằng R với mã nguồn đầy đủ.

names(lungcap)

## [1] "Age"    "FEV"    "Ht"     "Gender" "Smoke"

factor(lungcap$Gender)

##   [1] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
##  [38] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
##  [75] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
## [112] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
## [149] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
## [186] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
## [223] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F
## [260] F F F F F F F F F F F F F F F F F F F F M M M M M M M M M M M M M M M M M
## [297] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [334] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [371] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [408] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [445] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [482] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [519] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M
## [556] M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M M F F F
## [593] F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F M
## [630] M M M M M M M M M M M M M M M M M M M M M M M M M
## Levels: F M

plot(lungcap$FEV ~ lungcap$Age)

CHƯƠNG 2. MÔ HÌNH HỒI QUY TUYẾN TÍNH

2.1 Giới thiệu

Mô hình hồi quy tuyến tính là một trong những công cụ cơ bản và phổ biến nhất trong thống kê. Mô hình này giúp xác định và định lượng mối quan hệ giữa một biến phản hồi (liên tục) và một hoặc nhiều biến giải thích (liên tục hoặc phân loại).

Mô hình hồi quy tuyến tính là nền tảng của nhiều mô hình thống kê. Chúng mô tả mối quan hệ tuyến tính giữa một biến phản hồi liên tục và một hoặc nhiều biến giải thích.

2.2 Định nghĩa mô hình hồi quy tuyến tính

Mô hình hồi quy tuyến tính tổng quát biểu diễn mối liên hệ tuyến tính giữa biến phản hồi và các biến giải thích. Biến sai số ε đại diện cho nhiễu loạn ngẫu nhiên và giả định phân phối chuẩn có trung bình 0 và phương sai ^2.

Mô hình hồi quy tuyến tính có dạng như sau:

\[ y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} + \varepsilon_i,\quad \varepsilon_i \sim \mathcal{N}(0, \sigma^2) \]

2.3 Hồi quy tuyến tính đơn

Đây là trường hợp đặc biệt của mô hình hồi quy tuyến tính với chỉ một biến giải thích. Mục tiêu là ước lượng mối quan hệ tuyến tính giữa biến đầu vào x và biến đầu ra y.

library(GLMsData)
data(gestation)
x <- gestation$Age
y <- gestation$Weight
wts <- gestation$Births

2.3.1 Ước lượng bình phương nhỏ nhất (OLS)

Phương pháp này có dạng như sau:

\[ \min_{B_0, B_1} \sum_{i=1}^n (y_i - B_0 - B_1 x_i)^2 \]

beta0.A <- -0.9; beta1.A <- 0.1
mu.A <- beta0.A + beta1.A * x
SA <- sum(wts * (y - mu.A)^2); SA

## [1] 186.1106

2.3.2 Ước lượng hệ số:

Hệ số góc (slope) được tính theo công thức:

\[ \beta_1 = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n} (x_i - \bar{x})^2} \]

Hệ số chặn (intercept) là:

\[ \beta_0 = \bar{y} - \beta_1 \bar{x} \]

xbar <- weighted.mean(x, w=wts)
SSx <- sum(wts * (x - xbar)^2)
ybar <- weighted.mean(y, w=wts)
SSxy <- sum(wts * (x - xbar) * y)
beta1 <- SSxy / SSx
beta0 <- ybar - beta1 * xbar
mu <- beta0 + beta1 * x
RSS <- sum(wts * (y - mu)^2)
c(beta0=beta0, beta1=beta1, RSS=RSS)

##      beta0      beta1        RSS 
## -2.6783891  0.1537594 11.4198322

2.3.3 Ước lượng phương sai và phần dư:

Sau khi ước lượng được mô hình hồi quy tuyến tính đơn:

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \]

trong đó \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\), ta cần ước lượng phương sai của sai số ngẫu nhiên để đánh giá mức độ phân tán của dữ liệu so với mô hình.

Phương sai sai số được ước lượng theo công thức:

\[ \hat{\sigma}^2 = \frac{1}{n - 2} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]

Trong đó:

• \(\hat{\sigma}^2\): là ước lượng phương sai sai số (residual variance),
• \(y_i\): là giá trị quan sát,
• \(\hat{y}_i\): là giá trị dự đoán từ mô hình,
• \(n\): là số lượng quan sát,
• \(n - 2\): là bậc tự do, do có hai hệ số \(\beta_0\) và \(\beta_1\) được ước lượng.

df <- length(y) - 2
s2 <- RSS / df
c(df = df, s = sqrt(s2), s2 = s2)

##         df          s         s2 
## 19.0000000  0.7752701  0.6010438

Phương sai sai số càng nhỏ thì mô hình càng phù hợp với dữ liệu thực tế.

Ước lượng phương sai phần dư:

\[ s^2 = \frac{RSS}{n - p'} \] #### 2.3.4 Sai số chuẩn:

Sai số chuẩn của hệ số góc được tính theo công thức:

\[ SE(\hat{\beta}_1) = \sqrt{ \frac{\hat{\sigma}^2}{\sum_{i=1}^{n} (x_i - \bar{x})^2} } \]

Trong đó:

• \(\hat{\sigma}^2\): là ước lượng phương sai sai số,
• \(x_i\): là giá trị của biến độc lập,
• \(\bar{x}\): là trung bình của \(x\),
• \(SE(\hat{\beta}_1)\): là sai số chuẩn của hệ số hồi quy \(\hat{\beta}_1\).

var.b0 <- s2 * (1/sum(wts) + xbar^2 / SSx)
var.b1 <- s2 / SSx
sqrt(c(beta0 = var.b0, beta1 = var.b1))

##       beta0       beta1 
## 0.371172341 0.009493212

Sai số chuẩn càng nhỏ thì ước lượng \(\hat{\beta}_1\) càng chính xác.

2.4 Hồi quy tuyến tính bội

Khi có nhiều biến giải thích, mô hình trở thành hồi quy tuyến tính bội. Việc thêm biến cho phép mô hình giải thích tốt hơn sự biến thiên trong dữ liệu, nhưng cũng làm tăng nguy cơ đa cộng tuyến hoặc overfitting.

Tương tự hồi quy đơn nhưng có nhiều biến x. Phân tích bằng OLS vẫn áp dụng.

2.5 Dạng ma trận của mô hình hồi quy

Viết mô hình dưới dạng ma trận giúp biểu diễn ngắn gọn và thuận tiện hơn khi xử lý bằng máy tính. Việc ước lượng hệ số sử dụng công thức đại số tuyến tính là nền tảng của nhiều thuật toán hồi quy hiện đại.

Mô hình hồi quy tuyến tính nhiều biến có thể được viết gọn dưới dạng ma trận như sau:

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

Trong đó:

• \(\mathbf{y}\): là vector \(n \times 1\) chứa các giá trị biến phụ thuộc,
• \(\mathbf{X}\): là ma trận thiết kế \(n \times p\) (gồm 1 cột toàn 1 cho hệ số chặn và các cột biến độc lập),
• \(\boldsymbol{\beta}\): là vector hệ số hồi quy \(p \times 1\),
• \(\boldsymbol{\varepsilon}\): là vector sai số ngẫu nhiên \(n \times 1\).

Hệ số hồi quy được ước lượng theo công thức bình phương tối thiểu:

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \]

Điều này giả định rằng \(\mathbf{X}^T \mathbf{X}\) khả nghịch (nghĩa là không có đa cộng tuyến hoàn toàn).

library(GLMsData)
data(lungcap)
X <- model.matrix(~ Age + Ht + Gender + Smoke, data = lungcap)
y <- log(lungcap$FEV)
beta <- solve(t(X) %*% X) %*% t(X) %*% y
fitted <- X %*% beta

2.6 Ước lượng bằng R

Ví dụ Mô hình hồi quy trong R, có thể dùng lệnh sau:

Mô hình hồi quy đơn:

\[ model_simple <- lm(mpg ~ wt, data = mtcars) \] Mô hình hồi quy bội

\[ model_multi <- lm(mpg ~ wt + hp, data = mtcars) \]

LC.m1 <- lm(log(FEV) ~ Age + Ht + Gender + Smoke, data = lungcap)
summary(LC.m1)

## 
## Call:
## lm(formula = log(FEV) ~ Age + Ht + Gender + Smoke, data = lungcap)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63278 -0.08657  0.01146  0.09540  0.40701 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.943998   0.078639 -24.721  < 2e-16 ***
## Age          0.023387   0.003348   6.984  7.1e-12 ***
## Ht           0.042796   0.001679  25.489  < 2e-16 ***
## GenderM      0.029319   0.011719   2.502   0.0126 *  
## Smoke       -0.046068   0.020910  -2.203   0.0279 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1455 on 649 degrees of freedom
## Multiple R-squared:  0.8106, Adjusted R-squared:  0.8095 
## F-statistic: 694.6 on 4 and 649 DF,  p-value: < 2.2e-16

coef(LC.m1)

## (Intercept)         Age          Ht     GenderM       Smoke 
## -1.94399818  0.02338721  0.04279579  0.02931936 -0.04606754

2.7 Diễn giải hệ số hồi quy

Mỗi hệ số hồi quy biểu diễn sự thay đổi trung bình của biến phản hồi y khi biến giải thích tăng một đơn vị, trong khi các biến khác giữ nguyên. Việc diễn giải đúng giúp đưa ra kết luận có ý nghĩa thực tiễn.

Hệ số \(\beta_j\): đại diện cho ảnh hưởng biên của biến \(x_j\) đến \(y\).

Mỗi khi \(x_j\) thay đổi một đơn vị (giữ các biến khác không đổi), thì \(y\) thay đổi trung bình \(\beta_j\) đơn vị.

Hay nói cách khác:

Mỗi 1 đơn vị thay đổi trong \(x_j\) dẫn đến thay đổi trung bình \(\beta_j\) đơn vị trong \(y\), giả sử các biến khác được giữ nguyên.

Ví dụ hệ số của chiều cao là 0.0428 → log(FEV) tăng 0.0428 nếu chiều cao tăng 1 đơn vị (giữ các biến khác không đổi).

2.8 Suy luận thống kê

Sau khi ước lượng, cần đánh giá xem các hệ số hồi quy có ý nghĩa thống kê không. Ta sử dụng kiểm định giả thuyết và khoảng tin cậy để đưa ra kết luận về ý nghĩa và độ chính xác của các ước lượng.

Kiểm định:

\[ H_0 : \beta_j = 0 \quad \text{vs} \quad H_1 : \beta_j \ne 0 \]

Khoảng tin cậy:

\[ \hat{\beta}_j \pm t_{n - p, \alpha/2} \cdot SE(\hat{\beta}_j) \]

summary(LC.m1)

## 
## Call:
## lm(formula = log(FEV) ~ Age + Ht + Gender + Smoke, data = lungcap)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63278 -0.08657  0.01146  0.09540  0.40701 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.943998   0.078639 -24.721  < 2e-16 ***
## Age          0.023387   0.003348   6.984  7.1e-12 ***
## Ht           0.042796   0.001679  25.489  < 2e-16 ***
## GenderM      0.029319   0.011719   2.502   0.0126 *  
## Smoke       -0.046068   0.020910  -2.203   0.0279 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1455 on 649 degrees of freedom
## Multiple R-squared:  0.8106, Adjusted R-squared:  0.8095 
## F-statistic: 694.6 on 4 and 649 DF,  p-value: < 2.2e-16

2.9 Phân tích phương sai (ANOVA)

Phân tích phương sai giúp đánh giá mức độ phù hợp của mô hình bằng cách chia tổng phương sai thành phần giải thích được (SSR) và không giải thích được (SSE). Đây là bước quan trọng để kiểm định tổng thể mô hình.

Tách tổng phương sai:

• SSR: do mô hình

• SSE: phần dư

2.10 So sánh mô hình lồng nhau

Để so sánh hai mô hình trong R \[ - model1 <- lm(mpg ~ wt, data = mtcars) \\ - model2 <- lm(mpg ~ wt + hp, data = mtcars) \\ - anova(model1, model2) \]

2.11 So sánh mô hình không lồng nhau (AIC/BIC)

• AIC(model1, model2)
• BIC(model1, model2)

2.12 Chọn biến trong mô hình

• step(model2)

2.13 Nghiên cứu tình huống

Phần này trình bày một ví dụ thực tế áp dụng mô hình hồi quy tuyến tính vào phân tích dữ liệu. Qua đó, người học có thể rèn luyện kỹ năng mô hình hóa, kiểm định, và diễn giải kết quả mô hình

Phân tích tập dữ liệu thực tế, ước lượng, diễn giải và đánh giá mô hình.

2.14 Dự đoán với predict()

\[ predict(model_multi, newdata = data.frame(wt = 3, hp = 120)) \]

2.15 Tóm tắt

Chương này đã trình bày đầy đủ từ định nghĩa đến thực hành mô hình hồi quy tuyến tính. Người đọc nắm được công cụ lý thuyết và kỹ năng thực hành để áp dụng vào phân tích dữ liệu định lượng.

Mô hình hồi quy tuyến tính là nền tảng của GLM

OLS ước lượng hệ số tốt nhất

R cung cấp công cụ hiệu quả để ước lượng, kiểm định và dự đoán

Chương 3.MÔ HÌNH TUYẾN TÍNH TỔNG QUÁT (GLM)

3.1 Giới thiệu và tổng quan

GLM là sự mở rộng của hồi quy tuyến tính để xử lý các dạng dữ liệu như đếm, nhị phân, dương liên tục… GLM gồm hai phần: phân phối (EDM) và hàm liên kết.

3.2 Họ phân phối mũ (Exponential Family)

\[ P(y; \theta, \phi) = a(y, \phi) \exp\left\{ \frac{y \theta - \kappa(\theta)}{\phi} \right\} \]

Ví dụ: - Normal: \(\theta = \mu\) - Poisson: \(\theta = \log(\mu)\) - Binomial: \(\theta = \log(\mu / (1 - \mu))\) - Gamma: \(\theta = -1/\mu\)

3.3 Kỳ vọng và phương sai

\[ E[y] = \kappa'(\theta) = \mu, \quad \text{Var}[y] = \phi \kappa''(\theta) = \phi V(\mu) \]

3.4 Mô hình GLM

\[ g(\mu_i) = o_i + \beta_0 + \sum_{j=1}^{p} \beta_j x_{ji} \]

# Kiểm tra trước khi dùng
str(lungcap)

## 'data.frame':    654 obs. of  5 variables:
##  $ Age   : int  3 4 4 4 4 4 4 5 5 5 ...
##  $ FEV   : num  1.072 0.839 1.102 1.389 1.577 ...
##  $ Ht    : num  46 48 48 48 49 49 50 46.5 49 49 ...
##  $ Gender: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Smoke : int  0 0 0 0 0 0 0 0 0 0 ...

# Mô hình Gaussian
glm(FEV ~ Age + Ht + Gender + Smoke, family = gaussian(link = "identity"), data = lungcap)

## 
## Call:  glm(formula = FEV ~ Age + Ht + Gender + Smoke, family = gaussian(link = "identity"), 
##     data = lungcap)
## 
## Coefficients:
## (Intercept)          Age           Ht      GenderM        Smoke  
##    -4.45697      0.06551      0.10420      0.15710     -0.08725  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  649 Residual
## Null Deviance:       490.9 
## Residual Deviance: 110.3     AIC: 703.8

# Mô hình Poisson - tạo dữ liệu phù hợp
set.seed(1)
lungcap$y_pois <- rpois(nrow(lungcap), lambda = 2)
glm(y_pois ~ Age, family = poisson(link = "log"), data = lungcap)

## 
## Call:  glm(formula = y_pois ~ Age, family = poisson(link = "log"), data = lungcap)
## 
## Coefficients:
## (Intercept)          Age  
##    0.748181    -0.004938  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  652 Residual
## Null Deviance:       718.6 
## Residual Deviance: 718.3     AIC: 2224

# Mô hình logistic - tạo biến nhị phân
lungcap$y_bin <- ifelse(lungcap$Smoke == "yes", 1, 0)
glm(y_bin ~ Age, family = binomial(link = "logit"), data = lungcap)

## Warning: glm.fit: algorithm did not converge

## 
## Call:  glm(formula = y_bin ~ Age, family = binomial(link = "logit"), 
##     data = lungcap)
## 
## Coefficients:
## (Intercept)          Age  
##  -2.657e+01    1.613e-14  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  652 Residual
## Null Deviance:       0 
## Residual Deviance: 3.794e-09     AIC: 4

3.5 Ước lượng hệ số

Ước lượng hợp lý tối đa (MLE) qua thuật toán IRLS (iteratively reweighted least squares).

x1 <- lungcap$Age
x2 <- lungcap$Ht
y <- lungcap$y_pois
glm(y ~ x1 + x2, family = poisson)

## 
## Call:  glm(formula = y ~ x1 + x2, family = poisson)
## 
## Coefficients:
## (Intercept)           x1           x2  
##     0.18541     -0.02394      0.01228  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  651 Residual
## Null Deviance:       718.6 
## Residual Deviance: 715.9     AIC: 2224

3.6 Deviance

\[ D(y, \hat{\mu}) = 2 \sum w_i \left[ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right] \]

3.7 Quy trình ước lượng GLM

1. Chọn phân phối (EDM)
2. Chọn hàm liên kết
3. Xây dựng log-likelihood
4. Giải bằng IRLS
5. Đánh giá sai số chuẩn, p-value

# Gán biến trước khi dùng
x <- lungcap$Age
y <- lungcap$y_pois

fit <- glm(y ~ x, family = poisson)
summary(fit)

## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.748181   0.096641   7.742  9.8e-15 ***
## x           -0.004938   0.009367  -0.527    0.598    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 718.61  on 653  degrees of freedom
## Residual deviance: 718.33  on 652  degrees of freedom
## AIC: 2224.3
## 
## Number of Fisher Scoring iterations: 5

3.8 So sánh với hồi quy tuyến tính

• Hồi quy tuyến tính là trường hợp đặc biệt của GLM (Normal + identity link).
• GLM mở rộng để xử lý biến đếm, nhị phân, dương liên tục…
• Các khái niệm như phần dư, leverage vẫn được giữ lại thông qua IRLS.

3.9 Tóm tắt chương

GLM là khung mô hình hóa rất mạnh, áp dụng cho nhiều loại dữ liệu, dựa trên phân phối thuộc họ mũ và hàm liên kết.

Chương 4.ƯỚC LƯỢNG VÀ SUY LUẬN TRONG GLM

4.1 Giới thiệu và tổng quan

Sử dụng hợp lý tối đa (MLE) để ước lượng hệ số và tham số phân tán \(\phi\). Kiểm định Wald, LRT và Score được giới thiệu để kiểm định giả thuyết.

4.2 Ước lượng \(\phi\)

\tilde{\phi} = \frac{D(y, \hat{\mu})}{n - p'},\quad
\bar{\phi} = \frac{1}{n - p'} \sum \frac{w_i (y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)}

#3# 4.3 Ước lượng GLM bằng R

x1 <- lungcap$Age
x2 <- lungcap$Ht
y <- lungcap$y_pois
glm(y ~ x1 + x2, family = poisson)

## 
## Call:  glm(formula = y ~ x1 + x2, family = poisson)
## 
## Coefficients:
## (Intercept)           x1           x2  
##     0.18541     -0.02394      0.01228  
## 
## Degrees of Freedom: 653 Total (i.e. Null);  651 Residual
## Null Deviance:       718.6 
## Residual Deviance: 715.9     AIC: 2224

glm.out <- glm(y ~ x1 + x2, family = poisson, data = lungcap)
summary(glm.out)

## 
## Call:
## glm(formula = y ~ x1 + x2, family = poisson, data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.185413   0.376768   0.492    0.623
## x1          -0.023939   0.015577  -1.537    0.124
## x2           0.012275   0.007938   1.546    0.122
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 718.61  on 653  degrees of freedom
## Residual deviance: 715.94  on 651  degrees of freedom
## AIC: 2223.9
## 
## Number of Fisher Scoring iterations: 5

4.4 Diễn giải mô hình GLM

\log(\mu_i) = \beta_0 + \beta_1 x_i \Rightarrow \mu_i = \exp(\beta_0 + \beta_1 x_i)

4.5 Ước lượng trung bình và sai số chuẩn

pred <- predict(glm.out, newdata = data.frame(x1 = 2, x2 = 3), se.fit = TRUE, type = "link")
fit.mu <- family(glm.out)$linkinv(pred$fit)
se.mu <- pred$se.fit * family(glm.out)$mu.eta(pred$fit)

4.6 Ước lượng \(\phi\) với quasi-likelihood

Dùng deviance hoặc Pearson.

4.7 Kiểm định Wald

Z = \frac{\hat{\beta}_j - \beta_j^0}{se(\hat{\beta}_j)} \sim N(0, 1)

4.8 Kiểm định LRT

# Mô hình nhỏ hơn (ít biến hơn)
glm.fit1 <- glm(y_pois ~ Age, family = poisson, data = lungcap)

# Mô hình đầy đủ hơn (nhiều biến hơn)
glm.fit2 <- glm(y_pois ~ Age + Ht, family = poisson, data = lungcap)

# So sánh hai mô hình bằng kiểm định Chi-squared
anova(glm.fit1, glm.fit2, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y_pois ~ Age
## Model 2: y_pois ~ Age + Ht
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       652     718.33                     
## 2       651     715.94  1   2.3911    0.122

4.9 Kiểm định Score

# Load thư viện
library(statmod)

## Warning: package 'statmod' was built under R version 4.4.3

# Mô hình null (chỉ có intercept)
glm.fit.null <- glm(y_pois ~ 1, family = poisson, data = lungcap)

# Tạo biến cần kiểm định thêm vào
x1 <- model.matrix(~ Age, data = lungcap)[, -1]  # bỏ intercept

# Thực hiện score test
glm.scoretest(glm.fit.null, x1)

## [1] -0.5271859

4.10 Khoảng tin cậy

confint(glm.out)

## Waiting for profiling to be done...

##                    2.5 %      97.5 %
## (Intercept) -0.553651655 0.923310593
## x1          -0.054745874 0.006317257
## x2          -0.003284268 0.027831950

4.11 Vùng tin cậy nhiều/tham số đơn

(\hat{\zeta} - \zeta)^T I(\hat{\zeta})(\hat{\zeta} - \zeta) \le \chi^2_{q, 1 - \alpha}

4.12 So sánh mô hình không lồng nhau

AIC(glm.out)

## [1] 2223.914

BIC(glm.out)

## [1] 2237.364

4.13 Tóm tắt chương

GLM dùng MLE để suy luận; dùng Wald, LRT, score để kiểm định; và AIC/BIC để so sánh mô hình.

CHƯƠNG 5: MÔ HÌNH POISON

5.1 Giới thiệu

Mô hình Poisson là một dạng cụ thể của mô hình tuyến tính tổng quát (GLM) dùng cho dữ liệu đếm. Biến phản hồi là số lần xảy ra của sự kiện trong một đơn vị thời gian hoặc không gian.

5.2 Phân phối Poisson

Giả định rằng: \[ Y_i \sim \text{Poisson}(\mu_i), \quad \mu_i = E(Y_i) \] Phân phối Poisson có kỳ vọng và phương sai bằng nhau, tức là \(Var(Y_i) = \mu_i\).

Hàm xác suất: \[ P(Y_i = y_i) = \frac{e^{-\mu_i} \mu_i^{y_i}}{y_i!} \] Trong đó: - \(y_i\): số lượng sự kiện quan sát được - \(\mu_i\): giá trị kỳ vọng của số sự kiện

5.3 Hàm liên kết

Mô hình Poisson sử dụng hàm liên kết log: \[ \log(\mu_i) = \eta_i = x_i^T \beta \] Tức là: \[ \mu_i = \exp(x_i^T \beta) \]

Giải thích: - \(x_i\): vector các biến giải thích cho quan sát thứ i - \(\beta\): vector hệ số hồi quy - \(\eta_i\): predictor tuyến tính

5.4 Ước lượng và suy luận

Các hệ số được ước lượng bằng phương pháp tối đa hóa hàm hợp lý (maximum likelihood). Mô hình được kiểm định qua z-test hoặc kiểm định deviance.

5.5 Ước lượng bằng R

Ví dụ với dữ liệu đếm:

data(AirPassengers)
counts <- as.numeric(AirPassengers)
time <- time(AirPassengers)
month <- cycle(AirPassengers)

poisson_model <- glm(counts ~ month, family = poisson(link = "log"))
summary(poisson_model)

## 
## Call:
## glm(formula = counts ~ month, family = poisson(link = "log"))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 5.584364   0.010733 520.280  < 2e-16 ***
## month       0.007865   0.001442   5.454 4.94e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 7216.8  on 143  degrees of freedom
## Residual deviance: 7187.1  on 142  degrees of freedom
## AIC: 8253.9
## 
## Number of Fisher Scoring iterations: 4

5.6 Diễn giải hệ số hồi quy

Các hệ số hồi quy được diễn giải theo log tỉ lệ xảy ra (log-rate). Ví dụ: - Nếu \(\beta_j = 0.3\), thì biến liên quan làm tăng tỉ lệ xảy ra trung bình khoảng \(\exp(0.3) \approx 1.35\) lần.

5.7 Độ phân tán

Nếu phương sai lớn hơn trung bình, ta gặp hiện tượng phân tán vượt mức (overdispersion). Khi đó, mô hình Poisson có thể không phù hợp vì giả định Var = Mean bị vi phạm.

5.8 Kiểm tra phân tán và điều chỉnh

Dùng thống kê Pearson hoặc deviance để kiểm tra overdispersion. Có thể dùng mô hình quasi-Poisson hoặc negative binomial để điều chỉnh.

poisson_quasi <- glm(counts ~ month, family = quasipoisson(link = "log"))
summary(poisson_quasi)

## 
## Call:
## glm(formula = counts ~ month, family = quasipoisson(link = "log"))
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.584364   0.076954  72.567   <2e-16 ***
## month       0.007865   0.010340   0.761    0.448    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 51.40374)
## 
##     Null deviance: 7216.8  on 143  degrees of freedom
## Residual deviance: 7187.1  on 142  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

5.9 Dự đoán

Sử dụng predict() để dự đoán số lượng sự kiện:

predict(poisson_model, newdata = data.frame(month = 6), type = "response")

##        1 
## 279.0956

5.10 Tóm tắt

• Mô hình Poisson dùng để phân tích dữ liệu đếm
• Hàm liên kết log giúp chuyển đổi về tuyến tính
• Cần kiểm tra hiện tượng phân tán để chọn mô hình phù hợp hơn nếu cần

CHƯƠNG 6.MÔ HÌNH NHỊ THỨC ÂM

6.1 Giới thiệu

Mô hình nhị thức âm (Negative Binomial - NB) là một mở rộng của mô hình Poisson để xử lý hiện tượng phân tán vượt mức (overdispersion), tức là khi phương sai lớn hơn kỳ vọng.

6.2 Phân phối nhị thức âm

Giả định: \[ Y_i \sim NB(\mu_i, \theta) \] Trong đó: - \(\mu_i\): kỳ vọng của \(Y_i\) - \(\theta\): tham số phân tán (dispersion)

Phương sai: \[ Var(Y_i) = \mu_i + \frac{\mu_i^2}{\theta} \] Phân phối NB cho phép phương sai lớn hơn kỳ vọng và linh hoạt hơn so với Poisson.

6.3 Hàm liên kết

Hàm liên kết thường dùng là log: \[ \log(\mu_i) = x_i^T \beta \] Tương tự mô hình Poisson, mô hình NB thuộc họ GLM với hàm liên kết log.

6.4 Ước lượng và suy luận

Ước lượng các hệ số \(\beta\) và tham số phân tán \(\theta\) bằng phương pháp tối đa hóa hàm hợp lý. Việc ước lượng \(\theta\) thường dùng thuật toán lặp hoặc gói chuyên biệt.

6.5 Mô hình NB trong R

Dùng gói MASS và hàm glm.nb():

library(MASS)
data(Cars93, package = "MASS")
nb_model <- glm.nb(Price ~ Horsepower + Weight, data = Cars93)

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 29.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 37.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 40.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 29.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 25.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 7.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 47.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 35.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 36.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 32.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 61.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 21.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 24.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 28.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 22.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 22.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.700000

summary(nb_model)

## 
## Call:
## glm.nb(formula = Price ~ Horsepower + Weight, data = Cars93, 
##     init.theta = 38.78246161, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.421e+00  1.749e-01   8.126 4.45e-16 ***
## Horsepower  5.025e-03  7.254e-04   6.927 4.29e-12 ***
## Weight      2.455e-04  7.173e-05   3.423 0.000619 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(38.7825) family taken to be 1)
## 
##     Null deviance: 253.205  on 92  degrees of freedom
## Residual deviance:  77.021  on 90  degrees of freedom
## AIC: 560.3
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  38.8 
##           Std. Err.:  14.6 
## 
##  2 x log-likelihood:  -552.3

6.6 So sánh với mô hình Poisson

Khi có overdispersion, mô hình NB thường cho kết quả phù hợp hơn và sai số chuẩn chính xác hơn.

poisson_model <- glm(Price ~ Horsepower + Weight, data = Cars93, family = poisson())

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 29.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 37.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 40.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 29.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 25.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 7.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 12.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 47.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 35.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 34.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 36.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 32.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 31.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 61.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 15.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 21.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 13.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 16.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 20.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 14.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 17.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 24.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 28.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 11.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 10.900000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.500000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 8.600000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.800000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.400000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 18.200000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 22.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 9.100000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 19.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 23.300000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 22.700000

## Warning in dpois(y, mu, log = TRUE): non-integer x = 26.700000

dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / df.residual(poisson_model)
dispersion  # Nếu > 1 đáng kể → nên dùng NB

## [1] 1.572621

6.7 Dự đoán

predict(nb_model, newdata = data.frame(Horsepower = 150, Weight = 3000), type = "response")

##        1 
## 18.37969

Dự đoán trả về kỳ vọng số sự kiện xảy ra (giá trị trung bình của \(Y\)).

6.8 Tóm tắt

• Mô hình nhị thức âm mở rộng Poisson để xử lý phân tán vượt mức
• Phương sai phụ thuộc vào cả \(\mu\) và tham số \(\theta\)
• R hỗ trợ mô hình này thông qua hàm glm.nb() từ gói MASS

CHƯƠNG 7: CHẨN ĐOÁN VÀ KIỂM TRA MÔ HÌNH

7.1 Kiểm tra phần dư (Residuals)

Phần dư giúp đánh giá sự phù hợp của mô hình. Trong GLM, phần dư Pearson và phần dư deviance được sử dụng phổ biến.

• Residual deviance: so sánh log-likelihood giữa mô hình đã fit và mô hình hoàn hảo.

model <- glm(y_pois ~ Age + Ht, family = poisson, data = lungcap)
residuals(model, type = "deviance")[1:5]  # ví dụ một số phần dư

##           1           2           3           4           5 
## -0.76444044 -0.76521674  0.02014332  1.26597499 -0.78083895

7.2 Biểu đồ chuẩn đoán

Các biểu đồ thường dùng: - Residuals vs Fitted - Normal QQ plot - Cook’s distance

plot(model)  # tạo 4 biểu đồ chẩn đoán cơ bản

7.3 Leverage và ảnh hưởng

Leverage đo lường ảnh hưởng của điểm dữ liệu đến mô hình. Điểm có leverage cao và phần dư lớn có thể là outlier ảnh hưởng.

hatvalues(model)[1:5]  # các giá trị leverage đầu tiên

##           1           2           3           4           5 
## 0.012439811 0.009705965 0.009705965 0.009705965 0.008929733

7.4 Cook’s distance

Đo mức ảnh hưởng của từng quan sát đến toàn bộ mô hình.

cooks.distance(model)[1:5]  # Cook's distance đầu tiên

##            1            2            3            4            5 
## 2.032083e-03 1.579711e-03 1.345002e-06 6.884086e-03 1.506144e-03

7.5 DFBETAs và DFFITS

Giúp đánh giá ảnh hưởng của từng điểm đến từng hệ số ước lượng.

dfbeta(model)[1:5, ]

##     (Intercept)           Age            Ht
## 1 -0.0216342892  1.940289e-04  3.085303e-04
## 2 -0.0192458365  1.432404e-04  2.777418e-04
## 3  0.0005066211 -3.770613e-06 -7.311185e-06
## 4  0.0318403224 -2.369770e-04 -4.594962e-04
## 5 -0.0165293514  2.560891e-04  2.147991e-04

dffits(model)[1:5]

##            1            2            3            4            5 
## -0.082297165 -0.072566708  0.001909434  0.120140683 -0.070971060

7.6 Kiểm định và phân tích mô hình thay thế

Dùng so sánh mô hình (ANOVA) hoặc kiểm định score/likelihood để xem biến có nên đưa vào không.

model0 <- glm(y_pois ~ 1, family = poisson, data = lungcap)
model1 <- glm(y_pois ~ Age, family = poisson, data = lungcap)
anova(model0, model1, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y_pois ~ 1
## Model 2: y_pois ~ Age
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       653     718.61                     
## 2       652     718.33  1  0.27848   0.5977

CHƯƠNG 8: PHÂN TÁN QUÁ MỨC VÀ HÀM QUASI-LIKELIHOOD

8.1 Hiện tượng overdispersion

Overdispersion xảy ra khi phương sai lớn hơn kỳ vọng (Poisson giả định Var = Mean). Dễ gặp trong dữ liệu đếm.

Kiểm tra nhanh:

dispersion <- sum(residuals(model, type = "pearson")^2) / model$df.residual
dispersion  # Nếu > 1 nhiều thì có overdispersion

## [1] 0.9811308

8.2 Quasi-likelihood

Sử dụng khi không muốn giả định phân phối xác định nhưng vẫn mô hình hóa kỳ vọng và phương sai.

model_quasi <- glm(y_pois ~ Age + Ht, family = quasipoisson, data = lungcap)
summary(model_quasi)

## 
## Call:
## glm(formula = y_pois ~ Age + Ht, family = quasipoisson, data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.185413   0.373197   0.497    0.619
## Age         -0.023939   0.015430  -1.552    0.121
## Ht           0.012275   0.007862   1.561    0.119
## 
## (Dispersion parameter for quasipoisson family taken to be 0.9811308)
## 
##     Null deviance: 718.61  on 653  degrees of freedom
## Residual deviance: 715.94  on 651  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

8.3 Negative binomial

Thay vì dùng Poisson, có thể dùng phân phối negative binomial để xử lý overdispersion.

library(MASS)
model_nb <- glm.nb(y_pois ~ Age + Ht, data = lungcap)

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

summary(model_nb)

## 
## Call:
## glm.nb(formula = y_pois ~ Age + Ht, data = lungcap, init.theta = 13311.2109, 
##     link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.185408   0.376797   0.492    0.623
## Age         -0.023939   0.015578  -1.537    0.124
## Ht           0.012275   0.007938   1.546    0.122
## 
## (Dispersion parameter for Negative Binomial(13311.21) family taken to be 1)
## 
##     Null deviance: 718.51  on 653  degrees of freedom
## Residual deviance: 715.84  on 651  degrees of freedom
## AIC: 2225.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  13311 
##           Std. Err.:  149147 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -2217.917

8.4 So sánh mô hình

So sánh Poisson, Quasipoisson, và Negative Binomial để chọn mô hình tốt nhất.

AIC(model, model_nb)  # AIC nhỏ hơn là mô hình tốt hơn

##          df      AIC
## model     3 2223.914
## model_nb  4 2225.917

8.5 Lý do sử dụng quasi-likelihood

Giúp ước lượng độ lệch chuẩn đúng hơn khi không biết rõ phân phối, đặc biệt hữu ích trong dữ liệu sinh học, xã hội.

Ghi chú: Các mô hình cần kiểm tra thêm giả định để đảm bảo ý nghĩa thống kê và sự phù hợp.

CHƯƠNG 9.MÔ HÌNH CẤU TRÚC

9.1 Giới thiệu

Mô hình có cấu trúc là các mô hình mở rộng của GLM trong đó có mối liên hệ nội tại giữa các tham số hoặc dữ liệu. Các mô hình này được sử dụng khi dữ liệu có cấu trúc không độc lập như: lặp lại theo thời gian, theo cụm, hoặc phân cấp.

9.2 Dữ liệu phân cấp và lặp lại

Trong dữ liệu có cấu trúc phân cấp (hierarchical) hoặc lặp lại (repeated measures), các quan sát không độc lập hoàn toàn. Ví dụ: - Đo huyết áp nhiều lần trên cùng một bệnh nhân - Sinh viên trong cùng một lớp học

9.3 Mô hình hỗn hợp tuyến tính (Linear Mixed Models)

Mô hình hỗn hợp tuyến tính thêm các hiệu ứng ngẫu nhiên (random effects) vào mô hình tuyến tính:

\[ Y_{ij} = \beta_0 + \beta_1 x_{ij} + b_j + \varepsilon_{ij} \] Trong đó: - \(b_j \sim N(0, \sigma_b^2)\): hiệu ứng ngẫu nhiên của nhóm j - \(\varepsilon_{ij} \sim N(0, \sigma^2)\): sai số ngẫu nhiên

9.4 Ước lượng trong mô hình hỗn hợp

Các tham số được ước lượng bằng phương pháp REML (Restricted Maximum Likelihood) hoặc ML (Maximum Likelihood).

Trong R, dùng hàm lmer() từ gói lme4:

library(lme4)

## Warning: package 'lme4' was built under R version 4.4.3

## Loading required package: Matrix

model_lmm <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
summary(model_lmm)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

9.5 Mô hình tuyến tính tổng quát hỗn hợp (GLMM)

GLMM kết hợp giữa GLM và hiệu ứng ngẫu nhiên, thích hợp cho dữ liệu đếm hoặc nhị phân có cấu trúc:

Ví dụ:

data(cbpp, package = "lme4")
glmm_model <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), 
                    family = binomial, data = cbpp)
summary(glmm_model)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(incidence, size - incidence) ~ period + (1 | herd)
##    Data: cbpp
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     194.1     204.2     -92.0     184.1        51 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3816 -0.7889 -0.2026  0.5142  2.8791 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  herd   (Intercept) 0.4123   0.6421  
## Number of obs: 56, groups:  herd, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3983     0.2312  -6.048 1.47e-09 ***
## period2      -0.9919     0.3032  -3.272 0.001068 ** 
## period3      -1.1282     0.3228  -3.495 0.000474 ***
## period4      -1.5797     0.4220  -3.743 0.000182 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) perid2 perid3
## period2 -0.363              
## period3 -0.340  0.280       
## period4 -0.260  0.213  0.198

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

Hệ số cố định (fixed effects) được diễn giải như GLM. Các hiệu ứng ngẫu nhiên phản ánh sự biến thiên giữa các nhóm hoặc đơn vị.

9.7 So sánh mô hình có cấu trúc và GLM

GLMM và LMM phù hợp hơn GLM khi dữ liệu có cấu trúc nhóm, thời gian, hoặc đo lặp. Dùng AIC hoặc kiểm định likelihood ratio để so sánh mô hình.

9.8 Tóm tắt

• Mô hình có cấu trúc xử lý dữ liệu phụ thuộc hoặc phân nhóm
• Bao gồm LMM (dữ liệu liên tục) và GLMM (dữ liệu đếm, nhị phân)
• R hỗ trợ ước lượng qua lme4::lmer() và lme4::glmer()

CHƯƠNG 10: MÔ HÌNH CHO DỮ LIỆU NHỊ PHÂN

10.1 Dữ liệu nhị phân và hồi quy logistic

Dữ liệu nhị phân gồm hai giá trị (thành công/thất bại, 1/0). Hồi quy logistic dùng để mô hình xác suất thành công theo các biến giải thích.

Hàm logit: \[ \text{logit}(p) = \log\left(\frac{p}{1 - p}\right) \]

model_logit <- glm(y_bin ~ Age + Ht, family = binomial(link = "logit"), data = lungcap)

## Warning: glm.fit: algorithm did not converge

summary(model_logit)

## 
## Call:
## glm(formula = y_bin ~ Age + Ht, family = binomial(link = "logit"), 
##     data = lungcap)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.657e+01  1.903e+05       0        1
## Age          4.932e-15  7.727e+03       0        1
## Ht           8.500e-15  4.002e+03       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.0000e+00  on 653  degrees of freedom
## Residual deviance: 3.7942e-09  on 651  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

10.2 Hàm liên kết khác

Ngoài logit, còn có probit và cloglog:

$$ model_probit <- glm(y_bin ~ Age + Ht, family = binomial(link = “probit”), data = lungcap)\

model_cloglog <- glm(y_bin ~ Age + Ht, family = binomial(link = “cloglog”), data = lungcap) $$

10.3 Diễn giải hệ số

Hệ số hồi quy logistic biểu diễn log odds ratio. Để diễn giải dễ hơn, thường chuyển sang odds ratio:

exp(coef(model_logit))  # odds ratio

##  (Intercept)          Age           Ht 
## 2.900701e-12 1.000000e+00 1.000000e+00

10.4 Mô hình với biến phân loại

Nếu có biến phân loại như Gender hoặc Smoke, cần đảm bảo nó là factor:

lungcap$Gender <- as.factor(lungcap$Gender)
lungcap$Smoke <- as.factor(lungcap$Smoke)
model_cat <- glm(y_bin ~ Gender + Smoke, family = binomial, data = lungcap)

## Warning: glm.fit: algorithm did not converge

summary(model_cat)

## 
## Call:
## glm(formula = y_bin ~ Gender + Smoke, family = binomial, data = lungcap)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.657e+01  2.077e+04  -0.001    0.999
## GenderM      8.477e-14  2.794e+04   0.000    1.000
## Smoke1       6.212e-14  4.668e+04   0.000    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 0.0000e+00  on 653  degrees of freedom
## Residual deviance: 3.7942e-09  on 651  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

10.5 Kiểm định và độ phù hợp

Kiểm định tổng thể bằng deviance:

dev <- deviance(model_logit)
df <- df.residual(model_logit)
pchisq(dev, df, lower.tail = FALSE)  # p-value

## [1] 1

Ngoài ra có thể dùng Hosmer-Lemeshow test (trong gói ResourceSelection).

10.6 Đánh giá mô hình

Dùng ROC curve, AUC, hoặc bảng phân loại:

pred <- predict(model_logit, type = "response")
table(Predicted = pred > 0.5, Actual = lungcap$y_bin)

##          Actual
## Predicted   0
##     FALSE 654

Ghi chú: Mô hình nhị phân là một trong những ứng dụng phổ biến nhất của GLM trong khoa học xã hội, y tế và kinh tế lượng.

CHƯƠNG 11.DỮ LIỆU LIÊN TỤC DƯƠNG: GLM VỚI PHÂN PHỐI GAMMA VÀ INVERSE GAUSIAN

11.1 Giới thiệu

Dữ liệu liên tục dương phổ biến trong nhiều lĩnh vực (sinh học, kinh tế, kỹ thuật). Hai mô hình GLM thường dùng là: - Gamma GLM: phù hợp với dữ liệu có phương sai tăng theo bình phương trung bình. - Inverse Gaussian GLM: phù hợp khi dữ liệu lệch phải nhiều, phương sai tăng theo lập phương trung bình.

11.2 Mô hình hóa dữ liệu liên tục dương

Dữ liệu như FEV (Forced Expiratory Volume) là liên tục dương, cần mô hình phản ánh tính chất phân phối và mối quan hệ mean–variance.

11.3 Phân phối Gamma

• Hàm phương sai: V(μ) = μ²
• Dữ liệu có phân phối lệch nhẹ
• Thường dùng link "log"

library(GLMsData)
data(lungcap)
gamma_model <- glm(FEV ~ Age + Ht + Gender + Smoke,
                   family = Gamma(link = "log"),
                   data = lungcap)
summary(gamma_model)

## 
## Call:
## glm(formula = FEV ~ Age + Ht + Gender + Smoke, family = Gamma(link = "log"), 
##     data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.940198   0.076402 -25.395  < 2e-16 ***
## Age          0.022816   0.003253   7.013 5.86e-12 ***
## Ht           0.042985   0.001631  26.351  < 2e-16 ***
## GenderM      0.029409   0.011385   2.583   0.0100 *  
## Smoke       -0.040853   0.020315  -2.011   0.0447 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01997449)
## 
##     Null deviance: 70.791  on 653  degrees of freedom
## Residual deviance: 13.414  on 649  degrees of freedom
## AIC: 525.61
## 
## Number of Fisher Scoring iterations: 4

11.4 Phân phối Inverse Gaussian

• Hàm phương sai: V(μ) = μ³
• Mô hình tốt khi dữ liệu lệch mạnh
• Canonical link: 1/μ², nhưng "log" được dùng phổ biến

invgauss_model <- glm(FEV ~ Age + Ht + Gender + Smoke,
                      family = inverse.gaussian(link = "log"),
                      data = lungcap)
summary(invgauss_model)

## 
## Call:
## glm(formula = FEV ~ Age + Ht + Gender + Smoke, family = inverse.gaussian(link = "log"), 
##     data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.964597   0.076821 -25.574  < 2e-16 ***
## Age          0.021928   0.003532   6.209 9.54e-10 ***
## Ht           0.043535   0.001701  25.597  < 2e-16 ***
## GenderM      0.028638   0.011243   2.547   0.0111 *  
## Smoke       -0.039314   0.023159  -1.698   0.0901 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.008309335)
## 
##     Null deviance: 29.0674  on 653  degrees of freedom
## Residual deviance:  5.8086  on 649  degrees of freedom
## AIC: 574.71
## 
## Number of Fisher Scoring iterations: 4

11.5 Link function

• log: phổ biến và đảm bảo đầu ra dương
• inverse: dùng cho các mối quan hệ phi tuyến
• identity: chỉ dùng khi biến phản hồi không bị giới hạn
• 1/μ²: link chuẩn cho Inverse Gaussian

11.6 Ước lượng tham số phân tán φ

Dùng residual deviance hoặc Pearson residuals.

phi_gamma <- sum(residuals(gamma_model, type="pearson")^2) / df.residual(gamma_model)
phi_ig <- sum(residuals(invgauss_model, type="pearson")^2) / df.residual(invgauss_model)
phi_gamma; phi_ig

## [1] 0.01997449

## [1] 0.008309329

11.7 Nghiên cứu tình huống

Sách trình bày 2 case study, không dùng lungcap, nhưng quy trình áp dụng tương tự: - Mô hình hóa - Chẩn đoán - So sánh fitted values, residuals

11.8 Dùng R để fit mô hình

• glm(..., family = Gamma(link = "log"))
• glm(..., family = inverse.gaussian(link = "log"))

11.9 Tóm tắt

Mô hình	Hàm phương sai	Link phổ biến	Dữ liệu phù hợp
Gamma	V(μ) = μ²	log	Dữ liệu lệch nhẹ
Inverse Gaussian	V(μ) = μ³	log, 1/μ²	Dữ liệu lệch mạnh, thời gian

Tạo đồ thị chuẩn đoán dựa trên dataset trong sách

Nạp dữ liệu

library(GLMsData)
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.4.3

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

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

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

data(lungcap)
head(lungcap)

##   Age   FEV Ht Gender Smoke
## 1   3 1.072 46      F     0
## 2   4 0.839 48      F     0
## 3   4 1.102 48      F     0
## 4   4 1.389 48      F     0
## 5   4 1.577 49      F     0
## 6   4 1.418 49      F     0

Mô hình Gamma (log link)

gamma_model <- glm(FEV ~ Age + Ht + Gender + Smoke,
                   family = Gamma(link = "log"),
                   data = lungcap)
summary(gamma_model)

## 
## Call:
## glm(formula = FEV ~ Age + Ht + Gender + Smoke, family = Gamma(link = "log"), 
##     data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.940198   0.076402 -25.395  < 2e-16 ***
## Age          0.022816   0.003253   7.013 5.86e-12 ***
## Ht           0.042985   0.001631  26.351  < 2e-16 ***
## GenderM      0.029409   0.011385   2.583   0.0100 *  
## Smoke       -0.040853   0.020315  -2.011   0.0447 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01997449)
## 
##     Null deviance: 70.791  on 653  degrees of freedom
## Residual deviance: 13.414  on 649  degrees of freedom
## AIC: 525.61
## 
## Number of Fisher Scoring iterations: 4

Mô hình Inverse Gaussian (log link)

invgauss_model <- glm(FEV ~ Age + Ht + Gender + Smoke,
                      family = inverse.gaussian(link = "log"),
                      data = lungcap)
summary(invgauss_model)

## 
## Call:
## glm(formula = FEV ~ Age + Ht + Gender + Smoke, family = inverse.gaussian(link = "log"), 
##     data = lungcap)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.964597   0.076821 -25.574  < 2e-16 ***
## Age          0.021928   0.003532   6.209 9.54e-10 ***
## Ht           0.043535   0.001701  25.597  < 2e-16 ***
## GenderM      0.028638   0.011243   2.547   0.0111 *  
## Smoke       -0.039314   0.023159  -1.698   0.0901 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.008309335)
## 
##     Null deviance: 29.0674  on 653  degrees of freedom
## Residual deviance:  5.8086  on 649  degrees of freedom
## AIC: 574.71
## 
## Number of Fisher Scoring iterations: 4

Chẩn đoán mô hình Gamma

par(mfrow = c(2, 2))
plot(gamma_model)

Chẩn đoán mô hình Inverse Gaussian

par(mfrow = c(2, 2))
plot(invgauss_model)

Residuals vs Fitted cho Gamma

lungcap$gamma_resid <- residuals(gamma_model, type = "pearson")
lungcap$gamma_fitted <- fitted(gamma_model)

ggplot(lungcap, aes(x = gamma_fitted, y = gamma_resid)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Residuals vs Fitted (Gamma)", x = "Fitted values", y = "Pearson Residuals")

Residuals vs Fitted cho Inverse Gaussian

lungcap$ig_resid <- residuals(invgauss_model, type = "pearson")
lungcap$ig_fitted <- fitted(invgauss_model)

ggplot(lungcap, aes(x = ig_fitted, y = ig_resid)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Residuals vs Fitted (Inverse Gaussian)", x = "Fitted values", y = "Pearson Residuals")

CHƯƠNG 12: MÔ HÌNH CÓ HIỆU ỨNG NGẪU NHIÊN (RANDOM EFFECTS)

12.1 Mô hình hỗn hợp tuyến tính tổng quát (GLMM)

Khi dữ liệu có cấu trúc nhóm (ví dụ: bệnh nhân trong các bệnh viện), cần đưa hiệu ứng ngẫu nhiên vào mô hình để phản ánh sự phụ thuộc nội tại.

library(lme4)
lungcap$Ht_grp <- cut(lungcap$Ht, breaks = 3, labels = c("Low", "Medium", "High"))
library(GLMsData)
library(lme4)

data(lungcap)

# Tạo biến nhị phân
lungcap$y_bin <- ifelse(lungcap$FEV > 2.5, 1, 0)

# Mô hình GLMM với random intercept theo chiều cao (Ht)
model_glmm <- glmer(y_bin ~ Age + (1 | Ht), family = binomial, data = lungcap)
summary(model_glmm)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y_bin ~ Age + (1 | Ht)
##    Data: lungcap
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     507.1     520.6    -250.6     501.1       651 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.6916 -0.2560  0.0624  0.3733  4.6827 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ht     (Intercept) 3.515    1.875   
## Number of obs: 654, groups:  Ht, 56
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.65530    0.87149  -6.489 8.63e-11 ***
## Age          0.56599    0.08892   6.365 1.95e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age -0.926

12.2 Cú pháp mô hình hỗn hợp

• (1 | Group) chỉ ra intercept thay đổi theo nhóm.
• (Age | Group) cho phép slope và intercept thay đổi theo nhóm.

12.3 So sánh GLM và GLMM

GLMM cho phép xử lý dữ liệu có phụ thuộc nội tại (clustered/correlated). GLM giả định độc lập giữa các quan sát, nên không phù hợp khi có cấu trúc nhóm rõ ràng.

12.4 Ước lượng và hội tụ

GLMM thường dùng phương pháp tối đa hóa hợp lý xấp xỉ như Laplace hoặc quadrature. Việc hội tụ có thể chậm hoặc thất bại nếu mô hình quá phức tạp.

model_glmm2 <- glmer(y_bin ~ Age + (Age | Ht), family = binomial, data = lungcap)

KẾT LUẬN TỔNG QUAN VỀ QUYỂN SÁCH

Quyển “Generalized Linear Models with Examples in R” mang lại một cái nhìn toàn diện và thực hành về mô hình tuyến tính tổng quát (GLM). Các nội dung chính bao gồm:

• Giới thiệu lý thuyết GLM: liên kết, phân phối thuộc họ hàm mũ, hàm tuyến tính.
• Các mô hình quan trọng: hồi quy logistic, hồi quy Poisson, quasi-likelihood và negative binomial.
• Kiểm tra độ phù hợp mô hình và các công cụ chẩn đoán phần dư, ảnh hưởng.
• Mở rộng GLM với mô hình có hiệu ứng ngẫu nhiên (GLMM).

Điểm nổi bật là cách trình bày thực tế, tích hợp ví dụ cụ thể bằng ngôn ngữ R, giúp người học vừa hiểu lý thuyết vừa áp dụng thực hành. Đây là tài liệu quan trọng cho những ai làm việc với dữ liệu định lượng trong các lĩnh vực sinh học, y tế, xã hội, và kinh tế lượng.

Ghi chú cuối cùng: Mặc dù GLM là một công cụ mạnh mẽ, người dùng cần kiểm tra cẩn thận các giả định, cấu trúc dữ liệu và phù hợp mô hình để đưa ra suy luận chính xác.

TASK 2.THỐNG KÊ MÔ TẢ “SUPERMARKET TRANSACTIONS”

1. Đọc dữ liệu

data <- read.csv("C:/Users/Admin/Downloads/Supermarket Transactions.csv")
str(data)

## 'data.frame':    14059 obs. of  16 variables:
##  $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ PurchaseDate     : chr  "2007-12-18" "2007-12-20" "2007-12-21" "2007-12-21" ...
##  $ CustomerID       : int  7223 7841 8374 9619 1900 6696 9673 354 1293 7938 ...
##  $ Gender           : chr  "F" "M" "F" "M" ...
##  $ MaritalStatus    : chr  "S" "M" "M" "M" ...
##  $ Homeowner        : chr  "Y" "Y" "N" "Y" ...
##  $ Children         : int  2 5 2 3 3 3 2 2 3 1 ...
##  $ AnnualIncome     : chr  "$30K - $50K" "$70K - $90K" "$50K - $70K" "$30K - $50K" ...
##  $ City             : chr  "Los Angeles" "Los Angeles" "Bremerton" "Portland" ...
##  $ StateorProvince  : chr  "CA" "CA" "WA" "OR" ...
##  $ Country          : chr  "USA" "USA" "USA" "USA" ...
##  $ ProductFamily    : chr  "Food" "Food" "Food" "Food" ...
##  $ ProductDepartment: chr  "Snack Foods" "Produce" "Snack Foods" "Snacks" ...
##  $ ProductCategory  : chr  "Snack Foods" "Vegetables" "Snack Foods" "Candy" ...
##  $ UnitsSold        : int  5 5 3 4 4 3 4 6 1 2 ...
##  $ Revenue          : num  27.38 14.9 5.52 4.44 14 ...

2. Tổng quan dữ liệu

summary(data)

##        X         PurchaseDate         CustomerID       Gender         
##  Min.   :    1   Length:14059       Min.   :    3   Length:14059      
##  1st Qu.: 3516   Class :character   1st Qu.: 2549   Class :character  
##  Median : 7030   Mode  :character   Median : 5060   Mode  :character  
##  Mean   : 7030                      Mean   : 5117                     
##  3rd Qu.:10544                      3rd Qu.: 7633                     
##  Max.   :14059                      Max.   :10280                     
##  MaritalStatus       Homeowner            Children    AnnualIncome      
##  Length:14059       Length:14059       Min.   :0.00   Length:14059      
##  Class :character   Class :character   1st Qu.:1.00   Class :character  
##  Mode  :character   Mode  :character   Median :3.00   Mode  :character  
##                                        Mean   :2.53                     
##                                        3rd Qu.:4.00                     
##                                        Max.   :5.00                     
##      City           StateorProvince      Country          ProductFamily     
##  Length:14059       Length:14059       Length:14059       Length:14059      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##  ProductDepartment  ProductCategory      UnitsSold        Revenue     
##  Length:14059       Length:14059       Min.   :1.000   Min.   : 0.53  
##  Class :character   Class :character   1st Qu.:3.000   1st Qu.: 6.84  
##  Mode  :character   Mode  :character   Median :4.000   Median :11.25  
##                                        Mean   :4.081   Mean   :13.00  
##                                        3rd Qu.:5.000   3rd Qu.:17.37  
##                                        Max.   :8.000   Max.   :56.70

3. Thống kê mô tả các biến định lượng

data %>% 
  select(AnnualIncome, Children, UnitsSold, Revenue) %>%
  summary()

##  AnnualIncome          Children      UnitsSold        Revenue     
##  Length:14059       Min.   :0.00   Min.   :1.000   Min.   : 0.53  
##  Class :character   1st Qu.:1.00   1st Qu.:3.000   1st Qu.: 6.84  
##  Mode  :character   Median :3.00   Median :4.000   Median :11.25  
##                     Mean   :2.53   Mean   :4.081   Mean   :13.00  
##                     3rd Qu.:4.00   3rd Qu.:5.000   3rd Qu.:17.37  
##                     Max.   :5.00   Max.   :8.000   Max.   :56.70

4. Thống kê mô tả các biến định tính

categorical_vars <- c("Gender", "MaritalStatus", "Homeowner", "City", "StateorProvince", 
                      "Country", "ProductFamily", "ProductDepartment", "ProductCategory")

for (var in categorical_vars) {
  cat("\n\n###", var, "\n")
  print(table(data[[var]]))
}

## 
## 
## ### Gender 
## 
##    F    M 
## 7170 6889 
## 
## 
## ### MaritalStatus 
## 
##    M    S 
## 6866 7193 
## 
## 
## ### Homeowner 
## 
##    N    Y 
## 5615 8444 
## 
## 
## ### City 
## 
##      Acapulco    Bellingham Beverly Hills     Bremerton       Camacho 
##           383           143           811           834           452 
##   Guadalajara       Hidalgo   Los Angeles        Merida   Mexico City 
##            75           845           926           654           194 
##       Orizaba      Portland         Salem    San Andres     San Diego 
##           464           876          1386           621           866 
## San Francisco       Seattle       Spokane        Tacoma     Vancouver 
##           130           922           875          1257           633 
##      Victoria   Walla Walla        Yakima 
##           176           160           376 
## 
## 
## ### StateorProvince 
## 
##        BC        CA        DF  Guerrero   Jalisco        OR  Veracruz        WA 
##       809      2733       815       383        75      2262       464      4567 
##   Yucatan Zacatecas 
##       654      1297 
## 
## 
## ### Country 
## 
## Canada Mexico    USA 
##    809   3688   9562 
## 
## 
## ### ProductFamily 
## 
##          Drink           Food Non-Consumable 
##           1250          10153           2656 
## 
## 
## ### ProductDepartment 
## 
## Alcoholic Beverages         Baked Goods        Baking Goods           Beverages 
##                 356                 425                1072                 680 
##     Breakfast Foods        Canned Foods     Canned Products            Carousel 
##                 188                 977                 109                  59 
##            Checkout               Dairy                Deli                Eggs 
##                  82                 903                 699                 198 
##        Frozen Foods  Health and Hygiene           Household                Meat 
##                1382                 893                1420                  89 
##         Periodicals             Produce             Seafood         Snack Foods 
##                 202                1994                 102                1600 
##              Snacks       Starchy Foods 
##                 352                 277 
## 
## 
## ### ProductCategory 
## 
##         Baking Goods    Bathroom Products        Beer and Wine 
##                  484                  365                  356 
##                Bread      Breakfast Foods              Candles 
##                  425                  417                   45 
##                Candy     Canned Anchovies         Canned Clams 
##                  352                   44                   53 
##       Canned Oysters      Canned Sardines        Canned Shrimp 
##                   35                   40                   38 
##          Canned Soup          Canned Tuna Carbonated Beverages 
##                  404                   87                  154 
##    Cleaning Supplies        Cold Remedies                Dairy 
##                  189                   93                  903 
##        Decongestants               Drinks                 Eggs 
##                   85                  135                  198 
##           Electrical      Frozen Desserts       Frozen Entrees 
##                  355                  323                  118 
##                Fruit             Hardware        Hot Beverages 
##                  765                  129                  226 
##              Hygiene     Jams and Jellies     Kitchen Products 
##                  197                  588                  217 
##            Magazines                 Meat        Miscellaneous 
##                  202                  761                   42 
##  Packaged Vegetables       Pain Relievers       Paper Products 
##                   48                  192                  345 
##                Pizza     Plastic Products Pure Juice Beverages 
##                  194                  141                  165 
##              Seafood          Side Dishes          Snack Foods 
##                  102                  153                 1600 
##            Specialty        Starchy Foods           Vegetables 
##                  289                  277                 1728

5. Phân phối giới tính và tình trạng hôn nhân

table(data$Gender)

## 
##    F    M 
## 7170 6889

table(data$MaritalStatus)

## 
##    M    S 
## 6866 7193

6. Số lượng khách hàng theo thành phố (Top 10)

data %>%
  count(City, sort = TRUE) %>%
  head(10)

##             City    n
## 1          Salem 1386
## 2         Tacoma 1257
## 3    Los Angeles  926
## 4        Seattle  922
## 5       Portland  876
## 6        Spokane  875
## 7      San Diego  866
## 8        Hidalgo  845
## 9      Bremerton  834
## 10 Beverly Hills  811

7. Tổng doanh thu theo bang (Top 10)

data %>%
  group_by(StateorProvince) %>%
  summarise(TongDoanhThu = sum(Revenue, na.rm = TRUE)) %>%
  arrange(desc(TongDoanhThu)) %>%
  head(10)

## # A tibble: 10 × 2
##    StateorProvince TongDoanhThu
##    <chr>                  <dbl>
##  1 WA                    58420.
##  2 CA                    34730.
##  3 OR                    30138.
##  4 Zacatecas             17110.
##  5 BC                    11067.
##  6 DF                    10694.
##  7 Yucatan                8740.
##  8 Veracruz               6245.
##  9 Guerrero               5161.
## 10 Jalisco                 523.

8. Doanh thu theo nhóm sản phẩm

data %>%
  group_by(ProductFamily) %>%
  summarise(TongDoanhThu = sum(Revenue, na.rm = TRUE)) %>%
  arrange(desc(TongDoanhThu))

## # A tibble: 3 × 2
##   ProductFamily  TongDoanhThu
##   <chr>                 <dbl>
## 1 Food                132761.
## 2 Non-Consumable       34196.
## 3 Drink                15874.

9. Biểu đồ doanh thu theo phòng ban sản phẩm

data %>%
  group_by(ProductDepartment) %>%
  summarise(Revenue = sum(Revenue, na.rm = TRUE)) %>%
  ggplot(aes(x = reorder(ProductDepartment, Revenue), y = Revenue)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "Doanh thu theo phòng ban sản phẩm", x = "Phòng ban", y = "Doanh thu")

10. Kết luận

Dữ liệu giao dịch siêu thị cung cấp thông tin đa dạng từ thông tin khách hàng đến thông tin sản phẩm và doanh thu. Qua phân tích mô tả:

• Biến định lượng như thu nhập, số lượng bán và doanh thu có độ phân tán lớn, cho thấy sự khác biệt đáng kể giữa các khách hàng và sản phẩm.
• Biến định tính như giới tính, tình trạng hôn nhân và nhóm sản phẩm cho thấy sự đa dạng trong đặc điểm khách hàng và danh mục sản phẩm.
• Doanh thu tập trung nhiều ở một số bang và nhóm sản phẩm cụ thể, gợi ý tiềm năng trong việc phân khúc thị trường và tối ưu hóa danh mục sản phẩm.

Phân tích mô tả này là bước đầu quan trọng giúp hiểu rõ cấu trúc dữ liệu, làm cơ sở cho các phân tích sâu hơn như phân cụm khách hàng, dự đoán doanh thu hoặc phân tích hành vi mua hàng.