GLM dans R : Modèle linéaire généralisé et régression logistique

Qu'est-ce qu'un modèle linéaire généralisé (GLM) dans R ?

A Modèle linéaire généralisé (GLM) étend la régression linéaire ordinaire afin que la variable réponse puisse suivre une distribution autre que la distribution normale. R, vous en installez un avec le glm() fonction du package stats.

Chaque GLM est défini par trois composantes :

• Composante aléatoire : la distribution de probabilité de la variable de réponse, tirée de la famille exponentielle (binomiale, Poisson, Gamma, gaussienne et autres).
• Composante systématique : le prédicteur linéaire, c'est-à-dire la combinaison pondérée de vos variables explicatives.
• Fonction lien : la fonction qui relie la moyenne de la réponse au prédicteur linéaire, par exemple le lien logit pour les données binaires ou le lien logarithmique pour les comptages.

C’est cette structure qui confère au modèle son caractère « généralisé ». Au lieu d’imposer une distribution normale, on définit la distribution correcte par le biais de famille L'argument et R estiment les coefficients par maximum de vraisemblance. La régression logistique est simplement un modèle linéaire généralisé (GLM) avec une distribution binomiale et une fonction de lien logit ; c'est donc un point de départ naturel.

Qu'est-ce que la régression logistique dans R ?

La régression logistique est utilisée pour prédire une classe, c'est-à-dire une probabilité. La régression logistique peut prédire avec précision un résultat binaire.

Imaginez que vous souhaitiez prédire si un prêt est refusé/accepté en fonction de nombreux attributs. La régression logistique est de la forme 0/1. y = 0 si un prêt est rejeté, y = 1 s'il est accepté.

Un modèle de régression logistique diffère du modèle de régression linéaire de deux manières.

• Tout d'abord, la régression logistique n'accepte que les entrées dichotomiques (binaires) comme variable dépendante (c'est-à-dire un vecteur de 0 et 1).
• Deuxièmement, le résultat est modélisé par une fonction de lien probabiliste appelée sigmoïde (fonction logistique) due à sa forme en S :

La sortie de la fonction est toujours comprise entre 0 et 1. Vérifiez l'image ci-dessous

La fonction sigmoïde renvoie des valeurs de 0 à 1. Pour la tâche de classification, nous avons besoin d'une sortie discrète de 0 ou 1.

Pour convertir un flux continu en valeur discrète, nous pouvons fixer une limite de décision à 0.5. Toutes les valeurs supérieures à ce seuil sont classées 1

Maintenant que la fonction de lien est claire, comparez le modèle généralisé avec le modèle linéaire ordinaire que vous connaissez déjà.

GLM vs régression linéaire : principales différences dans R

Avant d'écrire du code, il est utile de savoir précisément quand la norme régression linéaire La fonction lm() n'est plus appropriée et la fonction glm() devrait la remplacer.

En résumé, choisissez lm() lorsque le résultat est une mesure distribuée normalement et glm() lorsque le résultat est une décision par oui ou par non, une tâche que vous pourriez autrement confier à un autre outil. arbre de décision classificateur, un décompte d'événements ou une quantité strictement positive dont la dispersion augmente avec sa moyenne.

Comment créer un modèle linéaire généralisé (GLM) dans R

La théorie étant établie, le reste de ce tutoriel applique un modèle linéaire généralisé binomial de bout en bout sur un jeu de données réel.

Utilisons le adulte Jeu de données illustrant la régression logistique. Le jeu de données « adulte » est idéal pour la tâche de classification. L’objectif est de prédire si le revenu annuel d’un individu, en dollars américains, dépassera 50 000 $. Le jeu de données contient 48 842 observations et dix variables.

• âge : âge de l’individu. Numérique
• éducation : Niveau d’éducation de l’individu. Facteur.
• état matrimonial : Maristatut réel de l’individu. Facteur c'est-à-dire jamais marié, marié-conjoint civil,…
• genre : Sexe de l’individu. Facteur, c'est-à-dire masculin ou féminin
• le revenu: Target variable. Revenu supérieur ou inférieur à 50K. Facteur c'est-à-dire >50K, <=50K

entre autres

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

Sortie :

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

Nous procéderons de la manière suivante :

• Étape 1 : Vérifiez les variables continues
• Étape 2 : Vérifier les variables factorielles
• Étape 3 : Ingénierie des fonctionnalités
• Étape 4 : Statistique récapitulative
• Étape 5 : Former/tester l’ensemble
• Étape 6 : Construire le modèle
• Étape 7 : Évaluer les performances du modèle
• Étape 8 : Améliorer le modèle

Votre tâche consiste à prédire quelle personne aura un revenu supérieur à 50 $.

Dans ce tutoriel, chaque étape sera détaillée pour effectuer une analyse sur un jeu de données réel.

Étape 1) Vérifiez les variables continues

