library(tidyverse)
library(epitools)
library(DescTools)
library(DT)
library(energy)
options(digits = 4)

1 Nhắc Lại Một Số Khái Niệm

Định nghĩa về thống kê

  • Thống kê (theo nghĩa hẹp): Đếm!

  • Thống kê (theo nghĩa rộng): Là một hệ thống các phương pháp được dùng để: Thu thập, tóm tắt, trình bày, phân tích dữ liệu, dự báo dự báo.

Phân loại thống kê - Thống kê mô tả

  • Thống kê suy diễn

Tổng thể (population): Là đối tượng (gồm nhiều phần tử) mà chúng ta cần nghiên cứu.

Mẫu (sample): Là một tập hợp con của tổng thể.

Quan sát (observation): Là một phần tử của mẫu.

Biến (thống kê - variable): Là một đặc điểm hay 1 tính chất của một phần tử trong tổng thể.

Dữ liệu (data) là những thông tin mà chúng ta có thể ghi nhận lại. Dữ liệu được chia làm 2 loại (nếu dựa vào tính chất):

  • Dữ liệu định tính Là những thông tin được ghi nhận lại, nhưng những thông tin này không là những con số. Ví dụ: Màu mắt: Xanh, Nâu; Loại xe: Số sàn, số Tự động; Nghề nghiệp: Công nhân, Giáo viên…

  • Dữ liệu định lượng Là những thông tin được ghi nhận lại và những thông tin này là những con số. Ví dụ: Chiều cao; Chi tiêu, Tốc độ,…

Các loại thang đo

Thang đo là một tiêu chuẩn. Trong quá trình nghiên cứu chúng ta cần “đo” đối tượng (quan sát - observation) nào thì chúng ta sẽ lấy đối tượng đó với tiêu chuẩn.

Có 4 loại thang đo cơ bản.

  • Thang đo Định danh.
  • Thang đo Thứ bậc.
  • Thang đo Khoảng.
  • Thang đo Tỉ lệ.

Các đặc trưng đo lường phổ biến:

  • Trung bình
  • Trung vị
  • Mode
  • Phương sai
  • Độ lệch chuẩn
  • Khoảng biến thiên
  • Tứ phân vị
data(cars)
str(cars)
datatable(cars)
summary(cars)

2 Phân Phối Poisson - Poisson Distribution

2.1 Định nghĩa

Biến ngẫu nhiên \(X\) có phân phối Poisson là biến ngẫu nhiên dùng để mô tả cho số lần xảy ra của một sự việc/biến cố (event) mà chúng ta quan tâm xảy ra trong một khoảng thời gian hoặc không gian cho trước.

Xác suất để biến ngẫu nhiên này nhận một giá trị cụ thể được tính bằng công thức:

\[ P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!} \]

2.2 Một vài ứng dụng (ví dụ thực tế):

  • Số tin nhắn quảng cáo nhận được mỗi giờ.
  • Số lần truy cập vào website thư viện mỗi tuần.
  • Số người bị đột quỵ ở mỗi phường trong 1 tuần.

2.3 Mô phỏng phân phối Poisson - Poisson distribution simulation

ps <- rpois(105,4.5)

ps <- as.data.frame(table(ps))

ps <- ps |> mutate(prop = Freq/sum(Freq))

ps |>  ggplot(aes(ps,prop)) +
  geom_point(color = 'red', size = 3)

ps |>  ggplot(aes(ps,prop)) +
  geom_col()

p2 <- data.frame(poi = rpois(300,2))

p6 <- data.frame(poi = rpois(300,6))

p2$lambda <- '2'

p6$lambda <- '6'

p2p6 <- rbind(p2,p6)

p2p6 |>
  ggplot(aes(poi, fill = lambda)) +
  geom_bar()

tmp <- p2p6 |>
  group_by(lambda,poi) |>
  summarise(n = n())
## `summarise()` has grouped output by 'lambda'. You can override using the
## `.groups` argument.
tmp |> ggplot(aes(poi, n, color = lambda)) + geom_point()

2.4 Kiểm Định Về Phân Phối Poisson

Kiểm định một dãy số có phân phối Poisson không? Với giả thuyết \(H_0: X\) có phân phối Poisson. Ví dụ:

X <- rpois(n = 100, lambda = 2.6)

poisson.tests(X, R = 3)
##       estimate statistic p.value      method
## M-CvM     2.67    0.4170  0.3333  M-CvM test
## M-AD      2.67    2.2628  0.3333   M-AD test
## E         2.67    0.7501  0.6667 Energy test

3 Phân Phối Nhị Thức - Binomial distribution

3.1 Định nghĩa

Biến ngẫu nhiên \(X\) có phân phối nhị thức là biến ngẫu nhiên dùng để mô tả cho số lần thành công của một dãy những sự việc (biến cố) có những tính chất sau: + Dãy gồm \(n\) sự việc độc lập nhau. + Xác suất Thành công là bằng nhau và bằng \(p\) cho từng sự việc.

Khi đó xác suất để \(X\) nhận một giá trị cụ thể là: \[P(X=k) = C_n^kp^k(1-p)^{n-k}\]

