GLM di R: Model Linier Umum dengan Contoh

Apa itu regresi logistik?

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 diukur dengan fungsi hubungan probabilistik berikut yang disebut sigmoid.dll karena bentuknya 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

Cara membuat Generalized Liner Model (GLM)

Mari gunakan dewasa kumpulan data untuk menggambarkan regresi logistik. "Dewasa" adalah kumpulan data yang bagus untuk tugas klasifikasi. Tujuannya adalah untuk memprediksi apakah pendapatan tahunan dalam dolar seseorang akan melebihi 50.000. Kumpulan data berisi 46,033 observasi dan sepuluh fitur:

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

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)

Penjelasan Kode

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

• 1: Gambarkan distribusi jam per minggu
• 2: 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 2 persen teratas

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

Penjelasan Kode

• quantile(data_adult$hours.per.week, .99): Hitung nilai 99 persen waktu kerja

Keluaran:

## 99% 
##  80

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

Penjelasan Kode

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

Penjelasan Kode

• 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 terampil. Anda ingin memplot diagram batang untuk setiap kolom dalam faktor bingkai data. Akan lebih mudah untuk mengotomatisasi proses, terutama jika 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)))

Penjelasan Kode

• 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

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:

Tingkat lama	Level baru
Prasekolah	keluar
10th	Keluar
11th	Keluar
12th	Keluar
1st-4	Keluar
5th-6th	Keluar
7th-8th	Keluar
9th	Keluar
Lulusan HS	Lulusan Tinggi
Beberapa perguruan tinggi	Komunitas
Assoc-acdm	Komunitas
Assoc-voc	Komunitas
Sarjana	Sarjana
Masters	Masters
Sekolah Prof	Masters
Gelar doktor	PhD

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

Penjelasan Kode

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

Tingkat lama	Level baru
Tidak pernah menikah	Belum nikah
Menikah-pasangan-absen	Belum nikah
Menikah dengan pasangan AF	Menikah
Menikah-Civ-pasangan
Terpisah	Terpisah
Bercerai
Janda	Janda

# 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. Distribusinya memiliki banyak pilihan berbeda. Hal ini mungkin dapat dijelaskan oleh jenis kontrak di 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()

Penjelasan Kode

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

Penjelasan Kode

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

Penjelasan Kode

• 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 memerlukan untuk membagi data antara set kereta dan set pengujian. Anda dapat menggunakan "fungsi" yang Anda buat di tutorial pembelajaran terawasi lainnya untuk membuat set latihan/ujian.

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 paket glm(). Itu 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)

Penjelasan Kode

• 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 diekstraksi dengan tanda $ diikuti dengan nama metriknya. Misalnya, Anda menyimpan model sebagai logit. Untuk mengekstrak kriteria AIC, Anda menggunakan:

logit$aic

Keluaran:

## [1] 27086.65

Langkah 7) Menilai kinerja model

Matriks Kebingungan

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

Penjelasan Kode

• 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 konfusi mewakili target sebenarnya, sedangkan setiap kolom mewakili target prediksi. Baris pertama matriks ini menganggap pendapatan lebih rendah dari 50k (kelas Salah): 6241 diklasifikasikan dengan benar sebagai individu dengan pendapatan lebih rendah dari 50k (Benar-benar negatif), sedangkan sisanya salah diklasifikasikan sebagai di atas 50k (Salah positif). Baris kedua menganggap penghasilan di atas 50k, kelas positifnya 1229 (Benar-benar positif), selagi Benar-benar negatif adalah 1074.

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

Penjelasan Kode

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

Keluaran:

## [1] 0.8277339

Model tersebut tampaknya mengalami satu masalah, yaitu melebih-lebihkan jumlah negatif palsu. Ini disebut paradoks tes akurasi. Kami menyatakan bahwa akurasi adalah rasio prediksi yang benar terhadap jumlah kasus. Kita dapat memiliki akurasi yang relatif tinggi tetapi modelnya tidak berguna. Hal ini terjadi ketika ada kelas dominan. Jika Anda melihat kembali matriks konfusi, Anda dapat melihat sebagian besar kasus diklasifikasikan sebagai benar-benar negatif. Bayangkan sekarang, model tersebut mengklasifikasikan semua kelas sebagai negatif (yaitu lebih rendah dari 50k). Anda akan mendapatkan akurasi 75 persen (6718/6718+2257). Model Anda berperforma lebih baik tetapi kesulitan membedakan hal positif dan negatif sebenarnya.

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

Penjelasan Kode

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

Penjelasan Kode

• 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

Ketika model menyatakan bahwa itu adalah individu dengan berat badan di atas 50 ribu, maka model tersebut benar hanya pada 54 persen kasus, dan dapat mengklaim individu dengan berat badan di atas 50 ribu pada 72 persen kasus.

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

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 memplot kurva ROC, kita perlu menginstal perpustakaan bernama RORC. Kita bisa menemukannya di conda perpustakaan. Anda dapat mengetikkan kode:

conda install -cr r-rocr –ya

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

Penjelasan Kode

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

Anda perlu menggunakan tes skor untuk membandingkan kedua model

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

Skornya sedikit lebih tinggi dari yang sebelumnya. Anda dapat terus mengerjakan data untuk mencoba mengalahkan skor.

Kesimpulan

Kita dapat meringkas fungsi untuk melatih regresi logistik dalam tabel di bawah ini:

Paket	Tujuan	fungsi	Argumen
-	Buat kumpulan data pelatihan/pengujian	buat_train_set()	data, ukuran, kereta
glm	Melatih Model Linear Umum	glm()	rumus, data, keluarga*
glm	Rangkum modelnya	ringkasan()	model terpasang
mendasarkan	Buat prediksi	meramalkan()	model yang dipasang, kumpulan data, tipe = 'respons'
mendasarkan	Buat matriks kebingungan	meja()	y, prediksi()
mendasarkan	Buat skor akurasi	jumlah(diag(tabel())/jumlah(tabel()
ROCR	Buat ROC : Langkah 1 Buat prediksi	ramalan()	prediksi(), y
ROCR	Buat ROC : Langkah 2 Buat kinerja	pertunjukan()	prediksi(), 'tpr', 'fpr'
ROCR	Buat ROC : Langkah 3 Plot grafik	merencanakan()	pertunjukan()

Yang lain GLM jenis modelnya adalah:

– binomial: (tautan = “logit”)

– gaussian: (tautan = “identitas”)

– Gamma: (tautan = “terbalik”)

– invers.gaussian: (link = “1/mu^2”)

– racun: (tautan = “log”)

– kuasi: (tautan = “identitas”, varians = “konstan”)

– kuasibinomial: (tautan = “logit”)

– kuasipoisson: (tautan = “log”)