Dans la première étape, vous pouvez voir la distribution des variables continues.

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

Code Explication

• continu <- select_if(data_adult, is.numeric) : utilisez la fonction select_if() de la bibliothèque dplyr pour sélectionner uniquement les colonnes numériques
• résumé (continu) : Imprimer la statistique récapitulative

Sortie :

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

Dans le tableau ci-dessus, vous pouvez voir que les données ont des échelles totalement différentes et que les heures par semaine présentent de grandes valeurs aberrantes (c'est-à-dire regardez le dernier quartile et la valeur maximale).

Vous pouvez y faire face en suivant deux étapes :

• Représentez graphiquement la distribution des heures par semaine
• Standardiser les variables continues

1. Tracer la distribution

Regardons de plus près la répartition des heures par semaine

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

Sortie :

La variable présente de nombreuses valeurs aberrantes et une distribution mal définie. Vous pouvez résoudre partiellement ce problème en supprimant les 0.01 pour cent des heures les plus élevées par semaine.

Syntaxe de base du quantile :

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.

Nous calculons le 99e centile des heures de travail hebdomadaires.

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

Code Explication

• quantile(data_adult$hours.per.week, .99): Calcule le 99e centile du temps de travail hebdomadaire

Sortie :

## 99% 
##  80

99 pour cent de la population travaille moins de 80 heures par semaine.

Vous pouvez déposer les observations au-dessus de ce seuil. Vous utilisez le filtre du déplyr bibliothèque.

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

Sortie :

## [1] 45537    10

2. Standardiser les variables continues

Vous pouvez standardiser chaque colonne pour améliorer les performances car vos données n'ont pas la même échelle. Vous pouvez utiliser la fonction mutate_if de la bibliothèque dplyr. La syntaxe de base est la suivante :

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

Vous pouvez standardiser les colonnes numériques comme suit :

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

Code Explication

• mutate_if(is.numeric, funs(scale)) : La condition est uniquement une colonne numérique et la fonction est une échelle

Sortie :

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

Étape 2) Vérifiez les variables factorielles

Cette étape a deux objectifs :

• Vérifiez le niveau dans chaque colonne catégorielle
• Définir de nouveaux niveaux

Nous diviserons cette étape en trois parties :

• Sélectionnez les colonnes catégorielles
• Stocker le graphique à barres de chaque colonne dans une liste
• Imprimer les graphiques

Nous pouvons sélectionner les colonnes de facteurs avec le code ci-dessous :

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

Code Explication

• data.frame(select_if(data_adult, is.factor)) : Nous stockons les colonnes de facteurs dans factor dans un type de trame de données. La bibliothèque ggplot2 nécessite un objet data frame.

Sortie :

## [1] 6

L'ensemble de données contient 6 variables catégorielles

La deuxième étape est plus complexe. Il s'agit de créer un graphique à barres pour chaque colonne du tableau de données « factor ». Il est plus pratique d'automatiser ce processus, surtout lorsqu'il y a de nombreuses colonnes.

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 Explication

• lapply() : Utilisez la fonction lapply() pour passer une fonction dans toutes les colonnes de l'ensemble de données. Vous stockez la sortie dans une liste
• function(x) : La fonction sera traitée pour chaque x. Ici x sont les colonnes
• ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)) : Créez un graphique à barres pour chaque élément x. Notez que pour renvoyer x sous forme de colonne, vous devez l'inclure dans get()