3.2 Một vài ứng dụng (ví dụ thực tế):

  • Tung đồng xu 10 lần, xác suất để có 3 lần mặt ngửa xuất hiện.
  • Mua vé số 10 lần, xác suất trúng số 2 lần.
  • Số mail quảng cáo/rác trong 30 mail nhận được.
  • Số người được nhận vào làm việc khi có 30 người nộp hồ sơ.

3.3 Mô phỏng phân phối Nhị Thức - Binomial distribution simulation

bs <- rbinom(100, 15, .2)
bs <- as.data.frame(table(bs))
bs <- bs |> mutate(prop = Freq/sum(Freq))
ggplot(bs, aes(bs, prop)) + geom_point(size = 3, color = 'red')

ggplot(bs, aes(bs, prop)) + geom_col()

4 Thống kê mô tả cho dữ liệu định tính

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

4.1.1 Bảng tần số (frequency table) - bảng tần suất

Trong phần này chúng ta sẽ sử dụng dữ liệu về thông mua hàng của khách hàng tại một số siêu thị như sau:

d <- read.csv('Supermarket Transactions.csv',header = T)

str(d)
## '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 ...
datatable(d)
table(d$Gender)
## 
##    F    M 
## 7170 6889
table(d$Gender)/sum(table(d$Gender))
## 
##    F    M 
## 0.51 0.49
table(d$StateorProvince)
## 
##        BC        CA        DF  Guerrero   Jalisco        OR  Veracruz        WA 
##       809      2733       815       383        75      2262       464      4567 
##   Yucatan Zacatecas 
##       654      1297

4.1.2 Đồ Thị Cột

d |> ggplot(aes(Gender)) +
  geom_bar()

