GLM di R: Model Linier Umum & Regresi Logistik

Apa itu Model Linier Umum (Generalized Linear Model/GLM) di R?

A Model Linier Umum (GLM) memperluas regresi linier biasa sehingga variabel respons dapat mengikuti distribusi selain distribusi normal. Dalam RAnda memasangnya dengan yang terintegrasi. glm() fungsi dari paket statistik.

Setiap GLM didefinisikan oleh tiga komponen:

• Komponen acak: Distribusi probabilitas variabel respons, diambil dari keluarga eksponensial (binomial, Poisson, Gamma, Gaussian, dan lainnya).
• Komponen sistematis: Prediktor linier, yaitu kombinasi tertimbang dari variabel penjelas Anda.
• Fungsi tautan: Fungsi yang menghubungkan nilai rata-rata respons dengan prediktor linier, misalnya tautan logit untuk data biner atau tautan log untuk data hitungan.

Struktur itulah yang membuat model tersebut "tergeneralisasi". Alih-alih memaksakan hasil yang terdistribusi normal, Anda menyatakan distribusi yang benar melalui keluarga Argumen dan R memperkirakan koefisien dengan metode kemungkinan maksimum. Regresi logistik hanyalah GLM dengan keluarga binomial dan tautan logit, jadi ini adalah titik awal yang alami.

Apa itu Regresi Logistik di R?

Regresi logistik digunakan untuk memprediksi suatu kelas, yaitu probabilitas. Regresi logistik dapat memprediksi hasil biner secara akurat.

Bayangkan Anda ingin memprediksi apakah suatu pinjaman ditolak/diterima berdasarkan banyak atribut. Regresi logistik berbentuk 0/1. y = 0 jika pinjaman ditolak, y = 1 jika diterima.

Model regresi logistik berbeda dari model regresi linier dalam dua hal.

• Pertama-tama, regresi logistik hanya menerima masukan dikotomis (biner) sebagai variabel terikat (yaitu vektor 0 dan 1).
• Kedua, hasilnya dipetakan melalui fungsi penghubung probabilistik yang disebut sigmoid.dll (fungsi logistik) karena bentuknya yang menyerupai huruf S:

Output dari fungsi selalu antara 0 dan 1. Periksa Gambar di bawah

Fungsi sigmoid mengembalikan nilai dari 0 hingga 1. Untuk tugas klasifikasi, kita memerlukan keluaran diskrit 0 atau 1.

Untuk mengubah aliran kontinu menjadi nilai diskrit, kita dapat menetapkan batas keputusan pada 0.5. Semua nilai di atas ambang batas ini diklasifikasikan sebagai 1

Setelah fungsi penghubung dipahami, bandingkan model umum dengan model linier biasa yang sudah Anda ketahui.

GLM vs Regresi Linier: Perbedaan Utama dalam R

Sebelum Anda menulis kode apa pun, ada baiknya mengetahui dengan tepat kapan standar tersebut berlaku. regresi linier Fungsi lm() sudah tidak sesuai lagi dan glm() harus menggantikannya.

Singkatnya, pilih lm() ketika hasilnya adalah pengukuran yang terdistribusi normal dan glm() setiap kali hasilnya adalah keputusan ya/tidak, tugas yang biasanya dapat Anda serahkan kepada seorang pohon keputusan pengklasifikasi, penghitung kejadian, atau kuantitas yang bernilai positif yang sebarannya meningkat seiring dengan nilai rata-ratanya.

Cara Membuat Model Linear Umum (GLM) di R

Setelah teori dipahami, sisa tutorial ini menerapkan GLM binomial secara menyeluruh pada dataset nyata.

Mari gunakan dewasa Kumpulan data ini digunakan untuk mengilustrasikan regresi logistik. Dataset "dewasa" sangat cocok untuk tugas klasifikasi. Tujuannya adalah untuk memprediksi apakah pendapatan tahunan seseorang dalam dolar AS akan melebihi 50,000. Dataset ini berisi 48,842 observasi dan sepuluh variabel:

