R言語におけるGLM：一般化線形モデルとロジスティック回帰

Rにおける一般化線形モデル（GLM）とは何ですか？

A 一般化線形モデル（GLM） 通常の線形回帰を拡張し、応答変数が正規分布以外の分布に従う場合にも対応できるようにする。 R内蔵のものを装着します glm() statsパッケージの関数。

すべてのGLMは、次の3つの要素によって定義されます。

• ランダム成分: 応答変数の確率分布は、指数型分布族（二項分布、ポアソン分布、ガンマ分布、ガウス分布など）から取得されます。
• 体系的な構成要素： 線形予測子とは、説明変数の重み付けされた組み合わせのことです。
• リンク機能: 応答の平均値を線形予測子に接続する関数。例えば、二値データの場合はロジットリンク関数、カウントデータの場合はログリンク関数など。

その構造がモデルを「一般化」するのです。正規分布の結果を強制するのではなく、正しい分布を宣言します。 家族 引数を指定すると、Rは最尤法によって係数を推定します。ロジスティック回帰は、二項分布族とロジットリンク関数を用いた一般化線形モデル（GLM）に過ぎないため、自然な出発点となります。

Rにおけるロジスティック回帰とは何ですか？

ロジスティック回帰は、クラス、つまり確率を予測するために使用されます。 ロジスティック回帰では、バイナリの結果を正確に予測できます。

多くの属性に基づいてローンが拒否されるか受け入れられるかを予測したいと想像してください。 ロジスティック回帰の形式は 0/1 です。 ローンが拒否された場合は y = 0、承認された場合は y = 1。

ロジスティック回帰モデルは、XNUMX つの点で線形回帰モデルとは異なります。

• まず第一に、ロジスティック回帰は従属変数として二値 (バイナリ) 入力のみを受け入れます (つまり、0 と 1 のベクトル)。
• 第二に、結果は確率的リンク関数と呼ばれるものを通してマッピングされます。 シグモイド S字型の形状による（ロジスティック）関数：

関数の出力は常に 0 から 1 の間です。下の画像を確認してください。

シグモイド関数は 0 から 1 までの値を返します。分類タスクでは、0 または 1 の離散出力が必要です。

連続フローを離散値に変換するには、決定限界を 0.5 に設定します。 このしきい値を超えるすべての値は 1 として分類されます。

リンク関数が明確になったところで、一般化モデルを、既に知っている通常の線形モデルと比較してみましょう。

GLMと線形回帰：R値の主な違い

コードを書く前に、標準がいつ適用されるかを正確に知っておくと役立ちます。 線形回帰 関数 lm() はもはや適切ではなく、代わりに glm() を使用する必要があります。

要するに、結果が正規分布に従う測定値の場合は lm() を選択し、結果がイエス/ノーの決定の場合は glm() を選択します。そうでなければ、このタスクは別の方法に任せることができます。 決定木 分類器、イベント数、または平均値とともにばらつきが大きくなる厳密に正の値。

Rで一般化線形モデル（GLM）を作成する方法

理論的な説明が終わったところで、このチュートリアルの残りの部分では、二項分布に基づく一般化線形モデル（GLM）を実際のデータセットに最初から最後まで適用する方法を解説します。

を使ってみましょう 成人 ロジスティック回帰を説明するためのデータセット。「成人」データセットは、分類タスクに最適なデータセットです。目的は、個人の年間所得（米ドル）が50,000ドルを超えるかどうかを予測することです。このデータセットには、48,842件の観測値と10個の変数が含まれています。

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

Code 説明

• 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. 分布をプロットする

週あたりの労働時間の分布を詳しく見てみましょう

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

私たちは週平均労働時間の99パーセンタイル値を算出します。

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

Code 説明

• quantile(data_adult$hours.per.week, .99): 週あたりの労働時間の99パーセンタイルを計算します。

出力：

## 99% 
##  80

人口の99パーセントは週80時間未満で働いています。

このしきい値を超える観測値は削除できます。 のフィルターを使用します。 dplyr としょうかん。

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

出力：

## [1] 45537    10

2. 連続変数を標準化する

データのスケールが同じではないため、各列を標準化してパフォーマンスを向上させることができます。 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)

Code 説明

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

Code 説明

• data.frame(select_if(data_social, is.factor)): データ フレーム タイプの要因に要因列を格納します。 ライブラリ ggplot2 にはデータ フレーム オブジェクトが必要です。

出力：

## [1] 6

データセットには 6 つのカテゴリ変数が含まれています