d |> ggplot(aes(Gender)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  ylab('Tỷ lệ %') + xlab('Giới tính')
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

d |> ggplot(aes(fct_infreq(StateorProvince))) +
  geom_bar(color = 'blue', fill = 'green')

d |> ggplot(aes(fct_infreq(StateorProvince))) +
  geom_bar(aes(y = (..count..)/sum(..count..)), color = 'blue', fill = 'green')

4.1.3 Đồ thị bánh

tmp <- d
tmp <- table(d$Gender)
tmp <- d |> group_by(Gender) |> summarise(freq = n()) |> mutate(tmp, per = freq/sum(freq))
tmp |> ggplot(aes(x = '', y = per, fill = Gender)) +
  geom_bar(stat = 'identity') + 
  coord_polar('y')

tmp <- d
tmp <- table(d$StateorProvince)
tmp <- d |> group_by(StateorProvince) |> summarise(freq = n()) |> mutate(tmp, per = freq/sum(freq))
tmp |> ggplot(aes(x = "", y = per, fill = StateorProvince)) + geom_col() + coord_polar('y')

4.2 Thống kê mô tả cho 2 biến

4.2.1 Đồ thị cột

ggplot(d, aes(Gender, fill = ProductFamily)) + geom_bar(position = 'dodge')

d |>group_by(Gender) |> summarise(M = mean(Revenue)) |> ggplot(aes(Gender,M)) + geom_col(fill = 'red')

d |> group_by(Gender) |> summarise(M = mean(UnitsSold, na.rm = T)) |> ggplot(aes(Gender, M)) + geom_col(fill = 'red')

4.2.2 Bảng tần số ( contingency tabble) - bảng tần suất

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

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

Chúng ta sẽ lập bảng ngẫu nhiên 2 chiều cho 2 biến Homeowner, ProductFamily như sau:

tmp <- table(d$Homeowner, d$ProductFamily)

tmp
##    
##     Drink Food Non-Consumable
##   N   504 4023           1088
##   Y   746 6130           1568

Hoặc chúng ta cũng có thể lập bảng ngẫu nhiên 2 chiều cho biến Homeoner và biến ProductFamily như sau:

tmp <- table(d$Homeowner, d$ProductFamily)

tmp
##    
##     Drink Food Non-Consumable
##   N   504 4023           1088
##   Y   746 6130           1568

Lưu ý: Biến nào làm tham số thứ nhất sẽ được xắp xếp theo hàng.

Nếu trong bảng ngẫu nhiên chúng ta tính tổng theo hàng và theo cột thì chúng ta gọi các tổng này là tần số biên (Marginal).

Chúng ta sẽ thực hiện như sau:

tmp <- table(d$Homeowner, d$ProductFamily)

addmargins(tmp)
##      
##       Drink  Food Non-Consumable   Sum
##   N     504  4023           1088  5615
##   Y     746  6130           1568  8444
##   Sum  1250 10153           2656 14059

Nếu chúng ta muốn tính bảng tần suất, tần suất biên tương ứng thì chúng ta thực hiện như sau:

tmp <- table(d$Homeowner, d$ProductFamily)

tmp <- prop.table(tmp)

tmp
##    
##       Drink    Food Non-Consumable
##   N 0.03585 0.28615        0.07739
##   Y 0.05306 0.43602        0.11153
addmargins(tmp)
##      
##         Drink    Food Non-Consumable     Sum
##   N   0.03585 0.28615        0.07739 0.39939
##   Y   0.05306 0.43602        0.11153 0.60061
##   Sum 0.08891 0.72217        0.18892 1.00000

Chúng ta cũng có thể mở rộng bảng ngẫu nhiên ra cho 3 biến, 4 biến và nhiều hơn nữa.

Ví dụ chúng ta có thể lập bảng ngẫu nhiên 3 chiều cho 3 biến: Gender, Homeowner, Country như sau:

tmp <- ftable(d$Gender, d$Homeowner, d$Country)

tmp
##      Canada Mexico  USA
##                        
## F N     136    866 1824
##   Y     237   1190 2917
## M N     184    607 1998
##   Y     252   1025 2823
addmargins(tmp)
##                     Sum
##     136  866 1824  2826
##     237 1190 2917  4344
##     184  607 1998  2789
##     252 1025 2823  4100
## Sum 809 3688 9562 14059

Hoặc chúng ta cũng có thể lập bảng tần suất tương ứng cho 3 biến Gender, Homeowner, Country như sau:

tmp <- ftable(d$Gender, d$Homeowner,d$Country)

tmp <- prop.table(tmp)

round(addmargins(tmp),4)
##                             Sum
##     0.0097 0.0616 0.1297 0.2010
##     0.0169 0.0846 0.2075 0.3090
##     0.0131 0.0432 0.1421 0.1984
##     0.0179 0.0729 0.2008 0.2916
## Sum 0.0575 0.2623 0.6801 1.0000

Hoặc chúng ta có thể lập các bảng ngẫu nhiên 2 chiều và bảng tần suất tương ứng cho 3 biến Gender, Homeowner, Country một cách khác như sau:

tmp <- table(d$Gender, d$Homeowner,d$Country)

tmp
## , ,  = Canada
## 
##    
##        N    Y
##   F  136  237
##   M  184  252
## 
## , ,  = Mexico
## 
##    
##        N    Y
##   F  866 1190
##   M  607 1025
## 
## , ,  = USA
## 
##    
##        N    Y
##   F 1824 2917
##   M 1998 2823
tmp <- prop.table(tmp)

tmp
## , ,  = Canada
## 
##    
##            N        Y
##   F 0.009674 0.016858
##   M 0.013088 0.017924
## 
## , ,  = Mexico
## 
##    
##            N        Y
##   F 0.061598 0.084643
##   M 0.043175 0.072907
## 
## , ,  = USA
## 
##    
##            N        Y
##   F 0.129739 0.207483
##   M 0.142115 0.200797
addmargins(tmp)
## , ,  = Canada
## 
##      
##              N        Y      Sum
##   F   0.009674 0.016858 0.026531
##   M   0.013088 0.017924 0.031012
##   Sum 0.022761 0.034782 0.057543
## 
## , ,  = Mexico
## 
##      
##              N        Y      Sum
##   F   0.061598 0.084643 0.146241
##   M   0.043175 0.072907 0.116082
##   Sum 0.104773 0.157550 0.262323
## 
## , ,  = USA
## 
##      
##              N        Y      Sum
##   F   0.129739 0.207483 0.337222
##   M   0.142115 0.200797 0.342912
##   Sum 0.271854 0.408279 0.680134
## 
## , ,  = Sum
## 
##      
##              N        Y      Sum
##   F   0.201010 0.308984 0.509994
##   M   0.198378 0.291628 0.490006
##   Sum 0.399388 0.600612 1.000000

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

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

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

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

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

4.4 Rủi ro tương đối (Relative Risk/Risk Ratio)

Ký hiệu \(\pi_i\) là tỷ lệ “thành công” của biến phụ thuộc (response variable) tương ứng với từng biểu hiện của biến độc lập.

Từ bảng tần xuất, chúng ta tính \(\frac{\pi_1}{\pi_2}\), phân số này gọi là Rủi ro tương đối (Relative risk) giữa 2 biểu hiện khác nhau của biến phụ thuộc.

Từ bảng

tmp <- table(d$MaritalStatus, d$Homeowner)

addmargins(tmp)
##      
##           N     Y   Sum
##   M    1719  5147  6866
##   S    3896  3297  7193
##   Sum  5615  8444 14059

Chúng ta tính rủi ro tương đối (relative risk).

## Lệnh này tính relativerisk cho vị trí 12 và 22
RelRisk(tmp) #package DescTools
## [1] 0.4622
## Lệnh này sẽ tính relativerisk cho vị trí 22 và 12
riskratio(tmp) #package epitool
## $data
##        
##            N    Y Total
##   M     1719 5147  6866
##   S     3896 3297  7193
##   Total 5615 8444 14059
## 
## $measure
##    risk ratio with 95% C.I.
##     estimate  lower  upper
##   M   1.0000     NA     NA
##   S   0.6114 0.5942 0.6292
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   M         NA           NA         NA
##   S          0   1.822e-277 3.663e-272
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
## Khi có thêm tham số rev = 'c' thì sẽ thực hiện việc
## đổi chỗ 2 cột trong bảngngẫu nhiên.
riskratio(tmp, rev = 'c') #package epitool
## $data
##        
##            Y    N Total
##   M     5147 1719  6866
##   S     3297 3896  7193
##   Total 8444 5615 14059
## 
## $measure
##    risk ratio with 95% C.I.
##     estimate lower upper
##   M    1.000    NA    NA
##   S    2.163 2.066 2.266
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   M         NA           NA         NA
##   S          0   1.822e-277 3.663e-272
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
## Câu lệnh này sẽ cung cấp nhiều thông tin hơn và thuận tiện hơn
## khi chúng ta biết thêm khái niệm oddratio trong phần tiếp theo
epitab(tmp, method = 'riskratio', rev = 'c')
## $tab
##    
##        Y     p0    N     p1 riskratio lower upper    p.value
##   M 5147 0.7496 1719 0.2504     1.000    NA    NA         NA
##   S 3297 0.4584 3896 0.5416     2.163 2.066 2.266 1.822e-277
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

4.5 Tỷ lệ chênh (Odd Ratio)

Nếu gọi xác suất “thành công” của biểu hiện thứ \(i\) của biến độc lập là \(\pi_i\) thì chúng ta kí hiệu Tỷ lệ cược (odd) của biểu hiện này là \(odd_i\) và được định nghĩa như sau: \[odd_i = \frac{\pi_i}{1-\pi_i}\] Nghĩa là chúng ta tính tỷ lệ thành công theo từng hàng trong bảng ngẫu nhiên.

Tỷ lệ chênh của biểu hiện thứ \(i\) và biểu hiện thứ \(j\) được kí hiệu là \(\theta_{ij}\) và được định nghĩa: \[\theta = \frac{odd_i}{odd_j} = \frac{\frac{\pi_i}{1-\pi_i}}{\frac{\pi_j}{1-\pi_j}} = \frac{\pi_i(1-\pi_j)}{\pi_j(1-\pi_i)}\] Ví dụ chúng sẽ tính odd cho bảng ngẫu nhiên sau với “thành công” chúng ta hiểu là đã có nhà.

tmp <- table(d$MaritalStatus, d$Homeowner)

tmp
##    
##        N    Y
##   M 1719 5147
##   S 3896 3297
## Câu lệnh này sẽ tính oddratio cho vị trí 11 và 21
OddsRatio(tmp) #package DescTools
## [1] 0.2826
## Câu lệnh này sẽ tính oddratio cho vị trí 22 và 12
oddsratio(tmp) #package epitools
## $data
##        
##            N    Y Total
##   M     1719 5147  6866
##   S     3896 3297  7193
##   Total 5615 8444 14059
## 
## $measure
##    odds ratio with 95% C.I.
##     estimate  lower  upper
##   M   1.0000     NA     NA
##   S   0.2827 0.2631 0.3036
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   M         NA           NA         NA
##   S          0   1.822e-277 3.663e-272
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
##Khi có thêm tham số rev = 'c' thì sẽ thực hiện việc
## đổi chỗ 2 cột trong bảngngẫu nhiên.
oddsratio(tmp, rev = 'c') #package epitools
## $data
##        
##            Y    N Total
##   M     5147 1719  6866
##   S     3297 3896  7193
##   Total 8444 5615 14059
## 
## $measure
##    odds ratio with 95% C.I.
##     estimate lower upper
##   M    1.000    NA    NA
##   S    3.538 3.294 3.801
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   M         NA           NA         NA
##   S          0   1.822e-277 3.663e-272
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
epitab(tmp, method = 'oddsratio') #package epitools
## $tab
##    
##        N     p0    Y     p1 oddsratio  lower  upper    p.value
##   M 1719 0.3061 5147 0.6095    1.0000     NA     NA         NA
##   S 3896 0.6939 3297 0.3905    0.2826 0.2631 0.3036 1.822e-277
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

5 Thống kê suy diễn cho dữ liệu định tính

5.1 Kiểm định tính độc lập cho 2 biến định tính

Giả thuyết \(H_0: X, Y\) độc lập.

Giá trị kiểm định:

Phương pháp Chi bình phương:

\[\chi^2 = { \underset{i,j}{\sum}}\frac{(n_{ij}-\hat{\mu}_{ij})^2}{\hat{\mu}_{ij}}\] Với \(n_{ij}\) là giá trị của ô \(i,j\), \(\hat{\mu}_{ij} = \frac{n_{i+}n_{+j}}{n}\)

Ví dụ 2: Chúng ta sẽ kiểm định sự độc lập giữa 2 biến MaritalStatus, Homeowner như sau:

tmp <- table(d$MaritalStatus, d$Homeowner)

tmp
##    
##        N    Y
##   M 1719 5147
##   S 3896 3297
chisq.test(tmp)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tmp
## X-squared = 1241, df = 1, p-value <2e-16

Qua kết quả kiểm định cho ta \(p-value < 2e-16 < 0.05\), nên bác bỏ \(H_0\), nghĩa là giữa việc sở hữu nhà và tình trạng hôn nhân là có liên quan với nhau.

Ví dụ 2: Chúng ta kiểm định sự độc lập giữa 2 biến MaritalStatus, ProductFamily

tmp <- table(d$MaritalStatus, d$ProductFamily)

tmp
##    
##     Drink Food Non-Consumable
##   M   628 4938           1300
##   S   622 5215           1356
chisq.test(tmp)
## 
##  Pearson's Chi-squared test
## 
## data:  tmp
## X-squared = 1.2, df = 2, p-value = 0.6

Kết quả \(p-value = 0.6\) nên chưa đủ chứng cứ/thông tin/cơ sở để nói rằng tình trạng hôn nhân độc lập với việc mua sắm.

5.2 Khoảng ước lượng cho tỉ lệ

Công thức ước lượng tỷ lệ (cho 1 tổng thể): \[\hat{p} - Z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le P \le \hat{p} - Z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\] Ước lượng tỷ lệ người mua nhiều hơn 4 đơn vị (UnitsSold) đồng thời kiểm định xem tỷ lệ (%) người mua nhiều hơn 4 sản phẩm có phải là 36% không (nghĩa là chúng ta kiểm định giả thuyết: \(H_0: p = 0.36\)).

tmp <- d[d$UnitsSold > 4,]

prop.test(length(tmp$UnitsSold), length(d$UnitsSold), p = 0.36)
## 
##  1-sample proportions test with continuity correction
## 
## data:  length(tmp$UnitsSold) out of length(d$UnitsSold), null probability 0.36
## X-squared = 3.7, df = 1, p-value = 0.05
## alternative hypothesis: true p is not equal to 0.36
## 95 percent confidence interval:
##  0.3598 0.3758
## sample estimates:
##      p 
## 0.3678

Ước lượng sự chênh lệch về tỷ lệ mua nhiều hơn 4 sản phẩm giữa nam và nữ.

\[ (\hat{p}_1-\hat{p}_2)- Z_{\alpha/2} \sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2}}\le p_1 - p_2 \le (\hat{p}_1-\hat{p}_2)+ Z_{\alpha/2} \sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2}}\]