• usia: usia individu. numerik
• pendidikan: Tingkat pendidikan individu. Faktor.
• status perkawinan: Maristatus tinggi individu. Faktor yaitu Belum Menikah, Menikah dengan WNI, …
• gender : Jenis kelamin individu. Faktor yaitu Laki-laki atau Perempuan
• penghasilan: Target variabel. Penghasilan di atas atau di bawah 50K. Faktor yaitu >50K, <=50K

di antara yang lain

library(dplyr)
data_adult <-read.csv("https://raw.githubusercontent.com/guru99-edu/R-Programming/master/adult.csv")
glimpse(data_adult)

Keluaran:

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

Kami akan melanjutkan sebagai berikut:

• Langkah 1: Periksa variabel kontinu
• Langkah 2: Periksa variabel faktor
• Langkah 3: Rekayasa fitur
• Langkah 4: Ringkasan statistik
• Langkah 5: Latihan/set pengujian
• Langkah 6: Bangun modelnya
• Langkah 7: Menilai kinerja model
• Langkah 8: Tingkatkan modelnya

Tugas Anda adalah memprediksi individu mana yang akan memiliki pendapatan lebih dari 50K.

Dalam tutorial ini, setiap langkah akan dirinci untuk melakukan analisis pada dataset nyata.

Langkah 1) Periksa variabel kontinu

Pada langkah pertama, Anda dapat melihat distribusi variabel kontinu.

continuous <-select_if(data_adult, is.numeric)
summary(continuous)

Code Penjelasan

• continuous <- select_if(data_adult, is.numeric): Gunakan fungsi select_if() dari pustaka dplyr untuk memilih kolom numerik saja
• ringkasan (berkelanjutan): Cetak statistik ringkasan

Keluaran:

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

Dari tabel di atas, Anda dapat melihat bahwa data memiliki skala yang benar-benar berbeda dan jam per minggu memiliki outlier yang besar (yaitu lihat kuartil terakhir dan nilai maksimum).

Anda dapat mengatasinya dengan mengikuti dua langkah berikut:

• Gambarkan distribusi jam per minggu.
• Standarisasi variabel kontinu

1. Plot distribusinya

Mari kita lihat lebih dekat distribusi jam per minggu

# Histogram with kernel density curve
library(ggplot2)
ggplot(continuous, aes(x = hours.per.week)) +
    geom_density(alpha = .2, fill = "#FF6666")

Keluaran:

Variabel tersebut memiliki banyak outlier dan distribusinya tidak terdefinisi dengan baik. Anda dapat mengatasi sebagian masalah ini dengan menghapus 0.01 persen jam teratas per minggu.

Sintaks dasar kuantil:

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.

Kami menghitung persentil ke-99 dari jam kerja mingguan.

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

Code Penjelasan

• quantile(data_adult$hours.per.week, .99): Menghitung persentil ke-99 dari waktu kerja mingguan

Keluaran:

## 99% 
##  80

99 persen penduduk bekerja kurang dari 80 jam per minggu.

Anda dapat menurunkan observasi di atas ambang batas ini. Anda menggunakan filter dari dplyr Perpustakaan.

data_adult_drop <-data_adult %>%
filter(hours.per.week<top_one_percent)
dim(data_adult_drop)

Keluaran:

## [1] 45537    10

2. Standarisasi variabel kontinu

Anda dapat menstandardisasi setiap kolom untuk meningkatkan performa karena data Anda tidak memiliki skala yang sama. Anda dapat menggunakan fungsi mutate_if dari perpustakaan dplyr. Sintaks dasarnya adalah:

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

Anda dapat membakukan kolom numerik sebagai berikut:

data_adult_rescale <- data_adult_drop %>%
	mutate_if(is.numeric, funs(as.numeric(scale(.))))
head(data_adult_rescale)

Code Penjelasan

• mutate_if(is.numeric, funs(scale)): Syaratnya hanya kolom numerik dan fungsinya skala

Keluaran:

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

Langkah 2) Periksa variabel faktor

Langkah ini memiliki dua tujuan:

• Periksa level di setiap kolom kategori
• Tentukan level baru

Kami akan membagi langkah ini menjadi tiga bagian:

• Pilih kolom kategorikal
• Simpan diagram batang setiap kolom dalam daftar
• Cetak grafiknya