La dernière étape est relativement simple. Vous souhaitez imprimer les 6 graphiques.

# Print the graph
graph

Sortie :

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Remarque : utilisez le bouton suivant pour accéder au graphique suivant.

Étape 3) Ingénierie des fonctionnalités

Deux variables catégorielles comportent plus de modalités que nécessaire dans le modèle. Vous allez les regrouper en catégories plus larges et mieux représentées.

Refondre l'éducation

Sur le graphique ci-dessus, vous pouvez voir que la variable éducation comporte 16 niveaux. C'est important et certains niveaux ont un nombre d'observations relativement faible. Si vous souhaitez améliorer la quantité d'informations que vous pouvez obtenir de cette variable, vous pouvez la redéfinir à un niveau supérieur. Autrement dit, vous créez des groupes plus grands avec un niveau d’éducation similaire. Par exemple, un faible niveau d’éducation se transformera en abandon scolaire. Les niveaux d'enseignement supérieurs seront remplacés par le master.

Voici le détail :

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 Explication

• Nous utilisons le verbe muter de la bibliothèque dplyr. Nous changeons les valeurs de l'éducation avec la déclaration ifelse

Dans le tableau ci-dessous, vous créez une statistique récapitulative pour voir, en moyenne, combien d'années d'études (valeur z) il faut pour atteindre le baccalauréat, le master ou le doctorat.

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

Sortie :

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

Refonte Maristatut-tal

Il est également possible de créer des niveaux inférieurs pour l'état civil. Dans le code suivant, vous modifiez le niveau comme suit :

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

Vous pouvez vérifier le nombre d'individus au sein de chaque groupe.

table(recast_data$marital.status)

Sortie :

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

Étape 4) Statistique récapitulative

Il est temps de vérifier quelques statistiques sur nos variables cibles. Dans le graphique ci-dessous, vous comptez le pourcentage d'individus gagnant plus de 50 $ en fonction de leur sexe.

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

Sortie :

Ensuite, vérifiez si l’origine de l’individu affecte ses revenus.

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

Sortie :

Le nombre d’heures travaillées par sexe.

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

Sortie :

Le diagramme en boîte confirme que la répartition du temps de travail correspond à différents groupes. Dans le box plot, les deux sexes n’ont pas d’observations homogènes.

Vous pouvez consulter la densité du temps de travail hebdomadaire par type d'enseignement. Les distributions présentent de nombreuses variations. Cela peut probablement s'expliquer par le type de formation.tract aux États-Unis.

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

• ggplot(recast_data, aes( x= hours.per.week)) : un tracé de densité ne nécessite qu'une seule variable
• geom_density(aes(color = education), alpha =0.5) : L'objet géométrique pour contrôler la densité

Sortie :

Pour confirmer vos pensées, vous pouvez effectuer un aller simple Test ANOVA:

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

Sortie :

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

Le test ANOVA confirme la différence de moyenne entre les groupes.

Non-linéarité

Avant d'exécuter le modèle, vous pouvez voir si le nombre d'heures travaillées est lié à l'âge.

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 Explication

• ggplot(recast_data, aes(x = age, y = hours.per.week)) : Définir l'esthétique du graphique
• geom_point(aes(color= revenue), size =0.5) : Construire le diagramme de points
• stat_smooth() : ajoutez la ligne de tendance avec les arguments suivants :
  • method='lm' : tracez la valeur ajustée si le régression linéaire
  • formule = y~poly(x,2) : Ajuster une régression polynomiale
  • se = TRUE : Ajouter l'erreur standard
  • aes(color= revenue) : décomposer le modèle par revenu

Sortie :

En un mot, vous pouvez tester les termes d'interaction dans le modèle pour détecter l'effet de non-linéarité entre le temps de travail hebdomadaire et d'autres caractéristiques. Il est important de détecter dans quelles conditions le temps de travail diffère.

Corrélation

La vérification suivante consiste à visualiser la corrélation entre les variables. Vous convertissez le type de niveau de facteur en numérique afin de pouvoir tracer une carte thermique contenant le coefficient de corrélation calculé avec la méthode 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 Explication