Thực hiện bài toán kiểm định giả thuyết sự bằng nhau về tỷ lệ mua nhiều hơn 4 sản của 2 tổng thể (nam và nữ), nghĩa là chúng ta thực hiện bài toán kiểm định \(H_0: p_1 = p_2\). Công thức:

tmpm <- d[d$Gender == 'M',]
tmpf <- d[d$Gender == 'F',]

tmpm3 <- tmpm[tmpm$UnitsSold > 4,]
tmpf3 <- tmpf[tmpf$UnitsSold > 4,]

a <- c(nrow(tmpm), nrow(tmpf))
b <- c(nrow(tmpm3), nrow(tmpf3))

prop.test(b,a)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  b out of a
## X-squared = 0.75, df = 1, p-value = 0.4
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.008922  0.023254
## sample estimates:
## prop 1 prop 2 
## 0.3715 0.3643

5.3 Các mô hình hồi quy

Mô hình hồi quy được chia làm 2 nhóm:

  • Mô hình hồi quy cổ điển \[\mu_i \equiv E[y_i] = \beta_0 + \beta_1x_{1i}+\dots + \beta_kx_{ki}\]

Để mô hình hồi quy cổ điển phải thỏa mãn 3 điều kiện sau:

  1. Linearity of the association between predictors and outcome variable.
  2. Gaussian distribution of responses.
  3. constant variance of response distribution.
  • Mô hình hồi quy tổng quát: \[\mu_i = g^{-1}(\beta_0 + \beta_1x_{1i}+\dots + \beta_kx_{ki})\] Hoặc \[ g(\mu_i) = \beta_0 + \beta_1x_{1i}+\dots + \beta_kx_{ki}\] Trong mô hình này hàm \(g(.)\) được gọi là hàm liên kết (link function) và \(g(.)\) phải là hàm đơn điệu.