Kita dapat memilih kolom faktor dengan kode di bawah ini:

# Select categorical column
factor <- data.frame(select_if(data_adult_rescale, is.factor))
	ncol(factor)

Code Penjelasan

• data.frame(select_if(data_adult, is.factor)): Kami menyimpan kolom faktor dalam faktor dalam tipe bingkai data. Perpustakaan ggplot2 memerlukan objek bingkai data.

Keluaran:

## [1] 6

Dataset berisi 6 variabel kategori

Langkah kedua lebih menantang. Anda ingin membuat grafik batang untuk setiap kolom dalam data frame faktor. Akan lebih mudah untuk mengotomatiskan proses ini, terutama dalam situasi di mana terdapat banyak kolom.

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 Penjelasan

• lapply(): Gunakan fungsi lapply() untuk meneruskan fungsi di semua kolom kumpulan data. Anda menyimpan hasilnya dalam daftar
• function(x): Fungsi akan diproses untuk setiap x. Di sini x adalah kolomnya
• ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): Buat diagram batang char untuk setiap elemen x. Catatan, untuk mengembalikan x sebagai kolom, Anda perlu memasukkannya ke dalam get()

Langkah terakhir relatif mudah. Anda ingin mencetak 6 grafik.

# Print the graph
graph

Keluaran:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Catatan: Gunakan tombol berikutnya untuk menavigasi ke grafik berikutnya

Langkah 3) Rekayasa fitur

Dua variabel kategorikal memiliki lebih banyak level daripada yang dibutuhkan model. Anda akan mengelompokkannya kembali ke dalam kategori yang lebih luas dan lebih lengkap.

Merombak pendidikan

Dari grafik di atas terlihat bahwa variabel pendidikan mempunyai 16 jenjang. Hal ini penting, dan beberapa tingkat memiliki jumlah observasi yang relatif rendah. Jika Anda ingin meningkatkan jumlah informasi yang dapat Anda peroleh dari variabel ini, Anda dapat menyusunnya kembali ke tingkat yang lebih tinggi. Yaitu, Anda membuat grup yang lebih besar dengan tingkat pendidikan yang sama. Misalnya, tingkat pendidikan yang rendah akan dikonversi menjadi putus sekolah. Jenjang pendidikan yang lebih tinggi akan diubah menjadi magister.

Berikut detailnya:

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 Penjelasan

• Kami menggunakan kata kerja bermutasi dari perpustakaan dplyr. Nilai-nilai pendidikan kita ubah dengan pernyataan ifelse

Pada tabel di bawah, Anda membuat ringkasan statistik untuk melihat, rata-rata, berapa tahun pendidikan (nilai z) yang diperlukan untuk mencapai gelar Sarjana, Magister, atau PhD.

recast_data %>%
	group_by(education) %>%
	summarize(average_educ_year = mean(educational.num),
		count = n()) %>%
	arrange(average_educ_year)

Keluaran:

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

Perombakan Maristatus tal

Anda juga dapat membuat level yang lebih rendah untuk status perkawinan. Dalam kode berikut, Anda mengubah level sebagai berikut:

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

Anda dapat memeriksa jumlah individu dalam setiap grup.

table(recast_data$marital.status)

Keluaran:

## ##     Married Not_married   Separated       Widow
##       21165       15359        7727        1286

Langkah 4) Ringkasan Statistik

Sekarang saatnya untuk memeriksa beberapa statistik tentang variabel target kita. Pada grafik di bawah, Anda menghitung persentase individu yang berpenghasilan lebih dari 50 ribu berdasarkan jenis kelamin mereka.

# Plot gender income
ggplot(recast_data, aes(x = gender, fill = income)) +
    geom_bar(position = "fill") +
    theme_classic()

Keluaran:

Selanjutnya, periksa apakah asal usul individu memengaruhi penghasilannya.

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

Keluaran:

Jumlah jam kerja berdasarkan jenis kelamin.

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

Keluaran:

Diagram kotak menegaskan bahwa distribusi waktu kerja sesuai dengan kelompok yang berbeda. Dalam diagram kotak, kedua jenis kelamin tidak memiliki pengamatan yang homogen.