2番目のステップはより複雑です。データフレームのfactorの各列に対して棒グラフを作成する必要があります。特に列数が多い場合は、このプロセスを自動化する方が便利です。

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 説明

• 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) 特徴量エンジニアリング

2つのカテゴリ変数には、モデルが必要とするレベル数よりも多くのレベルが含まれています。これらの変数を、より広範で、より多くのデータを含むカテゴリに再編成します。

リキャスト教育

上のグラフから、可変教育には 16 のレベルがあることがわかります。 これはかなりの量であり、一部のレベルでは観測値の数が比較的少なくなります。 この変数から取得できる情報量を増やしたい場合は、変数をより高いレベルに再キャストできます。 つまり、同じレベルの教育を受けたより大きなグループを作成します。 たとえば、教育レベルが低いと中退に変換されます。 高等教育レベルは修士に変更されます。

詳細は次のとおりです。

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 説明

• 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タルステータス

婚姻状況の下位レベルを作成することも可能です。次のコードでは、レベルを次のように変更します。

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

出力：

箱ひげ図は、労働時間の分布がさまざまなグループに当てはまることを表しています。箱ひげ図では、男女ともに均一な観察結果ではありません。

教育の種類別に週労働時間の密度を確認できます。分布には多くの異なるピークがあります。これはおそらく、教育の種類によって説明できます。trac米国ではt。

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

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

Code 説明

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

Code 説明

• 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) モデルを構築する

アルゴリズムの性能を確認するには、stats パッケージの 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)

Code 説明

• 式 <- 収入 ~ .: 適合するモデルを作成する
• 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"

各値はtrac$ 記号の後にメトリック名が続きます。たとえば、モデルを logit として保存した場合、tracAIC基準では、以下を使用します。

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

Code 説明

• 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	

混同行列の各行は実際のターゲットを表し、各列は予測されたターゲットを表します。この行列の最初の行は、50k未満の収入（ネガティブクラス）を考慮しています。6,310の観測値は、50k未満の収入を持つ個人として正しく分類されました（真陰性）、一方495件は50万以上と誤って分類された（偽陽性2行目は5万ドル以上の収入を考慮しています。1,229人が正しく識別されました（真陽性）一方、1,074人が見逃された（偽陰性).

モデルを計算できます 精度 観測値全体にわたる真の陽性と真の陰性を合計することによって

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

Code 説明

• sum(diag(table_mat)): 対角線の合計
• sum(table_mat): 行列の合計。

出力：

## [1] 0.8277339

このモデルには1つの問題があるようです。それは、偽陰性が多すぎるということです。これは 精度テストのパラドックス精度は、正しい予測の数と症例の総数の比率であると述べてきました。精度が比較的高くても、役に立たないモデルになることがあります。これは、支配的なクラスが存在する場合に起こります。混同行列を振り返ると、ほとんどの症例が真の陰性と分類されていることがわかります。ここで、モデルがすべての観測値を陰性（つまり、50k未満）と分類したと想像してみてください。それでも、精度は約75パーセント（6,805 / 9,108）になります。モデルのパフォーマンスは向上しますが、真の陽性と真の陰性を区別するのに苦労します。

このような状況では、より簡潔な指標を持つことが望ましいです。 次のことを確認できます。

• 精度=TP/(TP+FP)
• リコール=TP/(TP+FN)

精度と再現率

精度 ポジティブな予測の精度を調べます。 リコール 分類子によって正しく検出された陽性インスタンスの比率です。

これら XNUMX つのメトリクスを計算する XNUMX つの関数を構築できます。

1. 精度を構築する

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

Code 説明

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

Code 説明

• mat[1,1]: データ フレームの最初の列の最初のセル、つまり真陽性を返します。
• マット[2,1]; データ フレームの最初の列の XNUMX 番目のセル、つまり偽陰性を返します。

機能をテストできます

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

出力：

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

この2つの数値をよく読んでください。精度は0.71なので、モデルが個人の年収が5万ドル以上であると判断した場合、71%のケースで正しく判断しています。再現率は0.53なので、モデルは実際に5万ドル以上稼いでいる人のうち、53%しか検出できていません。

あなたは作成することができます 精度と再現率に基づいてスコアを付けます。 の はこれら XNUMX つのメトリクスの調和平均であり、低い値により多くの重みが与えられることを意味します。

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

出力：

## [1] 0.6103799

精度と再現率のトレードオフ

高精度と高再現率を両立させることは不可能です。

精度を高めると、正しい個人がより適切に予測されますが、多くの個人を見逃すことになります (再現率が低下します)。 状況によっては、再現率よりも高い精度が優先されます。 精度と再現率の間には凹型の関係があります。