Trong một mô hình hồi quy tuyến tính tổng quát sẽ có 3 thành phần:

  • Thành phần hệ thống (systematic component): Là thành tổ hợp tuyến tính của các biến độc lập.
  • Thành phần ngẫu nhiên (random component): Là phân phối xác suất của biến phụ thuộc (dependent/response variable).
  • Hàm liên kết (link function).

5.3.1 Một số ví dụ về hồi quy cổ điển

Vẽ đồ thị dạng scatter cho dữ liệu.

data(mtcars)

mtcars |> ggplot(aes(mpg, hp)) +
  geom_point(color = 'red')

Thêm đồ thị của hàm hồi quy vào đồ thị scatter.

data(mtcars)

mtcars |> ggplot(aes(mpg, hp)) + 
  geom_point(color = 'red') + 
  geom_smooth(formula = y ~ x, method = 'lm')

data(mtcars)

mtcars |> ggplot(aes(mpg, hp)) + 
  geom_point(color = 'red') + 
  geom_smooth(formula = y ~ poly(x,3), method = 'lm')

Ước lượng hàm hồi quy tổng quát

fit <- lm(mpg ~ hp, data = mtcars)

summary(fit)
## 
## Call:
## lm(formula = mpg ~ hp, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.712 -2.112 -0.885  1.582  8.236 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.0989     1.6339   18.42  < 2e-16 ***
## hp           -0.0682     0.0101   -6.74  1.8e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.86 on 30 degrees of freedom
## Multiple R-squared:  0.602,  Adjusted R-squared:  0.589 
## F-statistic: 45.5 on 1 and 30 DF,  p-value: 1.79e-07
fit <- lm(mpg ~ hp + wt, data = mtcars)