Anda dapat memeriksa kepadatan waktu kerja mingguan berdasarkan jenis pendidikan. Distribusi tersebut memiliki banyak puncak yang berbeda. Hal ini mungkin dapat dijelaskan oleh jenis pendidikan.tracdi AS.

# 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 Penjelasan

• ggplot(recast_data, aes( x= jam.per.minggu)): Plot kepadatan hanya memerlukan satu variabel
• geom_density(aes(color = education), alpha =0.5): Objek geometris untuk mengontrol kepadatan

Keluaran:

Untuk mengkonfirmasi pemikiran Anda, Anda dapat melakukan satu arah Uji ANOVA:

anova <- aov(hours.per.week~education, recast_data)
summary(anova)

Keluaran:

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

Uji ANOVA mengonfirmasi perbedaan rata-rata antar kelompok.

Non-linearitas

Sebelum menjalankan model, Anda dapat melihat apakah jumlah jam kerja berhubungan dengan usia.

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 Penjelasan

• ggplot(recast_data, aes(x = age, y = hours.per.week)): Mengatur estetika grafik
• geom_point(aes(color= pendapatan), size =0.5): Buatlah plot titik
• stat_smooth(): Tambahkan garis tren dengan argumen berikut:
  • metode='lm': Plot nilai yang dipasang jika regresi linier
  • rumus = y~poli(x,2): Cocok dengan regresi polinomial
  • se = TRUE: Tambahkan kesalahan standar
  • aes(color= income): Pecahkan model berdasarkan pendapatan

Keluaran:

Singkatnya, Anda dapat menguji istilah interaksi dalam model untuk mengetahui efek non-linearitas antara waktu kerja mingguan dan fitur lainnya. Penting untuk mendeteksi dalam kondisi apa waktu kerja berbeda.

Korelasi

Pemeriksaan selanjutnya adalah memvisualisasikan korelasi antar variabel. Anda mengonversi tipe tingkat faktor menjadi numerik sehingga Anda dapat memplot peta panas yang berisi koefisien korelasi yang dihitung dengan metode 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 Penjelasan

• data.frame(lapply(recast_data,as.integer)): Mengonversi data menjadi numerik
• ggcorr() memplot peta panas dengan argumen berikut:
  • metode: Metode untuk menghitung korelasi
  • nbreaks = 6: Jumlah istirahat
  • hjust = 0.8: Mengontrol posisi nama variabel di plot
  • label = TRUE: Tambahkan label di tengah jendela
  • label_size = 3: Label ukuran
  • color = “grey50”): Warna label

Keluaran:

Langkah 5) Set pelatihan/pengujian

Setiap diawasi Mesin belajar Tugas ini mengharuskan Anda untuk membagi data antara set pelatihan dan set pengujian. Anda dapat menggunakan "fungsi" yang telah Anda buat di tutorial pembelajaran terawasi lainnya untuk membuat set pelatihan/pengujian.

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)

Keluaran:

## [1] 36429     9

dim(data_test)

Keluaran:

## [1] 9108    9

Langkah 6) Bangun modelnya

Untuk melihat bagaimana kinerja algoritma, Anda menggunakan fungsi glm() dari paket stats. Model Linier Umum adalah kumpulan model. Sintaks dasarnya adalah:

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

Anda siap memperkirakan model logistik untuk membagi tingkat pendapatan antar serangkaian fitur.

formula <- income~.
logit <- glm(formula, data = data_train, family = 'binomial')
summary(logit)

Code Penjelasan

• rumus <- pendapatan ~ .: Buat model yang sesuai
• logit <- glm(rumus, data = data_train, family = 'binomial'): Cocokkan model logistik (family = 'binomial') dengan data data_train.
• ringkasan(logit): Cetak ringkasan model

Keluaran:

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

Ringkasan model kami mengungkapkan informasi menarik. Kinerja regresi logistik dievaluasi dengan metrik utama tertentu.

