GLM trong R: Mô hình tuyến tính tổng quát và hồi quy logistic

Mô hình tuyến tính tổng quát (GLM) trong R là gì?

A Mô hình tuyến tính tổng quát (GLM) Phương pháp này mở rộng hồi quy tuyến tính thông thường để biến phản hồi có thể tuân theo phân phối khác với phân phối chuẩn. Trong RBạn có thể lắp một cái có sẵn. glm() Hàm từ gói stats.

Mỗi mô hình GLM được xác định bởi ba thành phần:

• Thành phần ngẫu nhiên: Phân phối xác suất của biến phản hồi, được lấy từ họ phân phối mũ (nhị thức, Poisson, Gamma, Gaussian và các họ khác).
• Thành phần hệ thống: Mô hình dự báo tuyến tính, tức là sự kết hợp có trọng số của các biến giải thích của bạn.
• Chức năng liên kết: Hàm này kết nối giá trị trung bình của biến phản hồi với biến dự báo tuyến tính, ví dụ như hàm liên kết logit cho dữ liệu nhị phân hoặc hàm liên kết logarit cho dữ liệu đếm.

Cấu trúc đó chính là điều làm cho mô hình trở nên "tổng quát". Thay vì ép buộc một kết quả phân bố chuẩn, bạn khai báo phân bố chính xác thông qua... gia đình R ước tính các hệ số bằng phương pháp hợp lý tối đa. Hồi quy logistic đơn giản là một GLM với họ phân phối nhị thức và hàm liên kết logit, vì vậy nó là điểm khởi đầu tự nhiên.

Hồi quy Logistic trong R 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 thể hiện thông qua một hàm liên kết xác suất được gọi là... sigmoid Chức năng (logistics) do hình dạng chữ S của nó:

Đầ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

Giờ thì hàm liên kết đã rõ ràng, hãy so sánh mô hình tổng quát với mô hình tuyến tính thông thường mà bạn đã biết.

GLM so với Hồi quy tuyến tính: Những điểm khác biệt chính về R²

Trước khi viết bất kỳ đoạn mã nào, việc biết chính xác thời điểm sử dụng chuẩn là rất hữu ích. hồi quy tuyến tính Hàm lm() không còn phù hợp nữa và nên thay thế nó bằng hàm glm().

Tóm lại, hãy chọn lm() khi kết quả là một phép đo phân bố chuẩn và glm() bất cứ khi nào kết quả là một quyết định có/không, một nhiệm vụ mà bạn có thể giao cho một hàm khác. cây quyết định Bộ phân loại, số lượng sự kiện, hoặc một đại lượng dương có độ phân tán tăng theo giá trị trung bình của nó.

Hướng dẫn cách tạo mô hình tuyến tính tổng quát (GLM) trong R

Sau khi đã nắm vững lý thuyết, phần còn lại của hướng dẫn này sẽ áp dụng mô hình GLM nhị thức từ đầu đến cuối trên một tập dữ liệu thực tế.

Hãy sử dụng người lớn Bộ dữ liệu này được sử dụng để minh họa hồi quy Logistic. Bộ dữ liệu "người trưởng thành" là một bộ 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 của một cá nhân tính bằng đô la Mỹ có vượt quá 50,000 hay không. Bộ dữ liệu chứa 48,842 quan sát và mười biến:

• 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)

Code Giải thích

• 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:

• Vẽ biểu đồ phân bố số giờ làm việc mỗi tuần.
• Chuẩn hóa các biến liên tục

1. 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 vị thứ 99 của số giờ làm việc hàng tuần.

top_one_percent <- quantile(data_adult$hours.per.week, .99)
top_one_percent

Code Giải thích

• quantile(data_adult$hours.per.week, .99): Tính phân vị thứ 99 của thời gian làm việc hàng tuần

Đầu ra:

## 99% 
##  80

99% 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

2. 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)

Code Giải thích

• 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)

Code Giải thích

• 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 phức tạp hơn. Bạn cần vẽ biểu đồ cột cho mỗi cột trong khung dữ liệu factor. Việc tự động hóa quy trình sẽ thuận tiện hơn, đặc biệt là 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)))

Code Giải thích

• 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

Hai biến phân loại có nhiều cấp độ hơn mức cần thiết cho mô hình. Bạn sẽ nhóm chúng lại thành các danh mục rộng hơn, có nhiều cấp độ hơn.

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:

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")))))))

Code Giải thích

• 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:

# 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. Phân bố có nhiều điểm khác biệt. Điều này có thể được giải thích bởi loại hình đào tạo.tracở 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()

Code Giải thích

• 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()

Code Giải thích

• 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")

Code Giải thích

• 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ụ này yêu cầu bạn chia dữ liệu thành tập huấn luyện và tập kiểm tra. Bạn có thể sử dụng "hàm" mà bạn đã tạo trong các bài hướng dẫn về 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 hàm glm() từ gói stats. 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)

Code Giải thích

• 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 biểu thị bằng...tracđược lưu trữ với dấu $ theo sau là tên của các chỉ số. Ví dụ: bạn đã lưu mô hình dưới dạng logit. Để ví dụtracTheo 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

Code Giải thích

• 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 biểu thị một mục tiêu thực tế, trong khi mỗi cột biểu thị một mục tiêu dự đoán. Hàng đầu tiên của ma trận này xem xét mức thu nhập thấp hơn 50 (lớp tiêu cực): 6,310 quan sát đã được phân loại chính xác là những cá nhân có thu nhập thấp hơn 50 (Âm tính thật), trong khi 495 trường hợp bị phân loại sai là trên 50 (Dương tính giả). Hàng thứ hai xem xét thu nhập trên 50: 1,229 người đã được xác định chính xác (Đúng tích cực), trong khi 1,074 người khác bị bỏ sót (Sai âm).

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