• data.frame (lapply (recast_data, as.integer)): Convertir les données en numérique
• ggcorr() trace la carte thermique avec les arguments suivants :
  • method : Méthode pour calculer la corrélation
  • nbreaks = 6 : Nombre de pauses
  • hjust = 0.8 : Position de contrôle du nom de la variable dans le tracé
  • label = TRUE : Ajoute des étiquettes au centre des fenêtres
  • label_size = 3 : étiquettes de taille
  • color = « grey50 » : Couleur de l’étiquette

Sortie :

Étape 5) Entraîner/tester l’ensemble

Tout surveillé machine learning Cette tâche vous demande de répartir les données entre un ensemble d'entraînement et un ensemble de test. Vous pouvez utiliser la fonction que vous avez créée dans les autres tutoriels d'apprentissage supervisé pour créer ces ensembles.

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)

Sortie :

## [1] 36429     9

dim(data_test)

Sortie :

## [1] 9108    9

Étape 6) Construire le modèle

Pour observer les performances de l'algorithme, vous utilisez la fonction glm() du package stats. Modèle linéaire généralisé est une collection de modèles. La syntaxe de base est la suivante :

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

Vous êtes prêt à estimer le modèle logistique pour répartir le niveau de revenu entre un ensemble de fonctionnalités.

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

Code Explication

• formule <- revenu ~ .: Créer le modèle pour s'adapter
• logit <- glm(formula, data = data_train, family = 'binomial') : Ajustez un modèle logistique (family = 'binomial') avec les données data_train.
• summary(logit) : Imprimer le résumé du modèle

Sortie :

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

Le résumé de notre modèle révèle des informations intéressantes. Les performances d'une régression logistique sont évaluées avec des mesures clés spécifiques.

• AIC (Akaike Information Criteria) : C'est l'équivalent de R2 en régression logistique. Il mesure l’ajustement lorsqu’une pénalité est appliquée au nombre de paramètres. Plus petit AIC les valeurs indiquent que le modèle est plus proche de la vérité.
• Déviance nulle : s'adapte au modèle uniquement avec l'interception. Le degré de liberté est n-1. Nous pouvons l'interpréter comme une valeur du Chi carré (valeur ajustée différente du test d'hypothèse de valeur réelle).
• Déviance résiduelle : Modèle avec toutes les variables. Il est également interprété comme un test d’hypothèse du chi carré.
• Nombre d'itérations de Fisher Scoring : nombre d'itérations avant la convergence.

La sortie de la fonction glm() est stockée dans une liste. Le code ci-dessous montre tous les éléments disponibles dans la variable logit que nous avons construite pour évaluer la régression logistique.

# La liste est très longue, n'imprimez que les trois premiers éléments

lapply(logit, class)[1:3]

Sortie :

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

Chaque valeur peut être extracAjoutez le signe $ suivi du nom de la métrique. Par exemple, vous avez enregistré le modèle sous le nom logit. Pour en savoir plustracPour les critères AIC, vous utilisez :

logit$aic

Sortie :

## [1] 27086.65

Étape 7) Évaluer les performances du modèle

Matrice de confusion

Le matrice de confusion est un meilleur choix pour évaluer les performances de classification par rapport aux différentes mesures que vous avez vues auparavant. L'idée générale est de compter le nombre de fois où les instances vraies sont classées comme étant fausses.

Pour calculer la matrice de confusion, vous devez d’abord disposer d’un ensemble de prédictions afin qu’elles puissent être comparées aux cibles réelles.

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

Code Explication

• Predict(logit,data_test, type = 'response') : calcule la prédiction sur l'ensemble de test. Définissez type = 'response' pour calculer la probabilité de réponse.
• table(data_test$ Income, Predict > 0.5) : Calculez la matrice de confusion. prédire > 0.5 signifie qu'il renvoie 1 si les probabilités prédites sont supérieures à 0.5, sinon 0.

Sortie :

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

