Statistical Analysis

ANALISI DELLA REGRESSIONE LOGISTICA IN R

Useremo il dataset DEFAULT presente nel pacchetto ISLR

PACCHETTI PER CONDURRE L’ANALISI

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.0.3

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.0.4 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 4.0.4
## Warning: package 'tibble' was built under R version 4.0.3
## Warning: package 'dplyr' was built under R version 4.0.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ISLR)

## Warning: package 'ISLR' was built under R version 4.0.3

library(dplyr)
library(ggplot2)

Sintesi dei dati DEFAULT

summary(Default)

## default student balance income

## No :9667 No :7056 Min. : 0.0 Min. : 772
## Yes: 333 Yes:2944 1st Qu.: 481.7 1st Qu.:21340
## Median : 823.6 Median :34553
## Mean : 835.4 Mean :33517
## 3rd Qu.:1166.3 3rd Qu.:43808
## Max. :2654.3 Max. :73554
dati <- Default # i livelli No e Si sono ricodificati 1 e 0
dati <- dati %>%
select(default, balance,income,student) %>%
mutate(default = recode(default,
"No" = 0,

1
"Yes" = 1))
head(dati)

## default balance income student

## 1 0 729.5265 44361.625 No
## 2 0 817.1804 12106.135 Yes
## 3 0 1073.5492 31767.139 No
## 4 0 529.2506 35704.494 No
## 5 0 785.6559 38463.496 No
## 6 0 919.5885 7491.559 Yes

Adattiamo un modello di regressione lineare

modello_lineare <- lm(default ~ balance, data = dati)

RAPPRESENTAZIONE GRAFICA DEL MODELLO

ggplot(data = dati, aes(x = balance, y = default)) +

geom_point(aes(color = as.factor(default)), shape = 2) +
geom_smooth(method = "lm", color = "gray20", se = FALSE) +
theme_bw() +
labs(title = "Regressione lineare",
y = "Probabilità di deafult") +
theme(legend.position = "none")

## `geom_smooth()` using formula 'y ~ x'

Regressione lineare
1.00

0.75
Probabilità di deafult

0.50

0.25

0.00

0 1000 2000
balance

2
Trattandosi di una linea, se, ad esempio, la probabilità di default è prevista per qualcuno che ha un saldo di
10000, il valore ottenuto è maggiore di 1.
predict(object = modello_lineare, newdata = data.frame(balance = 10000))

## 1
## 1.22353
Per evitare questi problemi, la regressione logistica trasforma il valore restituito dalla regressione lineare
(bo + b1X) utilizzando una funzione il cui risultato è sempre compreso tra 0 e 1. Esistono diverse funzioni
che soddisfano questa descrizione, una delle più utilizzate è la funzione logistica noto anche come funzione
sigmoide)

ADATTIAMO UN MODELLO LOGISTICO

La regressione logistica in R è implementata con la funzione glm() Creiamo l’oggetto modello_logistico che
contiene i risultati della procedura glm() sui dati Default. L’opzione family=binomial produce un modello di
regressione logistica.
modello_logistico <- glm(default ~ balance, data = dati, family = "binomial")

RAPPRESENTAZIONE GRAFICA DEL MODELLO

I MODO
ggplot(data = dati, aes(x = balance, y = default)) +
geom_point(aes(color = as.factor(default)), shape = 1) +
stat_function(fun = function(x){predict(modello_logistico,
newdata = data.frame(balance = x),
type = "response")}) +
theme_bw() +
labs(title = "Regressione logistica",
y = "Probabilità di default") +
theme(legend.position = "none")

3
 Regressione logistica
1.00

0.75
Probabilità di default

0.50

0.25

0.00

0 1000 2000
balance

II MODO CON LA FUNZIONE geom_smooth CHE CONSENTE DI OTTENERE IL GRAFICO DIRET-

TAMENTE
ggplot(data = dati, aes(x = balance, y = default)) +
geom_point(aes(color = as.factor(default)), shape = 1) +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
color = "gray20",
se = FALSE) +
theme_bw() +
theme(legend.position = "none")

## `geom_smooth()` using formula 'y ~ x'

4
 1.00

0.75
default

0.50

0.25

0.00

0 1000 2000
balance

OUTPUT MODELLO DI REGRESSIONE LOGISTICA

summary(modello_logistico)

##
## Call:
## glm(formula = default ~ balance, family = "binomial", data = dati)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2697 -0.1465 -0.0589 -0.0221 3.7589
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.065e+01 3.612e-01 -29.49 <2e-16 ***
## balance 5.499e-03 2.204e-04 24.95 <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: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1596.5 on 9998 degrees of freedom
## AIC: 1600.5
##
## Number of Fisher Scoring iterations: 8

5
La funzione summary () restituisce la stima, gli errori standard, il punteggio z e ed i p-value per ciascuno
dei coefficienti. Fornisce anche la devianza nulla (la devianza solo per la media) e la devianza residua (la
devianza per il modello con tutti i predittori).
E’ possibile ottenere singole parti dell’oggetto modello_logistico.
coef(modello_logistico)

## (Intercept) balance
## -10.651330614 0.005498917
summary(modello_logistico)$coef

## Estimate Std. Error z value Pr(>|z|)

## (Intercept) -10.651330614 0.3611573721 -29.49221 3.623124e-191
## balance 0.005498917 0.0002203702 24.95309 1.976602e-137
#pvalue
summary(modello_logistico)$coef[,4]