Code Giải thích

• 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 này dường như mắc phải một vấn đề: nó tạo ra quá nhiều kết quả âm tính giả. Điều này được gọi là nghịch lý kiểm tra độ chính xácChúng ta đã nói rằng độ chính xác là tỷ lệ giữa số dự đoán đúng và tổng số trường hợp. Chúng ta có thể có độ chính xác tương đối cao nhưng mô hình lại vô dụng. Điều này xảy ra khi có một lớp chiếm ưu thế. Nếu bạn 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 được phân loại là âm tính thật. Bây giờ hãy tưởng tượng, mô hình phân loại mọi quan sát là âm tính (tức là thấp hơn 50k). Bạn vẫn sẽ đạt được độ chính xác khoảng 75% (6,805 / 9,108). 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 dương tính thật với âm tính thật.

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

1. 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))
}

Code Giải thích

• 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))
}

Code Giải thích

• 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

Hãy đọc kỹ hai con số này. Độ chính xác (Precision) là 0.71, nghĩa là khi mô hình dự đoán một cá nhân có thu nhập trên 50 đô la, thì dự đoán này chính xác trong 71% trường hợp. Độ nhạy (Recall) là 0.53, nghĩa là mô hình chỉ phát hiện được 53% số cá nhân thực sự có thu nhập trên 50 đô la.

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 gói có tên là ROCR. Có thể tìm thấy gói này trong conda. thư viện. Bạn có thể gõ mã:

conda install -c r r-rocr --yes

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))

Code Giải thích

• 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.

Sau đó, bạn so sánh điểm F1 của 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 F1 cao hơn một chút so với lần trước. Bạn có thể tiếp tục xử lý dữ liệu và cố gắng vượt qua điểm số này.

Cách diễn giải hệ số GLM và tỷ lệ chênh lệch trong R

Bảng tóm tắt được in ở Bước 6 báo cáo các hệ số về log-tỷ lệ cược thang đo này rất khó giải thích cho người không chuyên về kỹ thuật. Việc chuyển đổi chúng thành tỷ lệ chênh lệch (odds ratio) giúp mô hình dễ truyền đạt hơn nhiều.

Hãy làm theo bốn bước sau.

1. Lũy thừa các hệ số. Áp dụng hàm exp() cho mọi ước tính để logarit tỷ lệ cược trở thành tỷ lệ cược nhân.
2. Thêm khoảng tin cậy. Bọc hàm confint() trong hàm exp() để thu được khoảng tin cậy 95% trên cùng thang tỷ lệ cược.
3. So sánh từng giá trị với 1. Tỷ lệ chênh lệch (odds ratio) trên 1 làm tăng xác suất của lớp tích cực, giá trị dưới 1 làm giảm xác suất, và giá trị gần bằng 1 có nghĩa là yếu tố dự đoán này không đóng góp nhiều.
4. Kiểm tra ý nghĩa thống kê. Chỉ diễn giải các biến dự báo có giá trị p trong kết quả tóm tắt thấp hơn ngưỡng bạn đã chọn, thường là 0.05.

# Convert log-odds coefficients into odds ratios
odds_ratio <- exp(coef(logit))
round(odds_ratio, 3)

# Odds ratios with 95% confidence intervals
exp(cbind(OddsRatio = coef(logit), confint(logit)))

Đọc đầu ra. Hệ số của biến hours.per.week trong mô hình của chúng ta là 0.41942. Lũy thừa của nó ta được exp(0.41942) = 1.52, có nghĩa là việc tăng một độ lệch chuẩn về số giờ làm việc hàng tuần sẽ làm tăng khả năng kiếm được trên 50 đô la lên khoảng 1.5 lần, trong khi giữ nguyên tất cả các biến số khác.

Các hệ số âm hoạt động theo cách tương tự. marital.statusNot_married là -2.50346, vì vậy exp(-2.50346) = 0.08: những người chưa kết hôn có tỷ lệ khả năng xảy ra sự kiện chỉ bằng khoảng 8% so với người đã kết hôn. Vì các biến dự báo liên tục đã được chuẩn hóa ở Bước 1, hãy mô tả sự thay đổi theo đơn vị độ lệch chuẩn, chứ không phải theo số giờ thực tế.

Lưu ý dành cho các gia đình khác: Các hệ số lũy thừa chỉ là tỷ lệ chênh lệch trong họ phân phối nhị thức với liên kết logit. Với họ phân phối = “poisson” và liên kết logarit, các giá trị exp() tương tự được đọc là tỷ lệ thay vì tỷ lệ chênh lệch.

GLM trong R: Tham khảo nhanh các hàm

Hãy để bảng này bên cạnh bạn trong khi lập trình. Bảng này liệt kê mọi hàm được sử dụng trong tám bước ở trên, cùng với gói cung cấp hàm đó và các đối số mà hàm mong đợi.

Cai khac GLM Các gia đình có thể được lựa chọn thông qua lập luận về gia đình là:

• nhị thức: (liên kết = “logit”)
• Gaussian: (liên kết = “đồng nhất”)
• Gamma: (liên kết = “ngược”)
• inverse.gaussian: (link = “1/mu^2”)
• poisson: (link = “log”)
• quasi: (liên kết = “đồng nhất”, phương sai = “hằng số”)
• quasibinomial: (link = “logit”)
• quasipoisson: (liên kết = “log”)

Câu Hỏi Thường Gặp