• AIC (Kriteria Informasi Akaike): Ini setara dengan R2 dalam regresi logistik. Ini mengukur kecocokan ketika penalti diterapkan pada sejumlah parameter. Lebih kecil AIC nilai menunjukkan model tersebut mendekati kebenaran.
• Penyimpangan nol: Cocok dengan model hanya dengan intersep. Derajat kebebasannya adalah n-1. Kita dapat mengartikannya sebagai nilai Chi-square (nilai fit yang berbeda dengan nilai sebenarnya pengujian hipotesis).
• Residual Deviance: Model dengan semua variabel. Hal ini juga diartikan sebagai pengujian hipotesis Chi-square.
• Jumlah iterasi Fisher Scoring: Jumlah iterasi sebelum konvergen.

Output dari fungsi glm() disimpan dalam daftar. Kode di bawah ini menunjukkan semua item yang tersedia dalam variabel logit yang kami buat untuk mengevaluasi regresi logistik.

# Daftarnya sangat panjang, cetak hanya tiga elemen pertama

lapply(logit, class)[1:3]

Keluaran:

## $coefficients
## [1] "numeric"
## 
## $residuals
## [1] "numeric"
## 
## $fitted.values
## [1] "numeric"

Setiap nilai dapat berupatracdiawali dengan tanda $ diikuti nama metrik. Misalnya, Anda menyimpan model sebagai logit. Untuk contohtracBerdasarkan kriteria AIC, Anda menggunakan:

logit$aic

Keluaran:

## [1] 27086.65

Langkah 7) Menilai kinerja model

Matriks Kebingungan

The matriks kebingungan adalah pilihan yang lebih baik untuk mengevaluasi kinerja klasifikasi dibandingkan dengan berbagai metrik yang Anda lihat sebelumnya. Ide umumnya adalah menghitung berapa kali instance True diklasifikasikan menjadi False.

Untuk menghitung matriks konfusi, pertama-tama Anda harus memiliki serangkaian prediksi agar dapat dibandingkan dengan target sebenarnya.

predict <- predict(logit, data_test, type = 'response')
# confusion matrix
table_mat <- table(data_test$income, predict > 0.5)
table_mat

Code Penjelasan

• prediksi(logit,data_test, type = 'response'): Hitung prediksi pada set pengujian. Tetapkan type = 'response' untuk menghitung probabilitas respons.
• table(data_test$income, prediksi > 0.5): Hitung matriks kebingungan. prediksi > 0.5 berarti ia mengembalikan 1 jika probabilitas prediksi di atas 0.5, jika tidak 0.

Keluaran:

##        
##         FALSE TRUE
##   <=50K  6310  495
##   >50K   1074 1229	

Setiap baris dalam matriks kebingungan mewakili target sebenarnya, sedangkan setiap kolom mewakili target yang diprediksi. Baris pertama matriks ini mempertimbangkan pendapatan di bawah 50 ribu (kelas negatif): 6,310 observasi diklasifikasikan dengan benar sebagai individu dengan pendapatan di bawah 50 ribu (Benar-benar negatif), sementara 495 salah diklasifikasikan sebagai di atas 50 ribu (Salah positif). Baris kedua mempertimbangkan pendapatan di atas 50 ribu: 1,229 diidentifikasi dengan benar (Benar-benar positif), sedangkan 1,074 tidak terdeteksi (Salah negatif).

Anda dapat menghitung modelnya ketepatan dengan menjumlahkan positif benar + negatif benar atas total observasi

accuracy_Test <- sum(diag(table_mat)) / sum(table_mat)
accuracy_Test

Code Penjelasan

• sum(diag(table_mat)): Jumlah diagonalnya
• sum(table_mat): Jumlah matriks.

Keluaran:

## [1] 0.8277339

Model ini tampaknya memiliki satu masalah: menghasilkan terlalu banyak false negative. Ini disebut paradoks tes akurasiKami menyatakan bahwa akurasi adalah rasio prediksi yang benar terhadap jumlah total kasus. Kita dapat memiliki akurasi yang relatif tinggi tetapi model yang tidak berguna. Hal ini terjadi ketika ada kelas dominan. Jika Anda melihat kembali matriks kebingungan, Anda dapat melihat sebagian besar kasus diklasifikasikan sebagai negatif sejati. Bayangkan sekarang, model mengklasifikasikan setiap observasi sebagai negatif (yaitu kurang dari 50 ribu). Anda masih akan memperoleh akurasi sekitar 75 persen (6,805 / 9,108). Model Anda berkinerja lebih baik tetapi kesulitan membedakan positif sejati dengan negatif sejati.

