GLM in R: 例を含む一般化線形モデル
エブリン・クラーク
ロジスティック回帰とは何ですか?
ロジスティック回帰は、クラス、つまり確率を予測するために使用されます。 ロジスティック回帰では、バイナリの結果を正確に予測できます。
多くの属性に基づいてローンが拒否されるか受け入れられるかを予測したいと想像してください。 ロジスティック回帰の形式は 0/1 です。 ローンが拒否された場合は y = 0、承認された場合は y = 1。
ロジスティック回帰モデルは、XNUMX つの点で線形回帰モデルとは異なります。
- まず第一に、ロジスティック回帰は従属変数として二値 (バイナリ) 入力のみを受け入れます (つまり、0 と 1 のベクトル)。
- 第二に、結果は次の確率リンク関数によって測定されます。 シグモイド S字型のため。:
関数の出力は常に 0 から 1 の間です。下の画像を確認してください。
シグモイド関数は 0 から 1 までの値を返します。分類タスクでは、0 または 1 の離散出力が必要です。
連続フローを離散値に変換するには、決定限界を 0.5 に設定します。 このしきい値を超えるすべての値は 1 として分類されます。
一般化線形モデル (GLM) の作成方法
を使ってみましょう 成人 ロジスティック回帰を説明するためのデータセット。 「成人」は分類タスクに最適なデータセットです。 目的は、個人の年収がドルで 50.000 ドルを超えるかどうかを予測することです。 データセットには 46,033 個の観測値と XNUMX 個の特徴が含まれています。
- age: 個人の年齢。 数値
- 教育: 個人の教育レベル。 要素。
- 配偶者の有無： Mari個人の全体的なステータス。要素すなわち未婚、既婚配偶者、…
- 性別: 個人の性別。 要素、つまり男性または女性
- 所得： Target 変数。収入が 50K 以上か以下か。係数、つまり >50K、<=50K
とりわけ
library(dplyr)
data_adult <-read.csv("https://raw.githubusercontent.com/guru99-edu/R-Programming/master/adult.csv")
glimpse(data_adult)
出力：
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...
以下のように進めていきます。
- ステップ 1: 連続変数を確認する
- ステップ 2: 因子変数を確認する
- ステップ 3: 特徴量エンジニアリング
- ステップ 4: 統計の要約
- ステップ 5: トレーニング/テスト セット
- ステップ 6: モデルを構築する
- ステップ 7: モデルのパフォーマンスを評価する
- ステップ 8: モデルを改善する
あなたの仕事は、どの個人の収益が 50 を超えるかを予測することです。
このチュートリアルでは、実際のデータセットで分析を実行するための各ステップについて詳しく説明します。
ステップ 1) 連続変数を確認する
最初のステップでは、連続変数の分布を確認できます。
continuous <-select_if(data_adult, is.numeric) summary(continuous)
コードの説明
- Continuous <- select_if(data_social, is.numeric): dplyr ライブラリの関数 select_if() を使用して、数値列のみを選択します
- summary(continuous): 統計の概要を出力します。
出力：
## 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
上記の表から、データのスケールがまったく異なり、hours.per.weeks に大きな外れ値があることがわかります (つまり、最後の四分位値と最大値を見てください)。
次の 2 つの手順で対処できます。
- 1: 週あたりの労働時間の分布をプロットする
- 2: 連続変数を標準化する
- 分布をプロットする
週あたりの労働時間の分布を詳しく見てみましょう
# Histogram with kernel density curve
library(ggplot2)
ggplot(continuous, aes(x = hours.per.week)) +
geom_density(alpha = .2, fill = "#FF6666")
出力：
変数には外れ値が多く、分布が明確に定義されていません。この問題は、週あたりの時間の上位 0.01 パーセントを削除することで部分的に解決できます。
分位数の基本構文:
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.
上位 2% のパーセンタイルを計算します
top_one_percent <- quantile(data_adult$hours.per.week, .99) top_one_percent
コードの説明
- quantile(data_adult$hours.per.week, .99): 労働時間の99パーセントの値を計算する
出力：
## 99% ## 80
人口の98パーセントは週80時間未満で働いています。
このしきい値を超える観測値は削除できます。 のフィルターを使用します。 dplyr としょうかん。
data_adult_drop <-data_adult %>% filter(hours.per.week<top_one_percent) dim(data_adult_drop)
出力：
## [1] 45537 10
- 連続変数を標準化する
データのスケールが同じではないため、各列を標準化してパフォーマンスを向上させることができます。 dplyr ライブラリの関数 mutate_if を使用できます。 基本的な構文は次のとおりです。
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
次のように数値列を標準化できます。
data_adult_rescale <- data_adult_drop % > % mutate_if(is.numeric, funs(as.numeric(scale(.)))) head(data_adult_rescale)
コードの説明
- mutate_if(is.numeric, funs(scale)): 条件は数値列のみで、関数はscaleです
出力：
## 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
ステップ 2) 因子変数を確認する
このステップには XNUMX つの目的があります。
- 各カテゴリ列のレベルを確認する
- 新しいレベルを定義する
このステップは XNUMX つの部分に分けて説明します。
- カテゴリ列を選択します
- 各列の棒グラフをリストに保存する
- グラフを印刷する
以下のコードを使用して因子列を選択できます。
# Select categorical column factor <- data.frame(select_if(data_adult_rescale, is.factor)) ncol(factor)
コードの説明
- data.frame(select_if(data_social, is.factor)): データ フレーム タイプの要因に要因列を格納します。 ライブラリ ggplot2 にはデータ フレーム オブジェクトが必要です。
出力：
## [1] 6
データセットには 6 つのカテゴリ変数が含まれています
XNUMX 番目のステップはより熟練したものになります。 データ フレーム係数の列ごとに棒グラフをプロットするとします。 特に列がたくさんある状況では、プロセスを自動化する方が便利です。
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)))
コードの説明
- lagply(): 関数lapply()を使用して、データセットのすべての列に関数を渡します。 出力をリストに保存します
- function(x): 関数は x ごとに処理されます。 ここで x は列です
- ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): x 要素ごとに棒グラフを作成します。 x を列として返すには、それを get() 内に含める必要があることに注意してください。
最後のステップは比較的簡単です。 6 つのグラフを印刷したいとします。
# Print the graph graph
出力：
## [[1]]
## ## [[2]]
## ## [[3]]
## ## [[4]]
## ## [[5]]
## ## [[6]]
注: 次のグラフに移動するには、「次へ」ボタンを使用します。
ステップ 3) 特徴量エンジニアリング
リキャスト教育
上のグラフから、可変教育には 16 のレベルがあることがわかります。 これはかなりの量であり、一部のレベルでは観測値の数が比較的少なくなります。 この変数から取得できる情報量を増やしたい場合は、変数をより高いレベルに再キャストできます。 つまり、同じレベルの教育を受けたより大きなグループを作成します。 たとえば、教育レベルが低いと中退に変換されます。 高等教育レベルは修士に変更されます。
詳細は次のとおりです。
| 古いレベル | 新しいレベル |
|---|---|
| プレスクール | 脱落 |
| 10 | ドロップアウト |
| 11 | ドロップアウト |
| 12 | ドロップアウト |
| 1日〜4日 | ドロップアウト |
| 5th-6th | ドロップアウト |
| 7th-8th | ドロップアウト |
| 9 | ドロップアウト |
| 高校卒業 | 高等教育 |
| 某大学 | コミュニティ |
| Assoc-acdm | コミュニティ |
| Assoc-voc | コミュニティ |
| 学士 | 学士 |
| マスターズ | マスターズ |
| 教授学校 | マスターズ |
| 博士 | 博士 |
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")))))))
コードの説明
- dplyr ライブラリの動詞 mutate を使用します。 ifelse ステートメントで教育の価値観を変える
以下の表では、学士号、修士号、または博士号に到達するまでに平均して何年教育が必要か (Z 値) を確認するための要約統計を作成します。
recast_data % > % group_by(education) % > % summarize(average_educ_year = mean(educational.num), count = n()) % > % arrange(average_educ_year)
出力：
## # 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
リキャスト Mariタルステータス
婚姻状況の下位レベルを作成することも可能です。次のコードでは、レベルを次のように変更します。
| 古いレベル | 新しいレベル |
|---|---|
| 結婚したことがない | 結婚していない |
| 既婚配偶者不在 | 結婚していない |
| 既婚AF配偶者 | 既婚 |
| 既婚-市民-配偶者 | |
| 分離され | 分離され |
| 離婚した | |
| 寡婦 | 寡婦 |
# 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")))))
各グループ内の人数を確認できます。
table(recast_data$marital.status)
出力：
## ## Married Not_married Separated Widow ## 21165 15359 7727 1286
ステップ 4) 要約統計量
ターゲット変数に関する統計を確認してみましょう。 以下のグラフでは、性別に応じて 50 を超える収入を得ている個人の割合をカウントしています。
# Plot gender income
ggplot(recast_data, aes(x = gender, fill = income)) +
geom_bar(position = "fill") +
theme_classic()
出力：
次に、個人の出身地が収入に影響を与えるかどうかを確認します。
# 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))
出力：
性別別の労働時間数。
# 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()
出力：
箱ひげ図は、労働時間の分布がさまざまなグループに当てはまることを表しています。箱ひげ図では、男女ともに均一な観察結果ではありません。
教育の種類ごとに週の労働時間の密度を確認できます。 ディストリビューションには多くの異なるピックがあります。 おそらく、米国における契約の種類によって説明できるでしょう。
# Plot distribution working time by education
ggplot(recast_data, aes(x = hours.per.week)) +
geom_density(aes(color = education), alpha = 0.5) +
theme_classic()
コードの説明
- ggplot(recast_data, aes( x= hours.per.week)): 密度プロットには1つの変数のみが必要です
- geom_density(aes(color = education), alpha =0.5): 密度を制御する幾何学オブジェクト
出力：
自分の考えを確認するには、一方通行を実行できます。 分散分析テスト:
anova <- aov(hours.per.week~education, recast_data) summary(anova)
出力：
## 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
ANOVA テストにより、グループ間の平均の差を確認します。
非線形性
モデルを実行する前に、労働時間数が年齢に関連しているかどうかを確認できます。
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()
コードの説明
- ggplot(recast_data, aes(x = age, y = hours.per.week)): グラフの美観を設定する
- geom_point(aes(color=収入), size =0.5): ドットプロットを構築します
- stat_smooth(): 次の引数を使用してトレンド ラインを追加します。
- Method='lm': 次の場合に近似値をプロットします。 線形回帰
- Formula = y~poly(x,2): 多項式回帰を当てはめます
- se = TRUE: 標準エラーを追加します。
- aes(color= 収入): 収入によってモデルを分割します。
出力：
一言で言えば、モデル内の交互作用項をテストして、週の労働時間と他の特徴の間の非線形効果を検出できます。 どのような条件で作業時間が異なるのかを把握することが重要です。
相関
次のチェックは、変数間の相関関係を視覚化することです。 スピアマン法で計算された相関係数を含むヒート マップをプロットできるように、因子水準タイプを数値に変換します。
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")
コードの説明
- data.frame(lapply(recast_data,as.integer)): データを数値に変換します
- ggcorr() は次の引数を使用してヒートマップをプロットします。
- Method: 相関を計算するメソッド
- nbreaks = 6: ブレーク数
- hjust = 0.8: プロット内の変数名の制御位置
- label = TRUE: ウィンドウの中央にラベルを追加します
- label_size = 3: ラベルのサイズ
- color = “grey50”): ラベルの色
出力：
ステップ 5) トレーニング/テスト セット
任意の監視付き 機械学習 タスクでは、データをトレイン セットとテスト セットの間で分割する必要があります。 他の教師あり学習チュートリアルで作成した「関数」を使用して、トレーニング/テスト セットを作成できます。
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)
出力：
## [1] 36429 9
dim(data_test)
出力：
## [1] 9108 9
ステップ 6) モデルを構築する
アルゴリズムがどのように実行されるかを確認するには、glm() パッケージを使用します。 の 一般化線形モデル モデルのコレクションです。 基本的な構文は次のとおりです。
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")
これで、一連の特徴間で収入レベルを分割するためのロジスティック モデルを推定する準備が整いました。
formula <- income~. logit <- glm(formula, data = data_train, family = 'binomial') summary(logit)
コードの説明
- 式 <- 収入 ~ .: 適合するモデルを作成する
- logit <- glm(formula, data = data_train, family = 'binomial'): ロジスティック モデル (family = 'binomial') を data_train データで近似します。
- summary(logit): モデルの概要を出力します。
出力：
## ## 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
モデルの概要から興味深い情報が明らかになります。 ロジスティック回帰のパフォーマンスは、特定の主要な指標を使用して評価されます。
- AIC（赤池情報量基準）：これは以下と同等です。 R2 ロジスティック回帰では。 パラメータの数にペナルティを適用した場合の適合度を測定します。 より小さい AIC 値は、モデルが真実に近いことを示します。
- Null逸脱:切片のみを使用してモデルを適合させます。 自由度は n-1 です。 これをカイ二乗値 (実際の値の仮説検定とは異なる近似値) として解釈できます。
- 残差逸脱: すべての変数を含むモデル。 これは、カイ二乗仮説検定としても解釈されます。
- フィッシャー スコアリングの反復数: 収束するまでの反復数。
glm() 関数の出力はリストに保存されます。 以下のコードは、ロジスティック回帰を評価するために構築したロジット変数で使用できるすべての項目を示しています。
# リストは非常に長いので、最初の XNUMX つの要素のみを出力します
lapply(logit, class)[1:3]
出力：
## $coefficients ## [1] "numeric" ## ## $residuals ## [1] "numeric" ## ## $fitted.values ## [1] "numeric"
各値は、$ 記号の後ろにメトリクスの名前を付けることで抽出できます。 たとえば、モデルを logit として保存したとします。 AIC 基準を抽出するには、以下を使用します。
logit$aic
出力：
## [1] 27086.65
ステップ 7) モデルのパフォーマンスを評価する
混乱マトリックス
この 混同行列 これまでに確認したさまざまなメトリクスと比較して、分類パフォーマンスを評価するのに適した選択肢です。 一般的な考え方は、True インスタンスが False に分類される回数をカウントすることです。
混同行列を計算するには、実際のターゲットと比較できるように、まず一連の予測を用意する必要があります。
predict <- predict(logit, data_test, type = 'response') # confusion matrix table_mat <- table(data_test$income, predict > 0.5) table_mat
コードの説明
- detect(logit,data_test, type = 'response'): テスト セットの予測を計算します。 type = 'response' を設定して、応答確率を計算します。
- table(data_test$income, detect > 0.5): 混同行列を計算します。 予測 > 0.5 は、予測された確率が 1 を超える場合は 0.5 を返し、それ以外の場合は 0 を返すことを意味します。
出力：
## ## FALSE TRUE ## <=50K 6310 495 ## >50K 1074 1229
混同行列の各行は実際のターゲットを表し、各列は予測されたターゲットを表します。 この行列の最初の行は、収入が 50 未満であるとみなします (False クラス)。6241 人が収入が 50 未満の個人として正しく分類されました (真陰性)、残りの 50 つは誤って XNUMX を超えるものとして分類されました (偽陽性）。 50 行目は 1229 を超える収入を考慮しており、肯定的なクラスは XNUMX でした (真陽性）、 真陰性 1074でした。
モデルを計算できます 精度 観測値全体にわたる真の陽性と真の陰性を合計することによって
accuracy_Test <- sum(diag(table_mat)) / sum(table_mat) accuracy_Test
コードの説明
- sum(diag(table_mat)): 対角線の合計
- sum(table_mat): 行列の合計。
出力：
## [1] 0.8277339
このモデルには、偽陰性の数を過大評価するという XNUMX つの問題があるようです。 これを 精度テストのパラドックス。 精度とは、ケースの総数に対する正しい予測の比率であると述べました。 精度は比較的高いものの、役に立たないモデルを作成できます。 それは支配的なクラスがあるときに起こります。 混同マトリックスを振り返ると、ほとんどのケースが真陰性として分類されていることがわかります。 ここで、モデルがすべてのクラスをネガティブ (つまり、50k 未満) として分類したと想像してください。 精度は 75% (6718/6718+2257) になります。 モデルのパフォーマンスは向上しましたが、真の陽性と真の陰性を区別するのに苦労しています。
このような状況では、より簡潔な指標を持つことが望ましいです。 次のことを確認できます。
- 精度=TP/(TP+FP)
- リコール=TP/(TP+FN)
精度と再現率
精度 ポジティブな予測の精度を調べます。 リコール 分類子によって正しく検出された陽性インスタンスの比率です。
これら XNUMX つのメトリクスを計算する XNUMX つの関数を構築できます。
- 精度を構築する
precision <- function(matrix) {
# True positive
tp <- matrix[2, 2]
# false positive
fp <- matrix[1, 2]
return (tp / (tp + fp))
}
コードの説明
- mat[1,1]: データ フレームの最初の列の最初のセル、つまり真陽性を返します。
- マット[1,2]; データ フレームの XNUMX 列目の最初のセル、つまり偽陽性を返します。
recall <- function(matrix) {
# true positive
tp <- matrix[2, 2]# false positive
fn <- matrix[2, 1]
return (tp / (tp + fn))
}
コードの説明
- mat[1,1]: データ フレームの最初の列の最初のセル、つまり真陽性を返します。
- マット[2,1]; データ フレームの最初の列の XNUMX 番目のセル、つまり偽陰性を返します。
機能をテストできます
prec <- precision(table_mat) prec rec <- recall(table_mat) rec
出力：
## [1] 0.712877 ## [2] 0.5336518
モデルが 50 を超える個人であると主張する場合、それが正しいのはケースの 54 パーセントのみであり、ケースの 50 パーセントで 72 を超える個人であると主張できます。
あなたは作成することができます 
精度と再現率に基づいてスコアを付けます。 の はこれら XNUMX つのメトリクスの調和平均であり、低い値により多くの重みが与えられることを意味します。
f1 <- 2 * ((prec * rec) / (prec + rec)) f1
出力：
## [1] 0.6103799
精度と再現率のトレードオフ
高精度と高再現率を両立させることは不可能です。
精度を高めると、正しい個人がより適切に予測されますが、多くの個人を見逃すことになります (再現率が低下します)。 状況によっては、再現率よりも高い精度が優先されます。 精度と再現率の間には凹型の関係があります。
- 想像してみてください。患者が病気に罹患しているかどうかを予測する必要があるとします。 できるだけ正確に行う必要があります。
- 顔認識によって街頭で詐欺師の可能性のある人々を検出する必要がある場合、精度は低くても、詐欺師として分類された人々を多数捕まえる方がよいでしょう。 警察は詐欺行為をしていない人物を釈放できるだろう。
ROC曲線
この 受信機 Opera特性 カーブは、バイナリ分類で使用されるもう XNUMX つの一般的なツールです。 これは適合率/再現率曲線に非常に似ていますが、ROC 曲線は適合率と再現率をプロットする代わりに、真陽性率 (つまり再現率) と偽陽性率を示します。 偽陽性率は、誤って陽性として分類された陰性インスタンスの割合です。 これは、XNUMX から真陰性率を引いたものに等しくなります。 真陰性率とも呼ばれます 特異性。 したがって、ROC 曲線は次のようにプロットされます。 感度 (想起) 対 1 特異性
ROC 曲線をプロットするには、RORC というライブラリをインストールする必要があります。 condaで見つけることができます ライブラリ。 コードを入力できます:
conda install -cr r-rocr –yes
予測() 関数とパフォーマンス() 関数を使用して ROC をプロットできます。
library(ROCR) ROCRpred <- prediction(predict, data_test$income) ROCRperf <- performance(ROCRpred, 'tpr', 'fpr') plot(ROCRperf, colorize = TRUE, text.adj = c(-0.2, 1.7))
コードの説明
- 予測(predict, data_test$income): ROCR ライブラリは、入力データを変換する予測オブジェクトを作成する必要があります。
- Performance(ROCRpred, 'tpr','fpr'): グラフに生成する XNUMX つの組み合わせを返します。 ここでは、tpr と fpr が構築されます。 精度とリコールを一緒にプロットするには、「prec」、「rec」を使用します。
出力：
ステップ8） モデルを改善する
間の相互作用を使用して、モデルに非線形性を追加してみることができます。
- 年齢と週あたりの勤務時間
- 性別と週あたりの勤務時間。
両方のモデルを比較するにはスコアテストを使用する必要があります
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
出力：
## [1] 0.6109181
前回よりも若干スコアが上がっています。 スコアを更新するためにデータの作業を続けることができます。
まとめ
ロジスティック回帰をトレーニングする関数を以下の表にまとめます。
| パッケージ | DevOps Tools Engineer試験のObjective | 演算 | 引数 |
|---|---|---|---|
| – | トレーニング/テスト データセットの作成 | create_train_set() | データ、サイズ、トレイン |
| グルム | 一般化線形モデルをトレーニングする | glm() | 式、データ、ファミリー* |
| グルム | モデルを要約する | まとめ（） | フィットモデル |
| ベース | 予測する | predict（） | 近似モデル、データセット、タイプ = '応答' |
| ベース | 混同行列を作成する | テーブル（） | y、予測() |
| ベース | 精度スコアの作成 | sum(diag(table())/sum(table() | |
| ROCR | ROC の作成 : ステップ 1 予測の作成 | 予測（） | 予測()、y |
| ROCR | ROCの作成：ステップ2 パフォーマンスの作成 | パフォーマンス（） | 予測()、'tpr'、'fpr' |
| ROCR | ROC の作成 : ステップ 3 グラフをプロットする | プロット（） | パフォーマンス（） |
他の GLM モデルの種類は次のとおりです。
– 二項: (リンク = "ロジット")
– ガウス: (リンク = 「アイデンティティ」)
– ガンマ: (リンク = 「逆」)
– inverse.gaussian: (リンク = “1/mu^2”)
– ポアソン: (リンク = "ログ")
– 疑似: (リンク = 「同一性」、分散 = 「定数」)
– 準二項: (リンク = "ロジット")
– 準ポアソン: (リンク = "ログ")
