https://drive.google.com/file/d/1kDsZXF6eOfzqREUku5bgNgfw15fPYJZ1/view?usp=sharing
Phụ lục 3. Hồi quy bằng phần mềm R và Eviews
Chúng ta xét mô hình hồi quy hai biến số: \(𝑦 = 𝛽_1 + 𝛽_2𝑥 + 𝑢\)
Trong phần này, chúng ta lấy bộ số liệu CPS1988.csv làm ví dụ minh họa cho nội dung trên.
Trước hết chúng ta đánh giá mối quan hệ giữa biến education và biến wage.
Từ bộ số liệu ta nhận thấy khi education tăng lên thì wage cũng tăng, ta có đồ thị phân tán của wage theo education bằng câu lệnh sau:
plot(wage~education, pch = 16): đồ thi phân tán của của wage theo education.
CPS <- read.csv("CPS1988.csv", header = TRUE) # Gọi biến wage trong CPS1988 và đặt tên wage wage <- CPS1988\(wage # Gọi biến education trong CPS1988 và đặt tên education education <- CPs1988\)education plot(wage~education, pch = 16) plot
Dựa vào đồ thị phân tán ta cũng có kết quả như nhận xét ban đầu.
Tuy nhiên để đánh giá chính xác mối quan hệ tuyến tính giữa wage và education ta dùng hệ số tương quan để đánh giá.
Hệ số tương quan giữa wage và education được tính bằng hàm câu lệnh sau:
cor(wage, education): câu lệnh tính hệ số tương quan giữa wage và education.
cor(wage, education)
Kết quả câu lệnh cho thấy hệ số tương quan giữa wage và education là 0,301644.
Chúng ta có thể kiểm định giả thuyết tương quan giữa wage và education (tức là hệ số tương quan giữa wage và education bằng 0) với hàm cor.test () có sẵn trong R:
cor.test(wage, education)
Ta có p-value < \(2,2.10^(-16)\)+ < 0,05 , ta bác bỏ giả thuyết hệ số tương quan giữa wage và education bằng 0.
Để ước lượng mô hình hồi quy đơn bằng phần mềm R, chúng ta thực hiện câu lệnh sau:
reg <- lm(y~x, data = CPS1988): câu lệnh hồi quy tuyến tính 𝑦 theo 𝑥 và đặt tên reg; data = CPS1988: dữ liệu chính là đối tượng CPS1988.
Hàm lm (viết tắt từ linear model) trong R có thể tính toán các giá trị của \(\beta_1\) và \(\beta_2\) , cũng như \(\var\) một cách nhanh gọn. Chúng ta tiếp tục với ví dụ trên bằng R như sau:
reg <- lm(wage~education, data = CPS1988) reg
Trong lệnh lm(wage~education) cho ta thấy wage là một hàm của education. Kết quả tính toán của lm cho thấy \(𝛽_1\) = −12.83 và \(𝛽_2\) = 47.18. Nói cách khác với hai tham số này chúng ta có thể tính wage cho bất kỳ education nào trong mẫu dựa vào phương trình tuyến tính:
\(wage\) = −12,83 + 47,18.educaton
Thật ra, hàm lm còn cung cấp cho chúng ta nhiều thông tin khác và các thông tin đó được lưu trong đối tượng reg. Để xem các thông tin đó chúng ta thực hiện:
summary(reg)
Lệnh summary(reg) yêu cầu R liệt kê các tính toán trong reg gồm các phần sau:
Phần một mô tả chi tiết phần dư (Residuals) của hàm hồi quy
Phần hai trình bày ước số của \(𝛽_1\) và \(𝛽_2\) cùng với sai số chuẩn và giá trị của kiểm định t.
Giá trị kiểm định t cho \(𝛽_1\) là 53,085 với trị số 𝑝 < \(2 × 10^-12\), cho thấy \(𝛽_2\) không phải bằng 0.
Cột Coefficients: tên biến
Cột Estimate: giá trị của các hệ số hồi quy \(𝛽_i\)
Cột Std. Error: độ lệch chuẩn của các hệ số hồi quy \(𝛽_i\)
Cột t value (𝑡a = x2e)(x2): giá trị thống kê t, tiêu chuẩn kiểm định của bài toán kiểm định giả thuyết: 𝛽a = 0
Cột Pr(>|t|): p_value của bài toán kiểm định giả thuyết: \(𝛽_i\) = 0
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
Dấu hiệu cho biết các biến độc lập có nghĩa hay không, với các mức ý nghĩa cho trước: ‘*’: 1 ‰; ‘’: 1%; “0.05 “: 5%; “0.1”: 1%.
Ngoài ra, phần ba còn cho chúng ta một thông tin quan trọng, đó là trị số \(R^2\) hay hệ số xác định bội (coefficient of determination). Trị số R2 trong ví dụ này là 0.09099. Một hệ số cũng cần đề cập ở đây là hệ số điều chỉnh xác định bội (mà trong kết quả trên R gọi là “Adjusted R-squared”).
Với lệnh fitted() chúng ta có thể tính toán \(𝑦_i\) cho từng cá nhân như sau:
fitted(reg)
Chúng ta xét mô hình hồi quy k biến số:
\(𝑦 = 𝛽_1 + 𝛽_2𝑥_2 + ⋯ + 𝛽_i𝑥_i + ⋯𝛽_k𝑥_k + u\)
Để thực hiện hồi quy mô hình này ta thực hiện câu lệnh sau:
hoiquy <- lm(y~x2 + x3 + xi + ... + xk): câu lệnh hồi quy tuyến tính y theo các biến x2, x3, ..., xk.
Chúng ta tiếp tục sử dụng bộ số liệu CPS1988.csv để minh họa cho hồi quy đa biến. Ta thực hiện hồi quy wage theo education và experience ta được kết quả như sau:
reg <- lm(wage~education + experience, data = CPS1988) summary(reg)
Kết quả hàm hồi quy mẫu của wage theo education và experience:
\(wage\) = −385,0834 + 60,8964. education + 10,6057.experience
Giả sử chúng ta cần ước lượng mô hình sau:
𝑙𝑜𝑔 �1 −𝜋(𝜋𝑥()𝑥)Æ = 𝛼$ + 𝛽\(𝑥\) + 𝛽%𝑥% + ⋯ + 𝛽!𝑥!
Trong đó E(H) \(#E(H) được gọi là khả năng thành công và hàm 𝑙𝑜𝑔 ç\)#EE(H()H)é được gọi là hàm logit của 𝜋(𝑥) .
Việc ước lượng các hệ số của mô hình logit được thực hiện bằng phương pháp ước lượng hợp lí cực đại (maximum likelihood). Với R chúng ta có thể thực hiện ước lượng này lệnh glm (viết tắt của generalized linear models – các mô hình tuyến tính tổng quát) như sau:
logit <-glm(y~x1 + x2+ ...+ xk, family = binomial(link = “logit”), data = object): câu lệnh hồi quy logit của hàm y theo x1, x2, ..., xk và đặt tên logit.
Chúng ta lấy dữ liệu từ tập tin crab.csv minh họa cho hồi quy logit.
Trước tiên chúng ta mô tả các biến trong file crab.
# Đọc file crab.csv và đặt tên crab
t<-file.choose()
crab <- read.csv(t, header = T)Trong bộ dữ liệu gồm các biến sau:
sat: số vệ tinh xung quanh cua móng ngựa, nhận các giá trị từ 0 đến 15
𝑦: biến nhị phân nhận giá trị 1 nếu 𝑠𝑎𝑡 >= 1, nhận giá trị 0 nếu 𝑠𝑎𝑡 = 0
weight: trọng lượng cua móng ngựa
width: độ rộng mai cua móng ngựa
color: màu sắc cua móng ngựa
spine: số gai của cua móng ngựa
Để hồi quy mô hình Logit của biến 𝑦 theo các biến weight, width, color, spine bằng R ta thực hiện các câu lệnh sau:
crab <- read.csv("crab.csv" , header = TRUE): câu lệnh đọc file crab.csv và đặt tên crab.
logit <-glm(y~x1 + x2+ ...+ xk, family = binomial(link = “logit”), data = object): câu lệnh hồi quy mô hình Logit của biến y theo weight, width, color, spine và đặt tên logit; data = crab: dữ liệu chính là đối tượng crab.
summary(logit): câu lệnh hiển thị chi tiết kết quả hồi quy mô hình logit.
logit <-glm(y~weight+width+color+spine, family = binomial(link = "logit"), data = crab)
summary(logit)##
## Call:
## glm(formula = y ~ weight + width + color + spine, family = binomial(link = "logit"),
## data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1928 -0.9465 0.4919 0.8682 1.9892
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.0078436 3.8034845 -1.842 0.0654 .
## weight 0.0007949 0.0006917 1.149 0.2505
## width 0.2732525 0.1893304 1.443 0.1489
## color -0.5915286 0.2417241 -2.447 0.0144 *
## spine 0.2716837 0.2410407 1.127 0.2597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 225.76 on 172 degrees of freedom
## Residual deviance: 186.66 on 168 degrees of freedom
## AIC: 196.66
##
## Number of Fisher Scoring iterations: 4Từ bảng kết quả ta có mô hình hồi quy:
\(𝑙𝑜𝑔𝑖𝑡[𝜋°(𝑥)] = −7.5994 + 0,2733 . 𝑤𝑖𝑑𝑡ℎ + 0,7949. 𝑤𝑒𝑖𝑔ℎ𝑡 − 0,5915. 𝑐𝑜𝑙𝑜𝑟 + 0,2717. 𝑠𝑝𝑖𝑛𝑒\)
Chú ý: Đối với bài toán này, trên Eviews thao tác như sau:
B1. Nhập số liệu cho các biến Y, X1,…, Xm, trong đó biến đáp ứng Y là nhị phân (dữ liệu được gán 0 hoặc 1).
B2. H.quy ước lượng cho (1): → 𝑄𝑢𝑖𝑐𝑘 → 𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛 → (𝑔õ: 𝑦 𝑐 𝑥1 … 𝑥𝑚),𝑐ℎọ𝑛: 𝐵𝑖𝑛𝑎𝑟𝑦 → 𝑙𝑜𝑔𝑖𝑡(ℎ𝑜ặ𝑐 𝑝𝑟𝑜𝑏𝑖𝑡) → 𝑂𝐾
Trong ví dụ trên, Eviews cho kết quả:
Chúng ta có thể diễn giải kết quả trên như sau. Ta có hệ số của biến 𝑤𝑖𝑑𝑡ℎ là 0,2733, có nghĩa là độ rộng tăng 1 đơn vị thì tỷ lệ cược (Odds) được nhân lên \(𝑒^0,2733\) = 1,31423. Tương tự cho các hệ số của các biến còn lại, Tương tự, chúng ta có thể tính Odds cho các biến còn lại bằng câu lệnh sau:
exp(logit$coefficients): câu lệnh tính e mũ cho các hệ số hồi quy mẫu trong mô hình.
exp(logit$coef)## (Intercept) weight width color spine
## 0.0009047575 1.0007952608 1.3142320426 0.5534806048 1.3121719135Kết quả cho chúng ta các giá trị của \(𝑒Zf\) và \(𝑒x2\); chẳng hạn \(𝑒x% = 2,2143187226\).
Ngoài ra, dựa trên mô hình hồi quy logit ở trên chúng ta có thể dự báo xác suất có vệ tinh của con cua cái với 𝑤𝑖𝑑𝑡ℎ = 25, 𝑤𝑒𝑖𝑔ℎ𝑡 = 3, 𝑐𝑜𝑙𝑜𝑟 = 5 𝑣à 𝑠𝑝𝑖𝑛𝑒 = 3 như sau bằng câu lệnh sau:
dubao <- 1/(1+(exp(-logit\(coef[1]-logit\)coef[2]25-logit\(coef[3]*3- logit\)coef[4]5-logit$coef[5]*3))): câu lệnh thực hiện dự báo và đặt tên dubao
dubao: câu lệnh xem kết quả dự báo.
Kết quả sau khi thự hiện câu lệnh
1/(1+(exp(-logit$coef[1]-logit$coef[2]*25-logit$coef[3]*3- logit$coef[4]*5-logit$coef[5]*3)))## (Intercept)
## 0.0002457858Vậy xác suất con cua cái có vệ tinh với 𝑤𝑖𝑑𝑡ℎ = 25, 𝑤𝑒𝑖𝑔ℎ𝑡 = 3, 𝑐𝑜𝑙𝑜𝑟 = 5 𝑣à 𝑠𝑝𝑖𝑛𝑒 = 3 là 0, 371493.
Thống kê Wald: \(𝑍% = �e)x MxNÆ%\) (có phân phối xấp xỉ Chi – bình phương khi mẫu lớn). Kiểm định này gọi là kiểm định Wald.
Để thực hiện kiểm định Wald bằng phần mềm R ta làm như sau: B1. Cài đặt và gọi các packages sau:
library(mdscore)## Warning: package 'mdscore' was built under R version 3.6.3## Loading required package: MASSB2. Hồi quy mô hình logistic.
B3. Thực hiện kiểm định Wald bằng câu lệnh: wald.test(model , 𝑡𝑒𝑟𝑚); trong đó model (tên mô hình hồi quy), terms (số tham số trong mô hình)
Chúng ta sử dụng tập dữ liệu crab.csv tiến hành kiểm định Wald bằng R, ta được kết quả như sau:
# Đọc tập tin crab.csv và đặt ten crab
crab <- read.csv("crab.csv", header = TRUE)
# Hồi quy y theo width
logit <- glm(y~width, family = binomial(link="logit"), data = crab)
wald.test(logit,terms = 2)## $W
## [1] 23.88748
##
## $pvalue
## [1] 1.02134e-06
##
## attr(,"class")
## [1] "wald.test"Kết quả kiểm định Wald, giá trị thống kê \(𝑍^2\)= 23,88748, 𝑝 − 𝑣𝑎𝑙𝑢𝑒 = 1,02134 × \(10^-6\) < 0,05. Kết luận: với mức ý nghĩa 5%, sự xuất hiện của vệ tinh thực sự phụ thuộc vào chiều rộng của mai của con cua cái.
Đối với mô hình hồi quy logistic, để kiểm định giả thuyết \(𝐻_0\), người ta thường dùng thống kê tỷ số hợp lý (LR: 𝑃𝑠𝑒𝑢𝑑𝑜 − 𝑅%)
\[ 𝐿𝑅 = −2(𝐿_0 − 𝐿) \]
Dựa trên bộ dữ liệu crab.csv ta hồi quy logitics 𝑦 theo 𝑤𝑖𝑑𝑡ℎ, bài toán kiểm định giả thuyết \(𝐻_0:\) “Mô hình không phù hợp với dữ liệu điều tra” được thực hiện trên R như sau:
crab <- read.csv("crab.csv", header = TRUE)
logit <- glm(y~width, family = binomial(link="logit"), data = crab)Ta có kết quả: Thống kê 𝐿𝑅 = 31,306 và \(𝑝_value = 2,204 × 10\) < 0,05 : Bác bỏ \(𝐻_0\)
Kết luận, mô hình phù hợp với mức ý nghĩa 5%.
Để tìm chỉ số \(𝑃𝑠𝑒𝑢𝑑𝑜 − 𝑅^2\) bằng R ta thực hiện hồi quy logitics 𝑦 theo 𝑤𝑖𝑑𝑡ℎ từ tập dữ liệu crab.csv sau đó dùng lệnh pR2() ta được chỉ số cần tìm.
crab <- read.csv("crab.csv", header = TRUE)
fit <- glm(y~width, family = binomial(link="logit"), data = crab)Từ bảng kết quả ta có chỉ số Chỉ số \(Pseudo – R^2\) = 0,1386697
\[𝐵 = "∑" a9(𝑌a − 𝜋°a)% = ". 𝐷𝑒𝑣𝑖𝑎𝑛𝑐𝑒\]
Dựa trên mô hình hồi quy logitics 𝑦 theo 𝑤𝑖𝑑𝑡ℎ từ tập dữ liệu crab.csv, ta tìm hệ số Brier như sau:
crab <- read.csv("crab.csv", header = TRUE)
fit <- glm(y~width, family = binomial(link="logit"), data = crab)
pred.prob <- fitted(fit)
brierScore <- mean((pred.prob -crab$y)^2)
summary(brierScore)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1928 0.1928 0.1928 0.1928 0.1928 0.1928Từ bảng kết quả ta được hệ số Brier = 0,1928
Chú ý: Eviews đã cung cấp trên bảng hồi quy các chỉ số: Pseudo – R2, Deviance
Sử dụng gói “caret” chúng ta có thể nhanh chóng có được các thông tin về khả năng dự báo của mô hình bằng hàm confusionMatrix():
Để tìm ma trận nhầm lẫn trên R trước hết ta cài đặt và gọi packages “caret” và các packages liên quan.
library(survival)## Warning: package 'survival' was built under R version 3.6.3library(sandwich)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.6.3library(lattice)
library(zoo)## Warning: package 'zoo' was built under R version 3.6.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(caret)## Warning: package 'caret' was built under R version 3.6.3##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clusterSau khi cài đặt và gọi các packages trên ta tiếp tục các lệnh sau:
B1: Hồi quy mô hình logit.
B2: Dự báo xác suất cho 𝑦 dựa trên mô hình logit và làm tròn giá trị dự báo.
B3: Tìm ma trận nhầm lẫn bằng câu lệnh confusionMatrix( )
Chúng ta lấy dữ liệu từ tập tin crab.csv, để tìm ma trận nhầm lẫn theo các bước trên ta được kết quả như sau:
crab <- read.csv("crab.csv", header = TRUE)
logit <- glm(y~width, family = binomial(link="logit"), data = crab)
summary(logit)##
## Call:
## glm(formula = y ~ width, family = binomial(link = "logit"), data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0281 -1.0458 0.5480 0.9066 1.6942
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
## width 0.4972 0.1017 4.887 1.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 225.76 on 172 degrees of freedom
## Residual deviance: 194.45 on 171 degrees of freedom
## AIC: 198.45
##
## Number of Fisher Scoring iterations: 4pred <- fitted(logit)
yhat <- round(pred, 0)Từ bảng kết quả ta có:
Độ đặc hiệu: 0,7308; Độ nhạy: 0,6279; Độ chính xác toàn thể: 0,7052.
Chú ý: Để tìm ma trận nhầm lẫn trên Eviews thao tác như sau:
Sau khi chạy hồi quy:
Tìm cột dự báo \(𝑦° = @𝑟𝑜𝑢𝑛𝑑(𝑦 − 𝑟𝑒𝑠𝑖𝑑)\) cho Y:→ 𝑄𝑢𝑖𝑐𝑘 → 𝑆𝑒𝑟𝑖𝑒𝑠 𝑆𝑡𝑎𝑡𝑖𝑠𝑡𝑖𝑐𝑠 → 𝐻𝑖𝑠𝑡𝑜𝑔𝑟𝑎𝑚 𝑎𝑛𝑑 𝑆𝑡𝑎𝑡𝑠. Trong ô Series name: gõ @round(y – resid) → 𝑂𝐾, trong cửa sổ xuất hiện: → 𝑉𝑖𝑒𝑤 → 𝑆𝑝𝑟𝑒𝑎𝑑𝑆ℎ𝑒𝑒𝑑. Trong cửa sổ kết quả: chọn 𝑁𝑎𝑚𝑒 để lưu kết quả, chẳng hạn: 𝑦ℎ𝑎𝑡
Tiếp tục: Từ Workfile →view →Show: gõ: y yhat →OK. Xuất hiện cửa số Group: →view →N-way Tabulation →OK
Chúng ta cần ước lượng mô hình có dạng sau: \[𝑃𝑟𝑜𝑏𝑖𝑡(𝜋(𝑥)) = 𝛼 + 𝛽$𝑥$ + 𝛽%𝑥% + ⋯ + 𝛽!𝑥!\]
Để hồi quy mô hình Probit bằng R sử dụng câu lệnh sau:
probit<-glm(y~x1 + x2 + ... + xk, family = binomial(link = probit), data = file): câu lệnh hồi quy mô hình probit của y theo x1, x2, ..., xk và đặt tên probit.
Chúng ta lấy dữ liệu từ file crab.csv minh họa cho hồi quy mô hình Probit. Ta thực hiện các câu lệnh sau:
crab <- read.csv("crab.csv" , header = TRUE): câu lệnh đọc file crab.csv và đặt tên crab.
probit<-glm(y~width + color, family = binomial(link = “probit”), data = crab): câu lệnh hồi quy mô hình Probit của biến y theo width và color; data = crab: dữ liệu chính là đối tượng crab.
summary(probit): câu lệnh hiển thị kết quả hồi quy mô hình Probit
Chúng ta có kết quả sau khi thực hiện câu lệnh.
# Mô hình probit
crab<-read.csv(t, header = TRUE)
probit<-glm(y~width + color, family = binomial(link = "probit"), data = crab)
summary(probit)##
## Call:
## glm(formula = y ~ width + color, family = binomial(link = "probit"),
## data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1880 -0.9993 0.5418 0.8749 1.9665
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.80617 1.68232 -3.451 0.000558 ***
## width 0.27581 0.05939 4.644 3.41e-06 ***
## color -0.29188 0.13354 -2.186 0.028843 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 225.76 on 172 degrees of freedom
## Residual deviance: 189.07 on 170 degrees of freedom
## AIC: 195.07
##
## Number of Fisher Scoring iterations: 5Từ bảng kết quả, ta có mô hình Probit ước lượng như sau: \[𝑃𝑟𝑜𝑏𝑖𝑡𝜋°(𝑤𝑖𝑑𝑡ℎ, 𝑐𝑜𝑙𝑜𝑟) = −6,098048 + 0,275813. 𝑤𝑖𝑑𝑡ℎ − 0,291876. 𝑐𝑜𝑙𝑜𝑟\]
Suy ra:
\[𝜋°(𝑥, 𝐶𝑂𝐿𝑂𝑅) = Φ(−6,098048 + 0,275813. 𝑤𝑖𝑑𝑡ℎ − 0,291876. 𝑐𝑜𝑙𝑜𝑟)\]
Theo đó xác suất để một con cua cái có màu trung bình (𝑐𝑜𝑙𝑜𝑟 = 2), với độ rộng của mai là 𝑤𝑖𝑑𝑡ℎ = 25, ước tính là:
\[𝜋°(25; 2) = Φ(−6,098048 + 0,275813 × 25 − 0.291876 × 2) = 0,5832\]
Ký hiệu 𝜇(𝑥) là giá trị kỳ vọng cho biến ngẫu nhiên Poisson 𝑌 ứng với mức 𝑥 của biến 𝑋, tức là: 𝜇(𝑥) = 𝐸(𝑌|𝑋 = 𝑥). Mô hình loglinear Poisson có dạng:
\[𝑙𝑜𝑔[𝜇(𝑥)] = 𝛼 + 𝛽𝑥\]
Để hồi quy mô hình 𝑙𝑜𝑔[𝜇(𝑥)] = 𝛼 + 𝛽𝑥 bằng R, ta làm như sau:
poisson <- glm(μ(x) ~ x2 + x3 + + xk, family = poisson(link = “log”), data = object): câu lệnh hồi quy mô hình Poisson của biến y theo các biến x1, x2, ..., xk và đặt tên poisson.
Chúng ta lấy dữ liệu từ tập tin crab.csv minh họa cho hồi quy poisson, chúng ta thực hiện các câu lệnh sau:
crab <- read.csv("crab.csv" , header = TRUE): câu lệnh đọc file crab.csv và đặt tên crab.
poisson <- glm(sat ~width , family = poisson(link = log), data = crab): câu lệnh hồi quy mô hình Probit của biến sat theo width; data = crab: dữ liệu chính là đối tượng crab.
summary(poisson): câu lệnh hiển thị chi tiết kết quả hồi quy Poisson.
poisson <- glm(satell ~width , family = poisson(link = log), data = crab)
summary(poisson)##
## Call:
## glm(formula = satell ~ width, family = poisson(link = log), data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8526 -1.9884 -0.4933 1.0970 4.9221
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.30476 0.54224 -6.095 1.1e-09 ***
## width 0.16405 0.01997 8.216 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 632.79 on 172 degrees of freedom
## Residual deviance: 567.88 on 171 degrees of freedom
## AIC: 927.18
##
## Number of Fisher Scoring iterations: 6Chú ý: Để chạy mô hình Poisson trên Eviews, thao tác như sau:
B1. Nhập số liệu cho các biến Y, X, trong đó biến đáp ứng Y là Poisson (dữ liệu là các số tự nhiên).
B2. Hồi quy Poisson ước lượng: → 𝑄𝑢𝑖𝑐𝑘 → 𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛 → (𝑔õ: 𝑦 𝑐 𝑥), trong Method chọn Generalized Linear Model, trong Family: chọn Poisson, trong Link function: chọn log → 𝑂𝐾
Từ bảng kết quả, ta có mô hình ước lượng:
\[𝑙𝑜𝑔[𝜇̂(𝑥)] = 𝛼° + 𝛽{𝑥 = −3,304757 + 0,164045. 𝑥\]
Khi biến phản ứng 𝑌 có chỉ số (như kích thước quần thể) bằng t, thì dữ liệu tỷ lệ là các giá trị quan sát của 𝑌/𝑡. Giá trị kỳ vọng có điều kiện của tỷ lệ này là 𝜇(𝑥)/𝑡. Mô hình log - lin cho dữ liệu tỷ lệ có dạng:
\[𝑙𝑜𝑔[𝜇(𝑥)/𝑡] = 𝛼 + 𝛽𝑥\]
Chúng ta tiếp tục xem xét lại các dữ liệu cua móng ngựa. Mẫu này với 173 con cua cái có 66 giá trị chiều rộng 𝑋$ khác nhau. Ký hiệu là tổng số vệ tinh của t con cua cái có cùng một chiều rộng cụ thể, và cho 𝜇 = 𝐸(𝑌|𝑋$). Khi đó 𝜇/𝑡 là tỷ lệ của số vệ tinh kỳ vọng trên mỗi cua ở chiều rộng đó.
Nếu phân độ rộng 𝑋$ làm 8 nhóm, ta có bảng số liệu phân nhóm:
Nhập dữ liệu cho ba biến chiều rộng (X1) (nhập điểm giữa mỗi khoảng), số trường hợp (t) và số vệ tinh (Y) vào R và hồi quy Poisson cho tỷ lệ ta làm như sau:
X1 <- c(22.75, 23.75, 24.75, 25.75, 26.75, 27.75, 28.75, 29.75): câu lệnh tạo ta đối tượng X1.
Y <- c(14, 20, 67, 105, 63, 93, 71, 72): câu lệnh tạo ra đối tượng Y.
t <- c(14, 14, 28, 39, 22, 24, 18, 14): câu lệnh tạo ra đối tượng t. poisson <- glm(Y/t ~ X1, family = poisson(link="log")): câu lệnh hồi quy Poisson của Y/t theo X1
summary(poisson): câu lệnh hiển thị chi tiết kết quả hồi quy.
Thực hiện các câu lệnh trên ta có kết quả:
X1 <- c(22.75, 23.75, 24.75, 25.75, 26.75, 27.75, 28.75, 29.75)
Y <- c(14, 20, 67, 105, 63, 93, 71, 72)
t <- c(14, 14, 28, 39, 22, 24, 18, 14)
poisson <- glm(Y/t ~ X1, family = poisson(link="log"))## Warning in dpois(y, mu, log = TRUE): non-integer x = 1.428571## Warning in dpois(y, mu, log = TRUE): non-integer x = 2.392857## Warning in dpois(y, mu, log = TRUE): non-integer x = 2.692308## Warning in dpois(y, mu, log = TRUE): non-integer x = 2.863636## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.875000## Warning in dpois(y, mu, log = TRUE): non-integer x = 3.944444## Warning in dpois(y, mu, log = TRUE): non-integer x = 5.142857summary(poisson)##
## Call:
## glm(formula = Y/t ~ X1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.29899 -0.15919 -0.03878 0.17583 0.29244
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.17973 2.62234 -1.594 0.1110
## X1 0.19624 0.09597 2.045 0.0409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.74059 on 7 degrees of freedom
## Residual deviance: 0.29889 on 6 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4Chạy hồi quy cho mô hình tỷ lệ : 𝑙𝑜𝑔[𝜇(𝑥)/𝑡] = 𝛼 + 𝛽𝑥, nhận được kết quả hồi quy. Từ đó nhận được mô hình ước lượng:
\[𝑙𝑜𝑔 ø𝜇̂(𝑡𝑥)¿ = −4,179732 + 0,196243. 𝑥\]
Chú ý: Để chạy hồi quy Poisson tỷ lệ trên Eviews, thao tác như sau:
B1. Nhập số liệu các biến x, t, Y
B2. Hồi quy ước lượng cho (4): → 𝑄𝑢𝑖𝑐𝑘 → 𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛 → (𝑔õ: 𝑦/ 𝑡 𝑐 𝑥), trong ô Method chọn Generalized Linear Model, trong Family: chọn Poisson, trong Link function: chọn log → 𝑂𝐾
Ước lượng cho lượng vệ tinh bình quân cho một con cua cái có mức chiều rộng 𝑥 là:
\[𝜇̂(𝑥) = 𝑡. 𝑒𝑥𝑝{−4,179732 + 0,196243. 𝑥}\]
Với khoảng chiều rộng đầu tiên: 𝑥$ = 22,75, 𝑡$ = 14, ước lượng bằng R được thực hiện như sau:
14*exp(poisson\(coef[1]+poisson\)coef[2]22.75): câu lệnh dự báo số vẹ tinh bình quân cho con cua cái với 𝑥$ = 22,75, 𝑡$ = 14.
# Dự báo cho giá trị X1 = 22.75, t = 14
14*exp(poisson$coef[1]+poisson$coef[2]*22.75)## (Intercept)
## 18.61294Vậy ước lượng số lượng vệ tinh bình quân cho một con cua cái có mức chiều rộng 𝑥1 = 22,75 là: 18,61294
Trong mục này chúng ta đề cập đến mô hình hồi quy logistic mà trong đó biến đáp ứng Y có k + 1 giá trị khác nhau: 0, 1, 2, …, k (hay có các thuộc tính được gán số 0, 1, 2,…, k) với k > 1, và 𝑋 = (𝑋$, 𝑋%, … , 𝑋&) là tập biến giải thích.
Để xây dựng mô hình logistic đa thức, trước hết người ta chọn một thuộc tính nào đó của đáp ứng Y làm thuộc tính cơ sở, chẳng hạn thuộc tính Y = 0. Đặt 𝜋O(𝑥) = 𝑃(𝑌 = 𝑗|𝑋 = 𝑥), 𝑥 = (𝑥$, 𝑥%, … , 𝑥&), mô hình hồi quy logistic đa thức có dạng: \[𝐿𝑜𝑔 𝜋O(𝑥) 𝜋,(𝑥) = 𝛽O, + 𝛽O$. 𝑥$ + ⋯ + 𝛽O&. 𝑥&, 𝑗 = 1, 𝑘\]
Ví dụ. Sau khi tốt nghiệp phổ thông, học sinh Mĩ có ba lựa chọn để học tiếp. Một là vào các trường danh giá như Harvard (gọi là academic), hai là vào các trường nghề (vocation) và cuối cùng là trường đại cương (general). Cần nghiên cứu xem: điểm viết (write), tầng lớp xã hội (ses) ảnh hưởng như thế nào đến việc chọn trường (prog) của học sinh. Ta sử dụng bộ dữ liệu có tên mul.xlsx gồm 200 quan sát để phân tích. Đây là bộ dữ liệu lấy từ địa chỉ:
Để thực hiện hồi quy mô hình Logistic đa thức, trước hết cần cài đặt và gọi các gói sau:
install.packages("nnet")
install.packages("tidyverse")
library(foreign)
library(tidyverse)
library(nnet)
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
Tiếp theo, thực hiện các câu lệnh sau:
setwd(“D:/PhanmemR”): chỉ định R làm việc trên thư mục D:/PhanmemR.
library(xlsx): gọi packages xlsx TCC <- read.xlsx(“TCC.xlsx”, sheetIndex = 1, header = TRUE): đọc file có tên “TCC.xlsx”
mul <- read.xlsx("mul.xlsx", sheetIndex = 1, header = T): đọc tập tin
mul.xlsx và đặt tên mul; sheetIndex = 1: R chỉ đọc dữ liệu ở sheet 1.
Xem xu hướng chọn trường theo tầng lớp của họ thông qua câu lệnh.
with(mul, table(ses.group, prog.group)): câu lệnh tạo bảng ngẫu nhiên 2 chiều của biến ses và prog.
Kết quả câu lệnh
Tính điểm write trung bình cũng như độ lệch chuẩn của thí sinh ứng với từng nhóm trường mà họ thi vào bằng câu lệnh:
with(mul, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x))))): tính trung bình và độ lệch chuẩn của write ứng với từng nhóm thi vào.
Kết quả câu lệnh
mohinh <- multinom(prog.group ~ ses.group + write, data = mul): câu lệnh hồi quy mô hình logistic đa cấp độ. Kết qủa câu lệnh:
Xem xét hai biến định tính có bảng phân phối xác suất đồng thời là bảng 𝐼 × 𝐽 :
Trong đó 𝜋pO = 𝑃w𝐴 = 𝐴a, 𝐵 = 𝐵Ox; 𝜋aA = ∑O 𝜋pO ; 𝜋AO = ∑a 𝜋pO
Với mẫu ngẫu nhiên kích thước n, có bảng dữ liệu tương ứng:
Mô hình (X,Y): log 𝜇aO = 𝜆 + 𝜆a8 + 𝜆O=, ∀𝑖, 𝑗
Mô hình(X, Y, XY): log 𝜇aO = 𝜆 + 𝜆a8 + 𝜆O= + 𝜆aO 8=, ∀𝑖, 𝑗
Chạy mô hình (A, C, M): log 𝜇aO! = 𝜆 + 𝜆a: + 𝜆O¢ + 𝜆Ä! bằng các câu lệnh: loglin <-glm(count~A + C + M, family = poisson, data = Marijuana): câu lệnh hồi quy count (tần số) theo A (alcohol), C (cigarettes), M (marijuana) và đặt tên loglin.
Kết quả ước lượng mô hình thì R mặc định sẽ ước lượng cho tần số của ô đầu tiên (1,1) trong bảng, nghĩa là ứng với kết quả này thì:
\[𝑙𝑜𝑔(𝜇̂) = 4,17254 + 1,78511 + 0,64931 − 0,31542\]
Trước khi thực hiện việc ước lượng cho tần số kỳ vọng của các ô, R sẽ thực hiện các tính toán tương tự bảng trên cho các ô còn lại và sau đó thực hiện việc ước lượng tần số kỳ vọng cho các ô với câu lệnh sau:
uocluong <- fitted(loglin): ước lượng tần số kỳ vọng cho các ô và đặt tên uocluong.
Tương tự chúng ta tìm các mô hình hồi quy còn lại:
Mô hình (AC, AM, CM)
Mô hình (AC, M).
Mô hình (ACM).
Khi cỡ mẫu lớn, các thống kê likelihood-ratio 𝐺% và thống kê Pearson 𝜒%:
\[𝐺% = 2. s 𝑛aO! a,O,! . log ñ𝑛𝜇̂aO! aO!ó ; 𝜒% = s w𝑛aO!𝜇−̂aO!𝜇̂aO!x% a,O,!\]
Giả thuyết H0: Mô hình phù hợp.
Tiêu chuẩn bác bỏ giả thuyết H0 là:
\[𝜒% (ℎ𝑜ặ𝑐 𝐺%) ≥ 𝜒Z%(𝑑𝑓)\]
Hoặc đối với các thống kê này: 𝑃 − 𝑣𝑎𝑙𝑢𝑒 < 𝛼,
Để thực hiện kiểm định sự phù hợp của mô hình Logliniear bằng R, chúng ta thực hiện câu lệnh sau:
pchisq(deviance(model), df = df.residual(model), lower.tail = FALSE): trong đó deviance(model) là phần sai lệch của mô hình; df.residual(model) là bậc tự do của phần dư.
Ta lấy ví dụ minh họa cho kiểm định này như sau:
Marijuana <- read.table("Marijuana.txt", header = TRUE): câu lệnh đọc tập tin và đặt tên Marijuana.
loglin <- glm(count ~ (A + C + M)^2, family = poisson, data = Marijuana): câu lệnh hồi quy mô hình (AC, AM, CM) và đặt tên loglin.
pchisq(deviance(loglin), df = df.residual(loglin), lower.tail = F): câu lệnh kiểm định sự phù hợp mô hình.
Thực hiên các câu lệnh, ta được kết quả
Từ bảng kết quả trên ta có 𝑝 − 𝑣𝑎𝑢𝑙𝑒 = 0,5408396 > 0,05 ta chấp nhận H0. Vậy với mức ý nghĩa 5%, mô hình phù hợp.