summary(fit)
## 
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.941 -1.600 -0.182  1.050  5.854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879   23.28  < 2e-16 ***
## hp          -0.03177    0.00903   -3.52   0.0015 ** 
## wt          -3.87783    0.63273   -6.13  1.1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.59 on 29 degrees of freedom
## Multiple R-squared:  0.827,  Adjusted R-squared:  0.815 
## F-statistic: 69.2 on 2 and 29 DF,  p-value: 9.11e-12
fit <- lm(mpg ~ hp + wt + factor(am), data = mtcars)
summary(fit)
## 
## Call:
## lm(formula = mpg ~ hp + wt + factor(am), data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.422 -1.792 -0.379  1.225  5.532 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.00288    2.64266   12.87  2.8e-13 ***
## hp          -0.03748    0.00961   -3.90  0.00055 ***
## wt          -2.87858    0.90497   -3.18  0.00357 ** 
## factor(am)1  2.08371    1.37642    1.51  0.14127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.54 on 28 degrees of freedom
## Multiple R-squared:  0.84,   Adjusted R-squared:  0.823 
## F-statistic:   49 on 3 and 28 DF,  p-value: 2.91e-11
fit <- lm(am ~ mpg, data = mtcars)
summary(fit)
## 
## Call:
## lm(formula = am ~ mpg, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.620 -0.296 -0.105  0.252  0.847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.5915     0.2534   -2.33  0.02645 *  
## mpg           0.0497     0.0121    4.11  0.00029 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.406 on 30 degrees of freedom
## Multiple R-squared:  0.36,   Adjusted R-squared:  0.338 
## F-statistic: 16.9 on 1 and 30 DF,  p-value: 0.000285

Kết quả thu được có ý nghĩa không? Tại sao?

#fit <- lm(factor(am) ~ mpg, data = mtcars)
#summary(fit)

Câu lệnh báo lỗi! Lý do là gì?

6 Ước lượng mô hình hồi quy tổng quát

fit <- glm(mpg ~ hp,data = mtcars)

summary(fit)
## 
## Call:
## glm(formula = mpg ~ hp, data = mtcars)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.0989     1.6339   18.42  < 2e-16 ***
## hp           -0.0682     0.0101   -6.74  1.8e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 14.92)
## 
##     Null deviance: 1126.05  on 31  degrees of freedom
## Residual deviance:  447.67  on 30  degrees of freedom
## AIC: 181.2
## 
## Number of Fisher Scoring iterations: 2
fit <- glm(mpg ~ hp + wt + factor(am), data = mtcars)

summary(fit)
## 
## Call:
## glm(formula = mpg ~ hp + wt + factor(am), data = mtcars)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.00288    2.64266   12.87  2.8e-13 ***
## hp          -0.03748    0.00961   -3.90  0.00055 ***
## wt          -2.87858    0.90497   -3.18  0.00357 ** 
## factor(am)1  2.08371    1.37642    1.51  0.14127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.439)
## 
##     Null deviance: 1126.05  on 31  degrees of freedom
## Residual deviance:  180.29  on 28  degrees of freedom
## AIC: 156.1
## 
## Number of Fisher Scoring iterations: 2

6.1 Ước lượng hàm hồi quy cho dữ liệu nhị phân

Dữ liệu nhị phân là dữ liệu định tính chỉ nhận 2 giá trị: Đúng/Sai; Đồng ý/Không đồng ý; Có/Không; Thành công/Thất bại;….

Một số bài toán đặt ra:

  1. Chọn ngẫu nhiên 1 sinh viên, tính xác suất để chon được sinh viên giỏi.
  2. Chọn ngẫu nhiên 1 khách hàng, tính xác suất để chọn được người đã có gia đình.
  3. Chọn ngẫu nhiên 1 chiếc xe hơi, tính xác suất xe sử dụng hộp số tự động.