Chaque ligne d'une matrice de confusion représente une cible réelle, tandis que chaque colonne représente une cible prédite. La première ligne de cette matrice considère les revenus inférieurs à 50 000 (la classe négative) : 6 310 observations ont été correctement classées comme des individus ayant un revenu inférieur à 50 000 (Vrai négatif), tandis que 495 ont été classés à tort comme ayant un poids supérieur à 50k (Faux positif). La deuxième ligne prend en compte les revenus supérieurs à 50 000 : 1 229 ont été correctement identifiés (Vrai positif), tandis que 1 074 ont été manqués (Faux négatif).

Vous pouvez calculer le modèle précision en additionnant les vrais positifs + les vrais négatifs sur l'observation totale

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

Code Explication

• sum(diag(table_mat)) : Somme de la diagonale
• sum(table_mat) : Somme de la matrice.

Sortie :

## [1] 0.8277339

Le modèle semble souffrir d'un problème : il produit trop de faux négatifs. C'est ce qu'on appelle le paradoxe du test de précisionNous avons indiqué que la précision correspond au rapport entre les prédictions correctes et le nombre total de cas. On peut obtenir une précision relativement élevée sans pour autant avoir un modèle inutilisable. Cela se produit lorsqu'une classe est dominante. Si l'on examine la matrice de confusion, on constate que la plupart des cas sont classés comme vrais négatifs. Imaginons maintenant que le modèle classe chaque observation comme négative (c'est-à-dire inférieure à 50 000). On obtiendrait tout de même une précision d'environ 75 % (6 805 / 9 108). Le modèle est plus performant, mais peine à distinguer les vrais positifs des vrais négatifs.

Dans une telle situation, il est préférable d’avoir une métrique plus concise. Nous pouvons regarder :

• Précision=TP/(TP+FP)
• Rappel=TP/(TP+FN)

Précision vs rappel

La précision examine l’exactitude de la prédiction positive. Rappeler est le ratio d'instances positives correctement détectées par le classificateur ;

Vous pouvez construire deux fonctions pour calculer ces deux métriques

1. Construire avec précision

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

Code Explication

• mat[1,1] : Renvoie la première cellule de la première colonne de la trame de données, c'est à dire le vrai positif
• tapis[1,2]; Renvoie la première cellule de la deuxième colonne de la trame de données, c'est à dire le faux positif

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

Code Explication

• mat[1,1] : Renvoie la première cellule de la première colonne de la trame de données, c'est à dire le vrai positif
• tapis[2,1]; Renvoie la deuxième cellule de la première colonne de la trame de données, c'est à dire le faux négatif

Vous pouvez tester vos fonctions

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

Sortie :

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

Analysez attentivement ces deux chiffres. La précision est de 0.71, ce qui signifie que lorsque le modèle indique qu'une personne gagne plus de 50 000 €, il a raison dans 71 % des cas. Le rappel est de 0.53, ce qui signifie que le modèle ne détecte que 53 % des personnes qui gagnent réellement plus de 50 000 €.

Vous pouvez créer le score basé sur la précision et le rappel. Le est une moyenne harmonique de ces deux mesures, ce qui signifie qu'elle donne plus de poids aux valeurs les plus faibles.

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

Sortie :

## [1] 0.6103799

Compromis entre précision et rappel

Il est impossible d’avoir à la fois une haute précision et un haut rappel.

Si nous augmentons la précision, l’individu correct sera mieux prédit, mais nous en manquerons beaucoup (rappel plus faible). Dans certaines situations, nous préférons une précision supérieure au rappel. Il existe une relation concave entre précision et rappel.

• Imaginez, vous devez prédire si un patient souffre d'une maladie. Vous voulez être aussi précis que possible.
• Si vous devez détecter des personnes frauduleuses potentielles dans la rue grâce à la reconnaissance faciale, il serait préférable d'attraper de nombreuses personnes étiquetées comme frauduleuses, même si la précision est faible. La police pourra libérer l'individu non frauduleux.

La courbe ROC