Dalam situasi seperti ini, sebaiknya memiliki metrik yang lebih ringkas. Kita bisa melihat:

• Presisi=TP/(TP+FP)
• Penarikan=TP/(TP+FN)

Presisi vs Ingatan

Ketelitian melihat keakuratan prediksi positif. Mengingat kembali adalah rasio kejadian positif yang terdeteksi dengan benar oleh pengklasifikasi;

Anda dapat membuat dua fungsi untuk menghitung dua metrik ini

1. Membangun presisi

precision <- function(matrix) {
	# True positive
    tp <- matrix[2, 2]
	# false positive
    fp <- matrix[1, 2]
    return (tp / (tp + fp))
}

Code Penjelasan

• mat[1,1]: Mengembalikan sel pertama dari kolom pertama bingkai data, yaitu positif sebenarnya
• tikar[1,2]; Mengembalikan sel pertama dari kolom kedua bingkai data, yaitu positif palsu

recall <- function(matrix) {
# true positive
    tp <- matrix[2, 2]# false positive
    fn <- matrix[2, 1]
    return (tp / (tp + fn))
}

Code Penjelasan

• mat[1,1]: Mengembalikan sel pertama dari kolom pertama bingkai data, yaitu positif sebenarnya
• tikar[2,1]; Mengembalikan sel kedua dari kolom pertama bingkai data, yaitu negatif palsu

Anda dapat menguji fungsi Anda

prec <- precision(table_mat)
prec
rec <- recall(table_mat)
rec

Keluaran:

## [1] 0.712877
## [2] 0.5336518

Bacalah kedua angka ini dengan saksama. Presisi adalah 0.71, jadi ketika model mengatakan seseorang berpenghasilan di atas 50 ribu, model tersebut benar dalam 71 persen kasus. Recall adalah 0.53, jadi model hanya mendeteksi 53 persen dari individu yang benar-benar berpenghasilan di atas 50 ribu.

Anda dapat membuat file skor berdasarkan presisi dan recall. Itu adalah rata-rata harmonis dari kedua metrik ini, yang berarti memberikan bobot lebih pada nilai yang lebih rendah.

f1 <- 2 * ((prec * rec) / (prec + rec))
f1

Keluaran:

## [1] 0.6103799

Pengorbanan Presisi vs Perolehan

Tidak mungkin memiliki presisi tinggi dan recall yang tinggi.

Jika kita meningkatkan presisi, individu yang tepat akan dapat diprediksi dengan lebih baik, namun kita akan kehilangan banyak individu (recall lebih rendah). Dalam beberapa situasi, kami lebih memilih presisi yang lebih tinggi daripada perolehan kembali. Ada hubungan cekung antara presisi dan perolehan.

• Bayangkan, Anda perlu memprediksi apakah seorang pasien mengidap suatu penyakit. Anda ingin setepat mungkin.
• Jika Anda perlu mendeteksi calon penipu di jalan melalui pengenalan wajah, akan lebih baik jika menangkap banyak orang yang dicap penipu meski presisinya rendah. Polisi akan dapat membebaskan individu yang tidak melakukan penipuan.

Kurva ROC

The Penerima OperaKarakteristik kurva adalah alat umum lainnya yang digunakan dengan klasifikasi biner. Hal ini sangat mirip dengan kurva presisi/recall, namun alih-alih memplot presisi versus perolehan, kurva ROC menunjukkan tingkat positif sebenarnya (yaitu, perolehan kembali) terhadap tingkat positif palsu. Tingkat positif palsu adalah rasio kejadian negatif yang salah diklasifikasikan sebagai positif. Nilai ini sama dengan satu dikurangi tingkat negatif sebenarnya. Tingkat negatif sebenarnya juga disebut kekhususan. Oleh karena itu plot kurva ROC kepekaan (ingat) versus 1 kekhususan

Untuk membuat grafik kurva ROC, kita perlu menginstal paket bernama ROCR. Kita dapat menemukannya di conda. perpustakaan. Anda dapat mengetikkan kode:

conda install -c r r-rocr --yes

