GLM trong R: Mô hình tuyến tính tổng quát với ví dụ
Evelyn Clarke
Hồi quy logistic là gì?
Hồi quy logistic được sử dụng để dự đoán một lớp, tức là xác suất. Hồi quy logistic có thể dự đoán chính xác kết quả nhị phân.
Hãy tưởng tượng bạn muốn dự đoán liệu một khoản vay có bị từ chối/chấp nhận hay không dựa trên nhiều thuộc tính. Hồi quy logistic có dạng 0/1. y = 0 nếu khoản vay bị từ chối, y = 1 nếu được chấp nhận.
Mô hình hồi quy logistic khác với mô hình hồi quy tuyến tính ở hai điểm.
- Trước hết, hồi quy logistic chỉ chấp nhận đầu vào nhị phân (nhị phân) làm biến phụ thuộc (tức là vectơ 0 và 1).
- Thứ hai, kết quả được đo bằng hàm liên kết xác suất sau đây được gọi là sigmoid do nó có hình chữ S.:
Đầu ra của hàm luôn nằm trong khoảng từ 0 đến 1. Kiểm tra hình ảnh bên dưới
Hàm sigmoid trả về các giá trị từ 0 đến 1. Đối với nhiệm vụ phân loại, chúng ta cần một đầu ra rời rạc là 0 hoặc 1.
Để chuyển đổi một luồng liên tục thành giá trị rời rạc, chúng ta có thể đặt giới hạn quyết định ở mức 0.5. Tất cả các giá trị trên ngưỡng này được phân loại là 1
Cách tạo Mô hình lót tổng quát (GLM)
Hãy sử dụng người lớn tập dữ liệu để minh họa hồi quy logistic. “Người lớn” là một tập dữ liệu tuyệt vời cho nhiệm vụ phân loại. Mục tiêu là dự đoán liệu thu nhập hàng năm tính bằng đô la của một cá nhân có vượt quá 50.000 hay không. Bộ dữ liệu chứa 46,033 quan sát và mười tính năng:
- age: tuổi của cá nhân. Số
- : Trình độ học vấn của cá nhân. Nhân tố.
- tình trạng hôn nhân: Maritrạng thái của cá nhân. Yếu tố tức là Chưa từng kết hôn, Đã kết hôn-vợ/chồng,…
- giới tính: Giới tính của cá nhân. Yếu tố, tức là Nam hay Nữ
- thu nhập = earnings: Target biến. Thu nhập trên hoặc dưới 50K. Hệ số tức là >50K, <=50K
giữa những người khác
library(dplyr)
data_adult <-read.csv("https://raw.githubusercontent.com/guru99-edu/R-Programming/master/adult.csv")
glimpse(data_adult)
Đầu ra:
Observations: 48,842 Variables: 10 $ x <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,... $ age <int> 25, 38, 28, 44, 18, 34, 29, 63, 24, 55, 65, 36, 26... $ workclass <fctr> Private, Private, Local-gov, Private, ?, Private,... $ education <fctr> 11th, HS-grad, Assoc-acdm, Some-college, Some-col... $ educational.num <int> 7, 9, 12, 10, 10, 6, 9, 15, 10, 4, 9, 13, 9, 9, 9,... $ marital.status <fctr> Never-married, Married-civ-spouse, Married-civ-sp... $ race <fctr> Black, White, White, Black, White, White, Black, ... $ gender <fctr> Male, Male, Male, Male, Female, Male, Male, Male,... $ hours.per.week <int> 40, 50, 40, 40, 30, 30, 40, 32, 40, 10, 40, 40, 39... $ income <fctr> <=50K, <=50K, >50K, >50K, <=50K, <=50K, <=50K, >5...
Chúng ta sẽ tiến hành như sau:
- Bước 1: Kiểm tra các biến liên tục
- Bước 2: Kiểm tra các biến hệ số
- Bước 3: Kỹ thuật tính năng
- Bước 4: Thống kê tóm tắt
- Bước 5: Tập huấn luyện/kiểm tra
- Bước 6: Xây dựng mô hình
- Bước 7: Đánh giá hiệu quả của mô hình
- bước 8: Cải thiện mô hình
Nhiệm vụ của bạn là dự đoán cá nhân nào sẽ có doanh thu cao hơn 50K.
Trong hướng dẫn này, mỗi bước sẽ được trình bày chi tiết để thực hiện phân tích trên tập dữ liệu thực.
Bước 1) Kiểm tra các biến liên tục
Ở bước đầu tiên, bạn có thể thấy sự phân bố của các biến liên tục.
continuous <-select_if(data_adult, is.numeric) summary(continuous)
Giải thích mã
- liên tục <- select_if(data_adult, is.numeric): Sử dụng hàm select_if() từ thư viện dplyr để chỉ chọn các cột số
- tóm tắt (liên tục): In thống kê tóm tắt
Đầu ra:
## X age educational.num hours.per.week ## Min. : 1 Min. :17.00 Min. : 1.00 Min. : 1.00 ## 1st Qu.:11509 1st Qu.:28.00 1st Qu.: 9.00 1st Qu.:40.00 ## Median :23017 Median :37.00 Median :10.00 Median :40.00 ## Mean :23017 Mean :38.56 Mean :10.13 Mean :40.95 ## 3rd Qu.:34525 3rd Qu.:47.00 3rd Qu.:13.00 3rd Qu.:45.00 ## Max. :46033 Max. :90.00 Max. :16.00 Max. :99.00
Từ bảng trên, bạn có thể thấy rằng dữ liệu có thang đo hoàn toàn khác nhau và số giờ.mỗi tuần có các giá trị ngoại lệ lớn (.tức là hãy xem giá trị tối đa và tứ phân vị cuối cùng).
Bạn có thể giải quyết vấn đề này theo hai bước sau:
- 1: Vẽ biểu đồ phân bổ số giờ mỗi tuần
- 2: Chuẩn hóa các biến liên tục
- Vẽ đồ thị phân phối
Chúng ta hãy xem xét kỹ hơn sự phân bổ số giờ.mỗi tuần
# Histogram with kernel density curve
library(ggplot2)
ggplot(continuous, aes(x = hours.per.week)) +
geom_density(alpha = .2, fill = "#FF6666")
Đầu ra:
Biến có nhiều ngoại lệ và phân phối không được xác định rõ. Bạn có thể giải quyết một phần vấn đề này bằng cách xóa 0.01% số giờ cao nhất mỗi tuần.
Cú pháp cơ bản của lượng tử:
quantile(variable, percentile) arguments: -variable: Select the variable in the data frame to compute the percentile -percentile: Can be a single value between 0 and 1 or multiple value. If multiple, use this format: `c(A,B,C, ...) - `A`,`B`,`C` and `...` are all integer from 0 to 1.
Chúng tôi tính toán phần trăm 2 hàng đầu
top_one_percent <- quantile(data_adult$hours.per.week, .99) top_one_percent
Giải thích mã
- quantile(data_adult$hours.per.week, .99): Tính giá trị của 99% thời gian làm việc
Đầu ra:
## 99% ## 80
98% dân số làm việc dưới 80 giờ mỗi tuần.
Bạn có thể bỏ các quan sát trên ngưỡng này. Bạn sử dụng bộ lọc từ dplyr thư viện.
data_adult_drop <-data_adult %>% filter(hours.per.week<top_one_percent) dim(data_adult_drop)
Đầu ra:
## [1] 45537 10
- Chuẩn hóa các biến liên tục
Bạn có thể chuẩn hóa từng cột để cải thiện hiệu suất vì dữ liệu của bạn không có cùng tỷ lệ. Bạn có thể sử dụng hàm mutate_if từ thư viện dplyr. Cú pháp cơ bản là:
mutate_if(df, condition, funs(function)) arguments: -`df`: Data frame used to compute the function - `condition`: Statement used. Do not use parenthesis - funs(function): Return the function to apply. Do not use parenthesis for the function
Bạn có thể chuẩn hóa các cột số như sau:
data_adult_rescale <- data_adult_drop % > % mutate_if(is.numeric, funs(as.numeric(scale(.)))) head(data_adult_rescale)
Giải thích mã
- mutate_if(is.numeric, funs(scale)): Điều kiện chỉ là cột số và hàm là tỷ lệ
Đầu ra:
## X age workclass education educational.num ## 1 -1.732680 -1.02325949 Private 11th -1.22106443 ## 2 -1.732605 -0.03969284 Private HS-grad -0.43998868 ## 3 -1.732530 -0.79628257 Local-gov Assoc-acdm 0.73162494 ## 4 -1.732455 0.41426100 Private Some-college -0.04945081 ## 5 -1.732379 -0.34232873 Private 10th -1.61160231 ## 6 -1.732304 1.85178149 Self-emp-not-inc Prof-school 1.90323857 ## marital.status race gender hours.per.week income ## 1 Never-married Black Male -0.03995944 <=50K ## 2 Married-civ-spouse White Male 0.86863037 <=50K ## 3 Married-civ-spouse White Male -0.03995944 >50K ## 4 Married-civ-spouse Black Male -0.03995944 >50K ## 5 Never-married White Male -0.94854924 <=50K ## 6 Married-civ-spouse White Male -0.76683128 >50K
Bước 2) Kiểm tra các biến hệ số
Bước này có hai mục tiêu:
- Kiểm tra cấp độ trong từng cột phân loại
- Xác định cấp độ mới
Chúng ta sẽ chia bước này thành ba phần:
- Chọn các cột phân loại
- Lưu trữ biểu đồ thanh của từng cột trong danh sách
- In đồ thị
Chúng ta có thể chọn các cột yếu tố với mã bên dưới:
# Select categorical column factor <- data.frame(select_if(data_adult_rescale, is.factor)) ncol(factor)
Giải thích mã
- data.frame(select_if(data_adult, is.factor)): Chúng tôi lưu trữ các cột hệ số trong hệ số trong một loại khung dữ liệu. Thư viện ggplot2 yêu cầu đối tượng khung dữ liệu.
Đầu ra:
## [1] 6
Tập dữ liệu chứa 6 biến phân loại
Bước thứ hai khéo léo hơn. Bạn muốn vẽ biểu đồ thanh cho từng cột trong hệ số khung dữ liệu. Sẽ thuận tiện hơn khi tự động hóa quy trình, đặc biệt trong trường hợp có nhiều cột.
library(ggplot2)
# Create graph for each column
graph <- lapply(names(factor),
function(x)
ggplot(factor, aes(get(x))) +
geom_bar() +
theme(axis.text.x = element_text(angle = 90)))
Giải thích mã
- lapply(): Sử dụng hàm lapply() để truyền một hàm vào tất cả các cột của tập dữ liệu. Bạn lưu trữ đầu ra trong một danh sách
- function(x): Hàm sẽ được xử lý cho từng x. Ở đây x là các cột
- ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): Tạo biểu đồ thanh char cho mỗi phần tử x. Lưu ý, để trả về x dưới dạng một cột, bạn cần đưa nó vào trong get()
Bước cuối cùng tương đối dễ dàng. Bạn muốn in 6 biểu đồ.
# Print the graph graph
Đầu ra:
## [[1]]
## ## [[2]]
## ## [[3]]
## ## [[4]]
## ## [[5]]
## ## [[6]]
Lưu ý: Sử dụng nút tiếp theo để điều hướng đến biểu đồ tiếp theo
Bước 3) Kỹ thuật tính năng
Làm lại nền giáo dục
Từ biểu đồ trên, bạn có thể thấy rằng trình độ học vấn có thể thay đổi có 16 cấp độ. Điều này là đáng kể và một số cấp độ có số lượng quan sát tương đối thấp. Nếu bạn muốn cải thiện lượng thông tin bạn có thể nhận được từ biến này, bạn có thể chuyển nó lên cấp độ cao hơn. Cụ thể là bạn tạo các nhóm lớn hơn với trình độ học vấn tương tự. Ví dụ, trình độ học vấn thấp sẽ chuyển thành bỏ học. Các cấp học cao hơn sẽ được đổi thành thạc sĩ.
Đây là chi tiết:
| Cấp độ cũ | Cấp độ mới |
|---|---|
| Trường mầm non | rơi ra ngoài |
| 10th | Rơi ra ngoài |
| 11th | Rơi ra ngoài |
| 12th | Rơi ra ngoài |
| 1-4 | Rơi ra ngoài |
| 5th-6th | Rơi ra ngoài |
| 7th-8th | Rơi ra ngoài |
| 9th | Rơi ra ngoài |
| HS-Grad | Cao cấp |
| Một số trường đại học | Cộng đồng |
| PGS-acdm | Cộng đồng |
| PGS-Voc | Cộng đồng |
| Cử nhân | Cử nhân |
| Thạc sĩ | Thạc sĩ |
| trường chuyên nghiệp | Thạc sĩ |
| sĩ | Bằng Tiến sĩ |
recast_data <- data_adult_rescale % > %
select(-X) % > %
mutate(education = factor(ifelse(education == "Preschool" | education == "10th" | education == "11th" | education == "12th" | education == "1st-4th" | education == "5th-6th" | education == "7th-8th" | education == "9th", "dropout", ifelse(education == "HS-grad", "HighGrad", ifelse(education == "Some-college" | education == "Assoc-acdm" | education == "Assoc-voc", "Community",
ifelse(education == "Bachelors", "Bachelors",
ifelse(education == "Masters" | education == "Prof-school", "Master", "PhD")))))))
Giải thích mã
- Chúng tôi sử dụng động từ đột biến từ thư viện dplyr. Chúng ta thay đổi các giá trị của giáo dục bằng câu lệnh ifelse
Trong bảng bên dưới, bạn tạo một thống kê tóm tắt để xem trung bình cần bao nhiêu năm học (giá trị z) để đạt được Cử nhân, Thạc sĩ hoặc Tiến sĩ.
recast_data % > % group_by(education) % > % summarize(average_educ_year = mean(educational.num), count = n()) % > % arrange(average_educ_year)
Đầu ra:
## # A tibble: 6 x 3 ## education average_educ_year count ## <fctr> <dbl> <int> ## 1 dropout -1.76147258 5712 ## 2 HighGrad -0.43998868 14803 ## 3 Community 0.09561361 13407 ## 4 Bachelors 1.12216282 7720 ## 5 Master 1.60337381 3338 ## 6 PhD 2.29377644 557
Đúc lại Maritrạng thái tal
Bạn cũng có thể tạo các cấp độ thấp hơn cho trạng thái hôn nhân. Trong đoạn mã sau, bạn thay đổi cấp độ như sau:
| Cấp độ cũ | Cấp độ mới |
|---|---|
| Không bao giờ kết hôn | Chưa kết hôn |
| Đã kết hôn-vợ chồng-vắng mặt | Chưa kết hôn |
| Đã kết hôn-AF-vợ/chồng | Kết hôn |
| Đã kết hôn-công dân-vợ / chồng | |
| Ly thân | Ly thân |
| Разведенный | |
| góa phụ | Widow |
# Change level marry recast_data <- recast_data % > % mutate(marital.status = factor(ifelse(marital.status == "Never-married" | marital.status == "Married-spouse-absent", "Not_married", ifelse(marital.status == "Married-AF-spouse" | marital.status == "Married-civ-spouse", "Married", ifelse(marital.status == "Separated" | marital.status == "Divorced", "Separated", "Widow")))))
Bạn có thể kiểm tra số lượng cá nhân trong mỗi nhóm.
table(recast_data$marital.status)
Đầu ra:
## ## Married Not_married Separated Widow ## 21165 15359 7727 1286
Bước 4) Thống kê tóm tắt
Đã đến lúc kiểm tra một số số liệu thống kê về các biến mục tiêu của chúng ta. Trong biểu đồ bên dưới, bạn đếm phần trăm cá nhân kiếm được hơn 50 nghìn đô la theo giới tính của họ.
# Plot gender income
ggplot(recast_data, aes(x = gender, fill = income)) +
geom_bar(position = "fill") +
theme_classic()
Đầu ra:
Tiếp theo, hãy kiểm tra xem nguồn gốc của cá nhân có ảnh hưởng đến thu nhập của họ hay không.
# Plot origin income
ggplot(recast_data, aes(x = race, fill = income)) +
geom_bar(position = "fill") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90))
Đầu ra:
Số giờ làm việc theo giới tính.
# box plot gender working time
ggplot(recast_data, aes(x = gender, y = hours.per.week)) +
geom_boxplot() +
stat_summary(fun.y = mean,
geom = "point",
size = 3,
color = "steelblue") +
theme_classic()
Đầu ra:
Biểu đồ hộp xác nhận rằng sự phân bổ thời gian làm việc phù hợp với các nhóm khác nhau. Trong biểu đồ hộp, cả hai giới đều không có quan sát đồng nhất.
Bạn có thể kiểm tra mật độ thời gian làm việc hàng tuần theo loại hình giáo dục. Các bản phân phối có nhiều lựa chọn khác nhau. Có lẽ có thể giải thích được là do loại hợp đồng ở Mỹ.
# Plot distribution working time by education
ggplot(recast_data, aes(x = hours.per.week)) +
geom_density(aes(color = education), alpha = 0.5) +
theme_classic()
Giải thích mã
- ggplot(recast_data, aes( x=hours.per.week)): Biểu đồ mật độ chỉ yêu cầu một biến
- geom_d mật(aes(color = education), alpha =0.5): Đối tượng hình học để kiểm soát mật độ
Đầu ra:
Để xác nhận suy nghĩ của mình, bạn có thể thực hiện thao tác một chiều Kiểm định ANOVA:
anova <- aov(hours.per.week~education, recast_data) summary(anova)
Đầu ra:
## Df Sum Sq Mean Sq F value Pr(>F) ## education 5 1552 310.31 321.2 <2e-16 *** ## Residuals 45531 43984 0.97 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Kiểm định ANOVA xác nhận sự khác biệt về giá trị trung bình giữa các nhóm.
Phi tuyến tính
Trước khi chạy mô hình, bạn có thể xem số giờ làm việc có liên quan đến tuổi tác hay không.
library(ggplot2)
ggplot(recast_data, aes(x = age, y = hours.per.week)) +
geom_point(aes(color = income),
size = 0.5) +
stat_smooth(method = 'lm',
formula = y~poly(x, 2),
se = TRUE,
aes(color = income)) +
theme_classic()
Giải thích mã
- ggplot(recast_data, aes(x = age, y =hours.per.week)): Đặt tính thẩm mỹ của biểu đồ
- geom_point(aes(color= thu nhập), size =0.5): Xây dựng biểu đồ chấm
- stat_smooth(): Thêm đường xu hướng với các đối số sau:
- Method='lm': Vẽ biểu đồ giá trị phù hợp nếu hồi quy tuyến tính
- công thức = y~poly(x,2): Phù hợp với hồi quy đa thức
- se = TRUE: Thêm sai số chuẩn
- aes(color= thu nhập): Chia mô hình theo thu nhập
Đầu ra:
Tóm lại, bạn có thể kiểm tra các thuật ngữ tương tác trong mô hình để tìm ra hiệu ứng phi tuyến tính giữa thời gian làm việc hàng tuần và các tính năng khác. Điều quan trọng là phải phát hiện thời gian làm việc khác nhau trong điều kiện nào.
Tương quan
Kiểm tra tiếp theo là hình dung mối tương quan giữa các biến. Bạn chuyển đổi loại cấp độ yếu tố thành số để có thể vẽ bản đồ nhiệt chứa hệ số tương quan được tính toán bằng phương pháp Spearman.
library(GGally)
# Convert data to numeric
corr <- data.frame(lapply(recast_data, as.integer))
# Plot the graphggcorr(corr,
method = c("pairwise", "spearman"),
nbreaks = 6,
hjust = 0.8,
label = TRUE,
label_size = 3,
color = "grey50")
Giải thích mã
- data.frame(lapply(recast_data,as.integer)): Chuyển đổi dữ liệu thành số
- ggcorr() vẽ bản đồ nhiệt với các đối số sau:
- phương pháp: Phương pháp tính toán mối tương quan
- nbreaks = 6: Số lần ngắt
- hjust = 0.8: Vị trí điều khiển của tên biến trong đồ thị
- nhãn = ĐÚNG: Thêm nhãn vào giữa cửa sổ
- label_size = 3: Nhãn kích thước
- color=“grey50”): Màu sắc của nhãn
Đầu ra:
Bước 5) Tập huấn luyện/kiểm tra
Bất kỳ sự giám sát nào học máy nhiệm vụ yêu cầu phân chia dữ liệu giữa tập tàu và tập kiểm tra. Bạn có thể sử dụng “hàm” mà bạn đã tạo trong các hướng dẫn học có giám sát khác để tạo tập huấn luyện/kiểm tra.
set.seed(1234)
create_train_test <- function(data, size = 0.8, train = TRUE) {
n_row = nrow(data)
total_row = size * n_row
train_sample <- 1: total_row
if (train == TRUE) {
return (data[train_sample, ])
} else {
return (data[-train_sample, ])
}
}
data_train <- create_train_test(recast_data, 0.8, train = TRUE)
data_test <- create_train_test(recast_data, 0.8, train = FALSE)
dim(data_train)
Đầu ra:
## [1] 36429 9
dim(data_test)
Đầu ra:
## [1] 9108 9
Bước 6) Xây dựng mô hình
Để xem thuật toán hoạt động như thế nào, bạn sử dụng gói glm(). Các Mô hình tuyến tính tổng quát là một tập hợp các mô hình Cú pháp cơ bản là:
glm(formula, data=data, family=linkfunction() Argument: - formula: Equation used to fit the model- data: dataset used - Family: - binomial: (link = "logit") - gaussian: (link = "identity") - Gamma: (link = "inverse") - inverse.gaussian: (link = "1/mu^2") - poisson: (link = "log") - quasi: (link = "identity", variance = "constant") - quasibinomial: (link = "logit") - quasipoisson: (link = "log")
Bạn đã sẵn sàng ước tính mô hình logistic để phân chia mức thu nhập giữa một tập hợp các tính năng.
formula <- income~. logit <- glm(formula, data = data_train, family = 'binomial') summary(logit)
Giải thích mã
- công thức <- thu nhập ~.: Tạo mô hình cho phù hợp
- logit <- glm(công thức, data = data_train, family = 'binomial'): Điều chỉnh mô hình logistic (family = 'binomial') với dữ liệu data_train.
- summary(logit): In tóm tắt của mô hình
Đầu ra:
## ## Call: ## glm(formula = formula, family = "binomial", data = data_train) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.6456 -0.5858 -0.2609 -0.0651 3.1982 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 0.07882 0.21726 0.363 0.71675 ## age 0.41119 0.01857 22.146 < 2e-16 *** ## workclassLocal-gov -0.64018 0.09396 -6.813 9.54e-12 *** ## workclassPrivate -0.53542 0.07886 -6.789 1.13e-11 *** ## workclassSelf-emp-inc -0.07733 0.10350 -0.747 0.45499 ## workclassSelf-emp-not-inc -1.09052 0.09140 -11.931 < 2e-16 *** ## workclassState-gov -0.80562 0.10617 -7.588 3.25e-14 *** ## workclassWithout-pay -1.09765 0.86787 -1.265 0.20596 ## educationCommunity -0.44436 0.08267 -5.375 7.66e-08 *** ## educationHighGrad -0.67613 0.11827 -5.717 1.08e-08 *** ## educationMaster 0.35651 0.06780 5.258 1.46e-07 *** ## educationPhD 0.46995 0.15772 2.980 0.00289 ** ## educationdropout -1.04974 0.21280 -4.933 8.10e-07 *** ## educational.num 0.56908 0.07063 8.057 7.84e-16 *** ## marital.statusNot_married -2.50346 0.05113 -48.966 < 2e-16 *** ## marital.statusSeparated -2.16177 0.05425 -39.846 < 2e-16 *** ## marital.statusWidow -2.22707 0.12522 -17.785 < 2e-16 *** ## raceAsian-Pac-Islander 0.08359 0.20344 0.411 0.68117 ## raceBlack 0.07188 0.19330 0.372 0.71001 ## raceOther 0.01370 0.27695 0.049 0.96054 ## raceWhite 0.34830 0.18441 1.889 0.05894 . ## genderMale 0.08596 0.04289 2.004 0.04506 * ## hours.per.week 0.41942 0.01748 23.998 < 2e-16 *** ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 40601 on 36428 degrees of freedom ## Residual deviance: 27041 on 36406 degrees of freedom ## AIC: 27087 ## ## Number of Fisher Scoring iterations: 6
Bản tóm tắt mô hình của chúng tôi tiết lộ thông tin thú vị. Hiệu suất của hồi quy logistic được đánh giá bằng các số liệu chính cụ thể.
- AIC (Tiêu chí thông tin Akaike): Điều này tương đương với R2 trong hồi quy logistic. Nó đo lường sự phù hợp khi áp dụng hình phạt đối với số lượng tham số. Nhỏ hơn AIC giá trị cho thấy mô hình gần với sự thật hơn.
- Độ lệch null: Chỉ phù hợp với mô hình với phần chặn. Mức độ tự do là n-1. Chúng ta có thể hiểu nó là giá trị Chi bình phương (giá trị phù hợp khác với thử nghiệm giả thuyết giá trị thực tế).
- Độ lệch dư: Mô hình với tất cả các biến. Nó cũng được hiểu là kiểm tra giả thuyết Chi bình phương.
- Number of Fisher Scoring iterations: Số lần lặp trước khi hội tụ.
Đầu ra của hàm glm() được lưu trữ trong một danh sách. Mã bên dưới hiển thị tất cả các mục có sẵn trong biến logit mà chúng tôi đã xây dựng để đánh giá hồi quy logistic.
# Danh sách rất dài, chỉ in 3 phần tử đầu tiên
lapply(logit, class)[1:3]
Đầu ra:
## $coefficients ## [1] "numeric" ## ## $residuals ## [1] "numeric" ## ## $fitted.values ## [1] "numeric"
Mỗi giá trị có thể được trích xuất bằng ký hiệu $ theo sau là tên của số liệu. Ví dụ: bạn đã lưu mô hình dưới dạng logit. Để trích xuất tiêu chí AIC, bạn sử dụng:
logit$aic
Đầu ra:
## [1] 27086.65
Bước 7) Đánh giá hiệu quả của mô hình
Ma trận hỗn loạn
ma trận hỗn loạn là lựa chọn tốt hơn để đánh giá hiệu suất phân loại so với các số liệu khác nhau mà bạn đã thấy trước đây. Ý tưởng chung là đếm số lần các trường hợp Đúng được phân loại là Sai.
Để tính toán ma trận nhầm lẫn, trước tiên bạn cần có một bộ dự đoán để có thể so sánh chúng với các mục tiêu thực tế.
predict <- predict(logit, data_test, type = 'response') # confusion matrix table_mat <- table(data_test$income, predict > 0.5) table_mat
Giải thích mã
- dự đoán(logit,data_test, type='response'): Tính toán dự đoán trên tập kiểm tra. Đặt type = 'response' để tính xác suất phản hồi.
- table(data_test$thu nhập, dự đoán > 0.5): Tính toán ma trận nhầm lẫn. dự đoán > 0.5 có nghĩa là nó trả về 1 nếu xác suất dự đoán trên 0.5, nếu không thì 0.
Đầu ra:
## ## FALSE TRUE ## <=50K 6310 495 ## >50K 1074 1229
Mỗi hàng trong ma trận nhầm lẫn đại diện cho một mục tiêu thực tế, trong khi mỗi cột đại diện cho một mục tiêu được dự đoán. Hàng đầu tiên của ma trận này coi thu nhập thấp hơn 50k (loại Sai): 6241 được phân loại chính xác là cá nhân có thu nhập thấp hơn 50k (Âm tính thật), trong khi phần còn lại bị phân loại sai là trên 50k (Dương tính giả). Hàng thứ 50 xét thu nhập trên 1229k thì lớp dương là XNUMX (Đúng tích cực), trong khi Âm tính thật là 1074.
Bạn có thể tính toán mô hình chính xác bằng cách tính tổng giá trị dương thực + âm thực trên tổng số quan sát
accuracy_Test <- sum(diag(table_mat)) / sum(table_mat) accuracy_Test
Giải thích mã
- sum(diag(table_mat)): Tổng của đường chéo
- sum(table_mat): Tổng của ma trận.
Đầu ra:
## [1] 0.8277339
Mô hình dường như gặp phải một vấn đề, đó là nó đánh giá quá cao số lượng kết quả âm tính giả. Đây được gọi là nghịch lý kiểm tra độ chính xác. Chúng tôi đã tuyên bố rằng độ chính xác là tỷ lệ dự đoán đúng trên tổng số trường hợp. Chúng ta có thể có độ chính xác tương đối cao nhưng lại là một mô hình vô dụng. Nó xảy ra khi có một giai cấp thống trị. Nếu nhìn lại ma trận nhầm lẫn, bạn có thể thấy hầu hết các trường hợp đều được phân loại là âm tính thực sự. Bây giờ hãy tưởng tượng, mô hình đã phân loại tất cả các lớp là âm (tức là thấp hơn 50k). Bạn sẽ có độ chính xác là 75 phần trăm (6718/6718+2257). Mô hình của bạn hoạt động tốt hơn nhưng gặp khó khăn trong việc phân biệt mặt tích cực thực sự với mặt tiêu cực thực sự.
Trong tình huống như vậy, tốt nhất nên có số liệu ngắn gọn hơn. Chúng ta có thể nhìn vào:
- Độ chính xác=TP/(TP+FP)
- Thu hồi=TP/(TP+FN)
Độ chính xác và thu hồi
Độ chính xác xem xét tính chính xác của dự đoán tích cực. Nhớ lại là tỷ lệ các trường hợp tích cực được bộ phân loại phát hiện chính xác;
Bạn có thể xây dựng hai hàm để tính hai số liệu này
- Xây dựng độ chính xác
precision <- function(matrix) {
# True positive
tp <- matrix[2, 2]
# false positive
fp <- matrix[1, 2]
return (tp / (tp + fp))
}
Giải thích mã
- mat[1,1]: Trả về ô đầu tiên của cột đầu tiên của khung dữ liệu, tức là số dương thực sự
- chiếu[1,2]; Trả về ô đầu tiên của cột thứ hai của khung dữ liệu, tức là dương tính giả
recall <- function(matrix) {
# true positive
tp <- matrix[2, 2]# false positive
fn <- matrix[2, 1]
return (tp / (tp + fn))
}
Giải thích mã
- mat[1,1]: Trả về ô đầu tiên của cột đầu tiên của khung dữ liệu, tức là số dương thực sự
- chiếu[2,1]; Trả về ô thứ hai của cột đầu tiên của khung dữ liệu, tức là âm tính giả
Bạn có thể kiểm tra chức năng của mình
prec <- precision(table_mat) prec rec <- recall(table_mat) rec
Đầu ra:
## [1] 0.712877 ## [2] 0.5336518
Khi mô hình cho biết đó là một cá nhân trên 50k, nó chỉ đúng trong 54% trường hợp và có thể yêu cầu những cá nhân trên 50k trong 72% trường hợp.
Bạn có thể tạo 
cho điểm dựa trên độ chính xác và thu hồi. Các là giá trị trung bình hài hòa của hai số liệu này, có nghĩa là nó mang lại nhiều trọng số hơn cho các giá trị thấp hơn.
f1 <- 2 * ((prec * rec) / (prec + rec)) f1
Đầu ra:
## [1] 0.6103799
Sự cân bằng giữa độ chính xác và khả năng thu hồi
Không thể có cả độ chính xác cao và khả năng thu hồi cao.
Nếu chúng ta tăng độ chính xác, thì dự đoán đúng cá nhân sẽ tốt hơn, nhưng chúng ta sẽ bỏ lỡ rất nhiều cá thể trong số đó (tỷ lệ thu hồi thấp hơn). Trong một số trường hợp, chúng tôi muốn độ chính xác cao hơn mức thu hồi. Có một mối quan hệ lõm giữa độ chính xác và thu hồi.
- Hãy tưởng tượng, bạn cần dự đoán xem một bệnh nhân có mắc bệnh hay không. Bạn muốn càng chính xác càng tốt.
- Nếu bạn cần phát hiện những người có khả năng lừa đảo trên đường phố thông qua nhận dạng khuôn mặt, sẽ tốt hơn nếu bạn bắt được nhiều người bị gắn mác lừa đảo mặc dù độ chính xác thấp. Cảnh sát sẽ có thể trả tự do cho cá nhân không lừa đảo.
Đường cong ROC
Người nhận Operađặc điểm ting đường cong là một công cụ phổ biến khác được sử dụng với phân loại nhị phân. Nó rất giống với đường cong độ chính xác/thu hồi, nhưng thay vì vẽ đồ thị độ chính xác so với thu hồi, đường cong ROC hiển thị tỷ lệ dương thực sự (tức là thu hồi) so với tỷ lệ dương tính giả. Tỷ lệ dương tính giả là tỷ lệ các trường hợp âm tính được phân loại không chính xác thành dương tính. Nó bằng một trừ đi tỷ lệ âm thực sự. Tỷ lệ âm thực sự còn được gọi là tính cụ thể. Do đó, đồ thị đường cong ROC nhạy cảm (thu hồi) so với 1-đặc hiệu
Để vẽ đường cong ROC, chúng ta cần cài đặt thư viện có tên RORC. Chúng ta có thể tìm thấy trong conda thư viện. Bạn có thể gõ mã:
cài đặt conda -cr r-rocr –có
Chúng ta có thể vẽ biểu đồ ROC bằng các hàm dự đoán() và hiệu suất().
library(ROCR) ROCRpred <- prediction(predict, data_test$income) ROCRperf <- performance(ROCRpred, 'tpr', 'fpr') plot(ROCRperf, colorize = TRUE, text.adj = c(-0.2, 1.7))
Giải thích mã
- dự đoán(dự đoán, data_test$thu nhập): Thư viện ROCR cần tạo đối tượng dự đoán để chuyển đổi dữ liệu đầu vào
- performance(ROCRpred, 'tpr','fpr'): Trả về hai kết hợp để tạo ra trong biểu đồ. Ở đây, tpr và fpr được xây dựng. Tot vẽ độ chính xác và thu hồi cùng nhau, sử dụng “prec”, “rec”.
Đầu ra:
Bước 8) Cải thiện mô hình
Bạn có thể thử thêm tính phi tuyến tính vào mô hình bằng sự tương tác giữa
- tuổi và giờ.per.week
- giới tính và số giờ.per.week.
Bạn cần sử dụng bài kiểm tra điểm để so sánh cả hai mô hình
formula_2 <- income~age: hours.per.week + gender: hours.per.week + . logit_2 <- glm(formula_2, data = data_train, family = 'binomial') predict_2 <- predict(logit_2, data_test, type = 'response') table_mat_2 <- table(data_test$income, predict_2 > 0.5) precision_2 <- precision(table_mat_2) recall_2 <- recall(table_mat_2) f1_2 <- 2 * ((precision_2 * recall_2) / (precision_2 + recall_2)) f1_2
Đầu ra:
## [1] 0.6109181
Điểm số cao hơn một chút so với lần trước. Bạn có thể tiếp tục làm việc với dữ liệu và cố gắng đạt được điểm số.
Tổng kết
Chúng ta có thể tóm tắt hàm huấn luyện hồi quy logistic trong bảng dưới đây:
| Bưu kiện | Mục tiêu | Chức năng | Tranh luận |
|---|---|---|---|
| – | Tạo tập dữ liệu đào tạo/kiểm tra | create_train_set() | dữ liệu, kích thước, đào tạo |
| lém lỉnh | Huấn luyện một mô hình tuyến tính tổng quát | glm() | công thức, dữ liệu, họ* |
| lém lỉnh | Tóm tắt mô hình | bản tóm tắt() | mô hình vừa vặn |
| cơ sở | Đưa ra dự đoán | dự đoán () | mô hình, tập dữ liệu, loại = 'phản hồi' được trang bị |
| cơ sở | Tạo một ma trận nhầm lẫn | bàn() | y, dự đoán() |
| cơ sở | Tạo điểm chính xác | tổng(diag(bảng())/tổng(bảng() | |
| ROCR | Tạo ROC: Bước 1 Tạo dự đoán | sự dự đoán() | dự đoán(), y |
| ROCR | Tạo ROC: Bước 2 Tạo hiệu suất | hiệu suất() | dự đoán(), 'tpr', 'fpr' |
| ROCR | Tạo ROC: Bước 3 Vẽ đồ thị | âm mưu() | hiệu suất() |
Cai khac GLM loại mô hình là:
– nhị thức: (link=“logit”)
– gaussian: (link = “danh tính”)
– Gamma: (liên kết = “nghịch đảo”)
– inverse.gaussian: (link = “1/mu^2”)
– poisson: (link = “log”)
– gần đúng: (link = “danh tính”, phương sai = “hằng số”)
– quasibinomial: (link = “logit”)
– quasipoisson: (link = “log”)