## (Intercept) balance
## 3.623124e-191 1.976602e-137

AGGIUNGIAMO LE VARIABILI INCOME E STUDENT

modello_logistico_2 = glm(default ~ balance + student, data= dati, family=binomial)

summary(modello_logistico_2)

##
## Call:
## glm(formula = default ~ balance + student, family = binomial,
## data = dati)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4578 -0.1422 -0.0559 -0.0203 3.7435
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.075e+01 3.692e-01 -29.116 < 2e-16 ***
## balance 5.738e-03 2.318e-04 24.750 < 2e-16 ***
## studentYes -7.149e-01 1.475e-01 -4.846 1.26e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1571.7 on 9997 degrees of freedom
## AIC: 1577.7
##
## Number of Fisher Scoring iterations: 8

INTERPRETAZIONE DEL MODELLO

A differenza della regressione lineare, in cui b1 corrisponde alla variazione media della variabile dipendente
Y dovuta all’aumento di un’unità del predittore X, nella regressione logistica, b1 indica la variazione del

6
logaritmo dell’ ODDS a causa dell’aumento di un’unità X.
###STIMA INTERVALLI DI CONFIDENZA
confint(modello_logistico_2)

## Waiting for profiling to be done...

## 2.5 % 97.5 %
## (Intercept) -11.498144891 -10.049841998
## balance 0.005296591 0.006206095
## studentYes -1.007777857 -0.429075610

ODDS RATIO (OR) and 95% CI

Un concetto importante da comprendere, per interpretare i coefficienti beta logistici, è l’odds ratio. Un odds
ratio misura l’associazione tra una variabile predittiva (x) e la variabile risultato (y). Rappresenta il rapporto
tra le probabilità che si verifichi un evento (evento = 1) data la presenza del predittore x (x = 1), rispetto
alle probabilità che l’evento si verifichi in assenza di tale predittore (x = 0).
exp(cbind(OR = coef(modello_logistico_2), confint(modello_logistico_2)))

## Waiting for profiling to be done...

## OR 2.5 % 97.5 %
## (Intercept) 2.145622e-05 0.0000101489 4.319257e-05
## balance 1.005755e+00 1.0053106429 1.006225e+00
## studentYes 4.892520e-01 0.3650292262 6.511107e-01

CONVERTIRE IL LOGIT IN PROBABILITA’

Per convertire un logit (output glm) in probabilità, si possono seguire i 3 seguenti steps:
1) si prende il coefficiente di output glm (logit)
2) si prende l’exp del logit e si ottiene l’odds
3) si convertono gli odds in probabilità usando questa formula prob = odds / (1 + odds).
Ad esempio, se odds = 2/1, quindi la probabilità è 2 / (1 + 2) = 2/3 (~ .67)
odds <- exp(coef(modello_logistico_2)[3])
prob <- odds / (1 + odds)
prob

## studentYes
## 0.328522

REGOLE INTERPRETATIVE
Ricordiamo che exp(1) è circa 2.71 Cioè, se il logit è 1, l’odds sarà 2,7 quindi la probabilità (odds/(1+odds))
è 2,7 / 3,7, o circa 3/4, 75%.
Allo stesso modo importante, exp(0) = 1 Quindi, l’odds=1 saranno 1: 1, ovvero il 50%
Quindi, ogni volta che il logit è negativo, la probabilità associata è inferiore al 50% e v.v. (logit positivo,
probabilità superiore al 50%).

CONDIZIONI
1) Indipendenza: le osservazioni devono essere indipendenti l’una dall’altra.

7
 2) Relazione lineare tra il logaritmo naturale dell’odds ed una variabile continua: i modelli a forma di U
sono una chiara violazione di questa condizione.
3)La regressione logistica non richiede una distribuzione normale della variabile continua indipendente.
4)Numero di osservazioni: non esiste uno standard stabilito a questo riguardo, ma si raccomandano tra 50 e
100 osservazioni.

PREVISIONI
Si possono utilizzare le probabilità previste per comprendere il modello.
Le probabilità previste possono essere calcolate per variabili predittive sia categoriali che continue.
Per creare probabilità previste, dobbiamo prima creare un nuovo dataframe con i valori che vogliamo che le
variabili indipendenti assumano per creare le nostre previsioni.
Calcoliamo la probabilità di default per ogni valore di student, tenendo balance ed income ai valori medi.
Creiamo e visualizzazione il dataframe.
newdata1 <- with(dati, data.frame(balance = mean(balance), income = mean(income), student))
head(newdata1)

## balance income student

## 1 835.3749 33516.98 No
## 2 835.3749 33516.98 Yes
## 3 835.3749 33516.98 No
## 4 835.3749 33516.98 No
## 5 835.3749 33516.98 No
## 6 835.3749 33516.98 Yes
newdata1$studentp <- predict(modello_logistico_2, newdata = newdata1, type = "response")
head(newdata1)

## balance income student studentp

## 1 835.3749 33516.98 No 0.002583489
## 2 835.3749 33516.98 Yes 0.001265647
## 3 835.3749 33516.98 No 0.002583489
## 4 835.3749 33516.98 No 0.002583489
## 5 835.3749 33516.98 No 0.002583489
## 6 835.3749 33516.98 Yes 0.001265647