Mô hình logit \[logit(\pi) = log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1X_1 +\beta_2X_2+\dots+\beta_kX_k\] Với \(\pi\) là xác suất để biến phụ thuộc nhận giá trị “thành công”. Đối với hàm glm phạm trù thứ hai của biến nhị phân thể hiện cho “thành công” (chúng ta sử dụng hàm levels() để kiểm tra thứ tự của các phạm trù của một biến định tính.

fitted <- glm(factor(am) ~ hp, family = binomial(link = 'logit'), data = mtcars)

summary(fitted)
## 
## Call:
## glm(formula = factor(am) ~ hp, family = binomial(link = "logit"), 
##     data = mtcars)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.77661    0.91543    0.85     0.40
## hp          -0.00812    0.00607   -1.34     0.18
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43.230  on 31  degrees of freedom
## Residual deviance: 41.228  on 30  degrees of freedom
## AIC: 45.23
## 
## Number of Fisher Scoring iterations: 4
ggplot(data = fitted, aes(fitted$data$hp ,fitted$fitted.values)) + geom_point() + xlab('Biến hp') + ylab('Xác suất')

ggplot(data = fitted, aes(fitted$fitted.values)) +
  geom_histogram(binwidth = .001, color = 'green') +
  geom_density(color = 'red')

6.1.1 Giải thích kết quả sau khi chạy mô hình

Hàm hồi quy: \(logit(\pi) = log\left(\frac{\pi}{1-\pi}\right) = 0.776614 - 0.008117*hp\)

Từ mô hình logit ta suy ra được:

  • Nếu hồi qui đơn: \[\pi = \frac{e^{\beta_0 + \beta_1x_1}}{1 + e^{\beta_0 + \beta_1x_1}}\]

  • Nếu hồi quy bội: \[\pi = \frac{e^{\beta_0 + \beta_1x_1 + \beta_2x_2 +\dots+\beta_kx_k}}{1 + e^{\beta_0 + \beta_1x_1 + \beta_2x_2 +\dots+\beta_kx_k}}\]

Nếu chúng ta muốn sử dụng mô hình sau khi ước lượng để làm bài toán dự báo thì chúng ta thực hiện như sau: Tạo mới 1 bộ dữ liệu và thế vào hàm hồi quy vừa tìm được.

newdata = data.frame(hp=c(100, 90, 108, 90, 80, 90, 80, 90))

re <- predict(object = fitted, newdata = newdata)
round(re,4)
##       1       2       3       4       5       6       7       8 
## -0.0351  0.0461 -0.1000  0.0461  0.1272  0.0461  0.1272  0.0461
re <- predict(object = fitted, newdata = newdata, type = 'response')
round(re,4)
##      1      2      3      4      5      6      7      8 
## 0.4912 0.5115 0.4750 0.5115 0.5318 0.5115 0.5318 0.5115
fit <- glm(factor(am) ~ hp + mpg, family = binomial(link = 'logit'), data = mtcars)

summary(fit)
## 
## Call:
## glm(formula = factor(am) ~ hp + mpg, family = binomial(link = "logit"), 
##     data = mtcars)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -33.6052    15.0767   -2.23    0.026 *
## hp            0.0550     0.0269    2.05    0.041 *
## mpg           1.2596     0.5675    2.22    0.026 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43.230  on 31  degrees of freedom
## Residual deviance: 19.233  on 29  degrees of freedom
## AIC: 25.23
## 
## Number of Fisher Scoring iterations: 7
library(GLMsData)
data("turbines")

fit <- glm(data = turbines, formula = Fissures/Turbines ~ Hours, family = binomial(link = 'logit'), weights = Turbines)

fit
## 
## Call:  glm(formula = Fissures/Turbines ~ Hours, family = binomial(link = "logit"), 
##     data = turbines, weights = Turbines)
## 
## Coefficients:
## (Intercept)        Hours  
##   -3.923597     0.000999  
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:       113 
## Residual Deviance: 10.3  AIC: 49.8
Hour_n <- data.frame(Hours=c(1200,1350,1650))

predict(fit, Hour_n, type = 'response')
##       1       2       3 
## 0.06154 0.07079 0.09323
data("turbines")

glm(data = turbines, formula = cbind(Fissures, Turbines-Fissures) ~ Hours, family = binomial(link = 'logit'))
## 
## Call:  glm(formula = cbind(Fissures, Turbines - Fissures) ~ Hours, family = binomial(link = "logit"), 
##     data = turbines)
## 
## Coefficients:
## (Intercept)        Hours  
##   -3.923597     0.000999  
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:       113 
## Residual Deviance: 10.3  AIC: 49.8
fit
## 
## Call:  glm(formula = Fissures/Turbines ~ Hours, family = binomial(link = "logit"), 
##     data = turbines, weights = Turbines)
## 
## Coefficients:
## (Intercept)        Hours  
##   -3.923597     0.000999  
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:       113 
## Residual Deviance: 10.3  AIC: 49.8

Hàm liên kết cho dữ liệu nhị phân: logit, probit, cloglog, log, cauchit.

  1. \(logit(\pi) = log(\frac{\pi}{1-\pi})\)
  2. \(probit(\pi) = \Phi^{-1}(\pi)\)
  3. \(cloglog(\pi) = log(-log(1-\pi))\)

6.2 Các tiêu chí đánh giá mô hình

Để đánh giá các mô hình hồi cổ điển chúng ta thường dựa vào hệ số xác định mô hình (\(R^2\)), nhưng đối với mô các mô hình hồi quy tuyến tính tổng quát chúng ta sử dụng các tiêu chí sau:

6.2.1 AIC - Akaike Information Criterion

AIC được đề xuất bởi Akaike Hirotugu, một nhà thống kê học người Nhật. AIC là một tiêu chí được sử dụng một cách phổ biến để đánh giá một mô hình hồi quy được ước lượng bởi phương pháp Maximum Likekihood (ML). Một cách chung chung giá trị của AIC càng nhỏ thì mô hình càng tốt. AIC được tính bằng công thức sau: \[AIC = -2ln(L) + 2 k\] Với \(L\) là giá trị cực đại của hàm hợp lý (likelihood function) và \(k\) là số tham số của mô hình.

Khi thực hiện việc ước lượng mô hình hồi quy bằng lệnh glm thì chỉ số AIC đã được tính toán và thể hiện trên bảng kết quả ( bằng lệnh summary)

6.2.2 Deviance

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

lưu ý: Khi thực hiện việc ước lượng mô hình hồi quy bằng lệnh glm thì chỉ số AIC và Deviance đã được tính toán và thể hiện trên bảng kết quả ( bằng lệnh summary).

6.2.3 Brier Score

Là chỉ tiêu dùng để đánh giá mô hình hồi quy logistic, Brier Score được tính như sau:

\[B = \frac{1}{n}\sum_{i=1}^n(p_i - o_i)\] Trong đó: \(p_i, o_i\) lần lượt là giá trị xác suất quan sát được, và giá trị xác suất tính ra từ mô hình.

Ví dụ: Chúng ta thực hiện hàm hồi quy sau:

fit <- glm(factor(am) ~ hp + mpg, family = binomial(link = 'logit'), data = mtcars)

Sau đó chúng ta tính giá trị Brier Score như sau:

BrierScore(fit)
## [1] 0.1052

Giá trị của Brier Score càng nhỏ nghĩa là chênh lệch giữa xác suất thực tế và xác suất tính từ mô hình càng nhỏ, nghĩa là mô hình càng tốt.

6.2.4 Ma trận nhầm lẫn - Confusion matrix

Là ma trận thể hiện sự nhầm lẫn của mô hình.

Ví dụ chúng ta có mô hình sau:

fit <- glm(factor(am) ~ hp + mpg, family = binomial(link = 'logit'), data = mtcars)

và ma trận nhầm lẫn của mô hình này như sau:

library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following objects are masked from 'package:DescTools':
## 
##     MAE, RMSE
## The following object is masked from 'package:purrr':
## 
##     lift
confusionMatrix(table(predict(fit, type="response") >= 0.5,fit$data$am == '1'))
## Confusion Matrix and Statistics
## 
##        
##         FALSE TRUE
##   FALSE    16    3
##   TRUE      3   10
##                                         
##                Accuracy : 0.812         
##                  95% CI : (0.636, 0.928)
##     No Information Rate : 0.594         
##     P-Value [Acc > NIR] : 0.00757       
##                                         
##                   Kappa : 0.611         
##                                         
##  Mcnemar's Test P-Value : 1.00000       
##                                         
##             Sensitivity : 0.842         
##             Specificity : 0.769         
##          Pos Pred Value : 0.842         
##          Neg Pred Value : 0.769         
##              Prevalence : 0.594         
##          Detection Rate : 0.500         
##    Detection Prevalence : 0.594         
##       Balanced Accuracy : 0.806         
##                                         
##        'Positive' Class : FALSE         
##

Liên quan đến ma trận nhầm lẫn chúng ta có một số khái niệm sau:

6.3 Uớc lượng hàm hồi quy cho dữ liệu Poisson

Hồi quy Poisson có thể hồi quy cho: 1. Dữ liệu đếm (count data). 2. Dữ liệu dạng rate (rate data). 3. Tần số trong bảng ngẫu nhiên (contigency table).

Lưu ý rằng các loại dữ liệu này là dữ liệu có phân phối Poisson hoặc có liên quan đến phân phối Poisson.

Một số bài toán đặt ra:

  1. Chọn ngẫu nhiên 1 hồ sơ, tính xác suất để người này đã có \(k\) lần nhận thưởng.
  2. Chọn ngẫu nhiên 1 sinh viên, tính xác suất để người này đã rớt \(k\) môn học.
  3. Số lượng EMail rác phải nhận trong 1 tuần.
  4. Số hồ sơ công chứng không thành công trong 1 tháng.
  5. Số vụ cướp giật ở Sài Gòn trong 1 ngày.

Ví dụ sau chúng ta sẽ chạy mô hình hồi quy Poisson cho bộ dữ liệu về kết quả học tập, cụ thể chúng ta sẽ xây dựng mô hình thể hiện số lần trung bình nhận phần thưởng (count data) của 1 ứng cử viên ( thông qua một số đặc điểm).

dataset1 <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
# Dữ liệu này đã lưu thành file dataset1.csv
dataset1$prog <- factor(dataset1$prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
fit <- glm(data = dataset1, formula = num_awards ~ math, family = poisson(link = 'log'))
summary(fit)
## 
## Call:
## glm(formula = num_awards ~ math, family = poisson(link = "log"), 
##     data = dataset1)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.33353    0.59126   -9.02   <2e-16 ***
## math         0.08617    0.00968    8.90   <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: 287.67  on 199  degrees of freedom
## Residual deviance: 204.02  on 198  degrees of freedom
## AIC: 384.1
## 
## Number of Fisher Scoring iterations: 6

Ví dụ sau chúng sẽ chạy mô hình hồi quy Poisson cho bộ dữ liệu về cua móng ngựa. Mô hình này sẽ thể hiện số con cua đực trung bình (count data) bị thu hút bởi một con cua cái thông qua một số đặc điểm của con cua cái.

ds <- read.csv("Crab.csv", header = T)
model <- glm(data = ds, formula = satell ~ width + weight, family = poisson(link = 'log'))
summary(model)
## 
## Call:
## glm(formula = satell ~ width + weight, family = poisson(link = "log"), 
##     data = ds)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -1.295211   0.898896   -1.44   0.1496   
## width        0.046076   0.046750    0.99   0.3243   
## weight       0.000447   0.000159    2.82   0.0048 **
## ---
## 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: 559.90  on 170  degrees of freedom
## AIC: 921.2
## 
## Number of Fisher Scoring iterations: 6

Khi chạy mô hình hồi quy Poisson cho dữ liệu rate (rate data) chúng ta cần lưu ý:

  • \(Y\) = count (e.g., number violent acts).
  • \(t\) = index of the time or space (e.g., days in the community).
  • The sample rate of occurrence is \(Y/t\).
  • The expected value of the rate is \[E(Y_i/t) = \frac{1}{t}E(Y_i) = \mu_i/t\] Nên \[log(\mu/t) = \alpha + \beta x \\ log(\mu) - \log(t) = \alpha + \beta x \\ log(\mu) = \alpha + \beta x + log(t)\] \(log(t)\) gọi là độ dịch chuyển (offset).

Mô hình log - linear được viết lại như sau: \[log(\mu/t) = \alpha + \beta x \\ \mu/t = e^{\alpha}e^{\beta x} \\ \mu = te^{\alpha}e^{\beta x}\]