Le Cible OperaCaractéristique de réglage La courbe est un autre outil couramment utilisé avec la classification binaire. Elle est très similaire à la courbe précision/rappel, mais au lieu de tracer la précision par rapport au rappel, la courbe ROC montre le taux de vrais positifs (c'est-à-dire le rappel) par rapport au taux de faux positifs. Le taux de faux positifs est le rapport des instances négatives qui sont incorrectement classées comme positives. Il est égal à un moins le taux véritablement négatif. Le vrai taux négatif est également appelé spécificité. D'où la courbe ROC Un niveau de sensibilité élevée (rappel) versus 1-spécificité

Pour tracer la courbe ROC, il faut installer le paquet ROCR. On le trouve dans le répertoire conda. bibliothèque math. Vous pouvez taper le code :

conda install -c r r-rocr --yes

Nous pouvons tracer le ROC avec les fonctions prédiction() et performance().

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 Explication

• prédiction(predict, data_test$ Income) : la bibliothèque ROCR doit créer un objet de prédiction pour transformer les données d'entrée
• performance(ROCRpred, 'tpr','fpr') : Renvoie les deux combinaisons à produire dans le graphique. Ici, tpr et fpr sont construits. Pour tracer la précision et rappeler ensemble, utilisez « prec », « rec ».

Sortie :

Étape 8) Améliorer le modèle

Vous pouvez essayer d'ajouter de la non-linéarité au modèle avec l'interaction entre

• âge et heures par semaine
• sexe et heures par semaine.

Vous comparez ensuite le score F1 des deux modèles.

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

Sortie :

## [1] 0.6109181

Le score F1 est légèrement supérieur au précédent. Vous pouvez continuer à travailler sur les données et essayer de l'améliorer.

Comment interpréter les coefficients et les rapports de cotes des modèles linéaires généralisés (GLM) dans R ?

Le tableau récapitulatif imprimé à l'étape 6 présente les coefficients sur le log-cotes L'échelle est difficile à expliquer à un public non spécialisé. La conversion en rapports de cotes rend le modèle beaucoup plus facile à communiquer.

Suivez ces quatre étapes.

1. Exponentier les coefficients. Appliquez la fonction exp() à chaque estimation afin que les logarithmes des cotes deviennent des rapports de cotes multiplicatifs.
2. Ajouter un intervalle de confiance. Enveloppez confint() dans exp() pour obtenir l'intervalle de 95 % sur la même échelle de cotes.
3. Comparez chaque valeur à 1. Un rapport de cotes supérieur à 1 augmente la probabilité de la classe positive, une valeur inférieure à 1 la diminue et une valeur proche de 1 signifie que le prédicteur n'apporte que peu de choses.
4. Vérifier la signification statistique. N’interprétez que les prédicteurs dont la valeur p dans le résumé est inférieure au seuil que vous avez choisi, généralement 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)))

Lecture du résultat. Le coefficient pour les heures par semaine dans notre modèle est de 0.41942. En l'exponentiant, on obtient exp(0.41942) = 1.52, ce qui signifie qu'une augmentation d'un écart-type des heures de travail hebdomadaires multiplie par environ 1.5 les chances de gagner plus de 50 000 $, toutes choses égales par ailleurs.

Les coefficients négatifs fonctionnent de la même manière. marital.statusNot_married vaut -2.50346, donc exp(-2.50346) = 0.08 : les personnes célibataires ont environ 8 % des chances d’être mariées. Étant donné que les prédicteurs continus ont été standardisés à l’étape 1, décrivez les variations en unités d’écart type, et non en heures brutes.

Note aux autres familles : Les coefficients exponentiés représentent des rapports de cotes uniquement dans le cadre d'une loi binomiale avec une fonction de lien logit. Avec une loi de Poisson et une fonction de lien logarithmique, les mêmes valeurs exp() sont interprétées comme des rapports de taux.

GLM dans R : Guide rapide des fonctions

Gardez ce tableau à portée de main pendant que vous codez. Il répertorie chaque fonction utilisée dans les huit étapes ci-dessus, ainsi que le package qui la fournit et les arguments attendus.

L'autre GLM Les familles disponibles via l'argument familial sont :

• binomiale : (lien = « logit »)
• gaussienne : (lien = « identité »)
• Gamma : (lien = « inverse »)
• inverse.gaussienne : (lien = « 1/mu^2 »)
• poisson : (lien = « log »)
• quasi : (lien = « identité », variance = « constante »)
• quasi-binomial : (lien = « logit »)
• quasipoisson : (lien = « logarithme »)

FAQ