• 想像してみてください。患者が病気に罹患しているかどうかを予測する必要があるとします。 できるだけ正確に行う必要があります。
• 顔認識によって街頭で詐欺師の可能性のある人々を検出する必要がある場合、精度は低くても、詐欺師として分類された人々を多数捕まえる方がよいでしょう。 警察は詐欺行為をしていない人物を釈放できるだろう。

ROC曲線

その 受信機 Opera特性 カーブは、バイナリ分類で使用されるもう XNUMX つの一般的なツールです。 これは適合率/再現率曲線に非常に似ていますが、ROC 曲線は適合率と再現率をプロットする代わりに、真陽性率 (つまり再現率) と偽陽性率を示します。 偽陽性率は、誤って陽性として分類された陰性インスタンスの割合です。 これは、XNUMX から真陰性率を引いたものに等しくなります。 真陰性率とも呼ばれます 特異性。 したがって、ROC 曲線は次のようにプロットされます。 感度 (想起) 対 1 特異性

ROC曲線をプロットするには、ROCRというパッケージをインストールする必要があります。これはcondaで見つけることができます。 ライブラリ。 コードを入力できます:

conda install -c r 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))

Code 説明

• 予測(predict, data_test$income): ROCR ライブラリは、入力データを変換する予測オブジェクトを作成する必要があります。
• Performance(ROCRpred, 'tpr','fpr'): グラフに生成する XNUMX つの組み合わせを返します。 ここでは、tpr と fpr が構築されます。 精度とリコールを一緒にプロットするには、「prec」、「rec」を使用します。

出力：

ステップ 8) モデルを改善する

間の相互作用を使用して、モデルに非線形性を追加してみることができます。

• 年齢と週あたりの勤務時間
• 性別と週あたりの勤務時間。

次に、両モデルのF1スコアを比較します。

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

F1スコアは前回よりもわずかに高くなっています。引き続きデータを分析して、スコア更新を目指してください。

RにおけるGLM係数とオッズ比の解釈方法

ステップ 6 で印刷される要約表には、係数が報告されます。 対数オッズ スケールは、専門知識のない人には説明しにくい。それをオッズ比に変換することで、モデルの説明がはるかに容易になる。

以下の4つの手順に従ってください。

1. 係数を指数関数化する。 すべての推定値にexp()関数を適用して、対数オッズを乗法的なオッズ比にします。
2. 信頼区間を追加する。 confint() を exp() で囲むと、同じオッズスケールで 95 パーセント区間を取得できます。
3. それぞれの値を1と比較してください。 オッズ比が1より大きい場合は陽性となる確率が高くなり、1より小さい場合は確率が低くなり、1に近い場合は予測因子がほとんど役に立たないことを意味します。
4. 統計的有意性を確認してください。 要約出力におけるp値が、選択した閾値（通常は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)))

出力結果を読み取っています。 私たちのモデルにおける週あたりの労働時間の係数は0.41942です。これを指数関数化するとexp(0.41942) = 1.52となり、これは他のすべての変数を一定に保った場合、週あたりの労働時間が標準偏差1つ分増加すると、年収5万ドル以上を稼ぐ確率が約1.5倍になることを意味します。

負の係数も同様に機能します。marital.statusNot_married は -2.50346 なので、exp(-2.50346) = 0.08 となります。つまり、未婚の人は既婚の人の約 8% の確率で結婚していることになります。連続予測変数はステップ 1 で標準化されているため、変化は生の時間ではなく、標準偏差単位で記述してください。

他のご家族への注意： 指数化された係数は、ロジットリンクを持つ二項分布族の場合にのみオッズ比として扱われます。分布族が「ポアソン」で、対数リンクを使用する場合は、同じ exp() の値は代わりに率比として読み取られます。

RにおけるGLM：クイック関数リファレンス

コーディング中は、この表を手元に置いておいてください。この表には、上記の8つの手順で使用されているすべての関数、それらを提供するパッケージ、および関数が要求する引数が記載されています。

他の GLM 家族間の議論を通じて利用できる家族構成は以下のとおりです。

• 二項分布: (リンク = “logit”)
• ガウス分布: (リンク = “identity”)
• ガンマ: (リンク = “逆数”)
• 逆ガウス分布: (リンク = “1/mu^2”)
• ポアソン: (リンク = “log”)
• 準: (リンク = “恒等写像”, 分散 = “定数”)
• 準二項分布: (リンク = “logit”)
• 準ポアソン: (リンク = “log”)

よくあるご質問