Kita dapat memplot ROC dengan fungsi prediksi() dan kinerja().

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 Penjelasan

• prediksi(prediksi, data_test$pendapatan): Perpustakaan ROCR perlu membuat objek prediksi untuk mengubah data masukan
• performance(ROCRpred, 'tpr','fpr'): Mengembalikan dua kombinasi yang akan dihasilkan dalam grafik. Di sini, tpr dan fpr dibangun. Untuk memplot presisi dan mengingat secara bersamaan, gunakan “prec”, “rec”.

Keluaran:

Langkah 8) Tingkatkan modelnya

Anda dapat mencoba menambahkan non-linearitas pada model dengan interaksi antar keduanya

• usia dan jam per minggu
• jenis kelamin dan jam per minggu.

Kemudian, Anda membandingkan skor F1 dari kedua model tersebut.

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

Keluaran:

## [1] 0.6109181

Skor F1 sedikit lebih tinggi dari sebelumnya. Anda dapat terus mengolah data dan mencoba untuk mengalahkan skor tersebut.

Cara Menginterpretasikan Koefisien GLM dan Rasio Peluang di R

Tabel ringkasan yang dicetak pada Langkah 6 melaporkan koefisien pada log-odds skala, yang sulit dijelaskan kepada audiens non-teknis. Mengubahnya menjadi rasio peluang membuat model jauh lebih mudah dikomunikasikan.

Ikuti keempat langkah ini.

1. Pangkatkan koefisiennya. Terapkan fungsi exp() pada setiap estimasi sehingga log-odds menjadi rasio odds multiplikatif.
2. Tambahkan interval kepercayaan. Bungkus fungsi confint() dengan fungsi exp() untuk mendapatkan interval 95 persen pada skala peluang yang sama.
3. Bandingkan setiap nilai dengan 1. Rasio odds di atas 1 meningkatkan probabilitas kelas positif, nilai di bawah 1 menurunkannya, dan nilai mendekati 1 berarti prediktor tersebut memberikan sedikit tambahan.
4. Periksa signifikansi statistik. Hanya interpretasikan prediktor yang nilai p-nya dalam output ringkasan berada di bawah ambang batas yang Anda pilih, biasanya 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)))

Membaca hasilnya. Koefisien untuk jam per minggu dalam model kita adalah 0.41942. Dengan memangkatkannya, kita peroleh exp(0.41942) = 1.52, yang berarti bahwa peningkatan satu standar deviasi dalam jam kerja mingguan akan melipatgandakan peluang untuk mendapatkan penghasilan di atas 50 ribu sekitar 1.5 kali, dengan asumsi semua variabel lain tetap konstan.

Koefisien negatif bekerja dengan cara yang sama. marital.statusNot_married adalah -2.50346, jadi exp(-2.50346) = 0.08: individu yang belum menikah memiliki peluang sekitar 8 persen dibandingkan individu yang sudah menikah. Karena prediktor kontinu telah distandarisasi pada Langkah 1, jelaskan perubahan dalam satuan deviasi standar, bukan jam mentah.

Catatan untuk keluarga lain: Koefisien yang dipangkatkan hanya merupakan rasio peluang di bawah keluarga binomial dengan tautan logit. Dengan keluarga = “poisson” dan tautan log, nilai exp() yang sama dibaca sebagai rasio laju.

GLM di R: Referensi Fungsi Cepat

Simpan tabel ini di samping Anda saat Anda melakukan pengkodean. Tabel ini mencantumkan setiap fungsi yang digunakan dalam delapan langkah di atas, beserta paket yang menyediakannya dan argumen yang diharapkan.

Yang lain GLM Keluarga yang tersedia melalui proses pengajuan permohonan bantuan keluarga adalah:

• binomial: (tautan = “logit”)
• gaussian: (tautan = “identitas”)
• Gamma: (tautan = “kebalikan”)
• invers.gaussian: (tautan = “1/mu^2”)
• poisson: (tautan = “log”)
• quasi: (link = “identitas”, varians = “konstan”)
• kuasibinomial: (tautan = “logit”)
• quasipoisson: (tautan = “log”)

Pertanyaan Umum Demo Slot