Modello statistico per prevedere il peso dei neonati
Aristide De Sanctis
2024-11-30
setwd("C:/Users/arist/Desktop/CORSI/PROFESSION_AI/Progetti/Progetto_3_STAT_INFERENZIALE")L’obiettivo di questo progetto è identificare quali, delle variabili presenti nel dataset, sono più predittive del peso alla nascita dei neonati, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione che potrebbero indicare nascite premature.
Carichiamo il dataset:
data = read.csv(file = "neonati.csv",
header = TRUE,
stringsAsFactors = TRUE)
n = nrow(data)
print(paste("Nel dataset sono presenti", n, "unità statistiche"))## [1] "Nel dataset sono presenti 2500 unità statistiche"Indichiamo che tipo di variabili sono contenute nel dataset:
Anni.madre -> età della madre espressa in anni (quantitativa discreta)
N.gravidanze -> n. gravidanze della madre (quantitativa discreta)
Fumatrici -> Variabile binaria, 0 indica non fumatrici, 1 fumatrici (qualitativa discreta)
Gestazione -> numero di settimane di gestazione (quantitativa discreta)
Peso -> indica il peso del neonato alla nascita espresso in grammi (quantitativa continua)
Lunghezza -> lunghezza del neonato (quantitativa continua)
Cranio -> diametro craniale, misurabili anche durante la gravidanza tramite ecografie (quantitativa continua)
Tipo.parto -> variabile binaria, può essere naturale Nat o cesareo Ces (qualitativa discreta)
Ospedale -> in quale ospedale è nato il bambino (qualitativa discreta) può essere osp1, osp2 o osp3
Sesso -> variabile binaria che indica il sesso del neonato (qualitativa discreta), può assumere M o F.
Analisi Preliminare
Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.
Un focus particolare sarà dato alla relazione tra le variabili materne (età, fumo, gravidanze precedenti, ecc…) e il peso del neonato.
Per svolgere questo task visualizzerò gli istogrammi e i boxplot con alcuni indici di posizione e variabilità.
Valuteremo cosa fare con eventuali valori anomali e in particolare per la variabile risposta vedremo se si distribuisce normalmente.
Effettuare queste verifiche ci serve per:
Test parametrici (come t-test, ANOVA, regressione lineare) spesso assumono che i dati siano normalmente distribuiti o che gli errori dei residui lo siano. Dati normalmente distribuiti garantiscono affidabilità nei risultati.
Inoltre la normalità facilita il calcolo di intervalli di confidenza, probabilità e test di ipotesi, utilizzando modelli teorici noti.
Se i dati non sono normali, potremmo dover usare test non parametrici o trasformazioni.
Influenza sulle statistiche: Gli outlier possono distorcere la media, la varianza e i risultati di test statistici, portando a conclusioni errate.
Comprensione del fenomeno: Un valore anomalo potrebbe indicare errori di misurazione, input errati o un fenomeno interessante da investigare (ad esempio, un evento raro).
Per la vasualizzazione grafica userò queste due funzioni scritte da me che danno in output rispettivamente un istogramma e un boxplot della variabile di interesse.
Queste verranno usate per variabili quantitative continue.
Userò barplot per variabili qualitative o quantitative discrete con poche modalità.
In seguito userò visualizzazioni e test di ipotesi per interpretare la relazione tra la variabile target (peso) e le altre variabili esplicative.
Funzione per Istogramma:
Un istogramma è un tipo di grafico utilizzato per rappresentare la distribuzione di una variabile quantitativa. Mostra come i dati sono distribuiti suddividendo i valori in intervalli (o bin) e contando il numero di osservazioni che rientrano in ciascun intervallo. È una forma di rappresentazione grafica che aiuta a visualizzare la densità e la forma della distribuzione dei dati.
histogram = function(variable,
x_legend = 1,
y_legend = 1,
cex = 1,
ylim_max = 1) {
# Crea un istogramma con annotazioni visive per la media e la mediana.
#
# Arguments:
# - variable: vettore numerico per cui calcolare e disegnare l'istogramma.
# - x_legend: coordinata X della posizione della legenda (default = 1).
# - y_legend: coordinata Y della posizione della legenda (default = 1).
# - cex: fattore di scala per la dimensione del testo nella legenda (default = 1).
# - ylim_max: lunghezza visualizzata dell'asse y
# Dettagli:
# - Disegna un istogramma con densità relativa (freq = FALSE) e sovrappone una curva di densità stimata.
# - Aggiunge linee verticali per indicare:
# - La mediana (linea rossa tratteggiata).
# - La media (linea blu tratteggiata).
# - La legenda mostra i valori calcolati per la media e la mediana.
# Esempio:
# set.seed(123)
# dati <- rnorm(100, mean = 50, sd = 10)
# histogram(variable = dati, x_legend = 60, y_legend = 0.04, cex = 0.8)
# Valori calcolati:
# - mean_var: Media della variabile.
# - median_var: Mediana della variabile.
# - min_var: Minimo della variabile.
# - max_var: Massimo della variabile.
# - q1: Primo quartile.
# - q3: Terzo quartile.
# Note:
# - La funzione assume che la variabile fornita sia numerica.
# - Per variabili con valori mancanti (`NA`), la funzione ignora automaticamente questi valori nei calcoli.
# calculate position indices
mean_var = mean(variable)
median_var = median(variable)
min_var = min(variable)
max_var = max(variable)
q1 <- quantile(variable, 0.25, na.rm = TRUE)
q3 <- quantile(variable, 0.75, na.rm = TRUE)
# Calcola densità stimata
dens <- density(variable)
# Trova i valori della densità per la media e la mediana
density_mean <- approx(dens$x,
dens$y,
xout = mean_var)$y
density_median <- approx(dens$x,
dens$y,
xout = median_var)$y
# Ottieni il nome della variabile
var_name <- deparse(substitute(variable))
# Se la variabile è una colonna di un data.frame, estrai il nome della colonna senza il prefisso 'data$'
if (grepl("\\$", var_name)) {
var_name = sub(".*\\$", "", var_name)
}
# crea istogramma
hist(variable,
main = paste("histogram of", var_name),
xlab = var_name,
ylab = "Relative frequencies",
col = "lightcyan",
border = "black",
breaks = 10,
freq = FALSE,
ylim = c(0, ylim_max)
)
lines(density(variable), col="red", lwd=1.5)
# Disegna segmenti verticali per media e mediana che si fermano alla curva di densità
segments(x0 = median_var, y0 = 0,
x1 = median_var, y1 = density_median,
col = "red",
lwd = 3,
lty = 2)
segments(x0 = mean_var, y0 = 0,
x1 = mean_var, y1 = density_mean,
col = "blue", lwd = 3, lty = 2)
# legend
legend(x = x_legend, y = y_legend,
legend = c(paste("Median =", round(median_var, 2)),
paste("Mean =", round(mean_var, 2))),
col = c("red", "blue"),
lwd = 2, lty = 2, title="Legend",
cex = cex, xpd = TRUE)
}Funzione per il BoxPlot:
Un boxplot, o diagramma a scatola e baffi, è un grafico usato per riassumere una variabile quantitativa, mostrando la distribuzione, la mediana, e i possibili outlier. È particolarmente utile per confrontare più distribuzioni.
Elementi principali di un boxplot:
Scatola (box): rappresenta il range interquartile (IQR), cioè il 50% centrale dei dati (dal primo quartile q1 al terzo quartile q3).
Una scatola larga indica una distribuzione più disperasa con maggiore variabilità nei dati
Una scatola stretta suggerisce meno variabilità
Linea centrale nella scatola: indica la mediana (il valore centrale dei dati).
Una mediana centrata indica simmetria della distribuzione
Una mediana spostata suggerisce una distribuzione asimmetrica
Baffi (whiskers): si estendono fino ai non considerati outlier. Generalmente:
- Si fermano a 1.5 volte l’IQR dal q1 o dal q3.
Outlier: sono rappresentati come punti individuali fuori dai baffi.
Più boxplot affiancati permettono di confrontare la mediana, la dispersione, e la presenza di outlier tra i gruppi.
# boxplot
boxplot_1 = function(variable,
x_legend = 0,
y_legend = 1.5,
cex = 1) {
# Crea un boxplot orizzontale con annotazioni visive dei principali valori statistici.
# Arguments:
# - variable: vettore numerico per cui calcolare e disegnare il boxplot.
# - x_legend: coordinata X della posizione della legenda (default = 0).
# - y_legend: coordinata Y della posizione della legenda (default = 1.5).
# - cex: fattore di scala per la dimensione del testo nella legenda (default = 1).
# Dettagli:
# - Disegna un boxplot orizzontale e aggiunge linee per il primo quartile (verde),
# la mediana (rosso) e il terzo quartile (blu).
# - Punti distintivi sono utilizzati per il minimo e il massimo (nero).
# Esempio:
# set.seed(123)
# dati <- rnorm(100, mean = 50, sd = 10)
# boxplot_1(variable = dati, x_legend = 55, y_legend = 1.2, cex = 0.8)
# Calcolare gli indici di posizione
mean_var = mean(variable)
median_var = median(variable)
min_var = min(variable)
max_var = max(variable)
q1 <- quantile(variable, 0.25, na.rm = TRUE)
q3 <- quantile(variable, 0.75, na.rm = TRUE)
# Ottieni il nome della variabile
var_name <- deparse(substitute(variable))
# Se il nome della variabile è del tipo "data$Anni.madre", estrai solo "Anni.madre"
if (grepl("\\$", var_name)) {
var_name <- sub(".*\\$", "", var_name)
}
# Crea il boxplot
boxplot(variable,
main = paste("Boxplot of", var_name),
xlab = var_name,
col = "lightcyan",
border = "black",
horizontal = TRUE # Boxplot orizzontale
)
# add lines for quartiles
segments(q1, 0.8, q1, 1.2, col = "green", lwd = 4)
segments(median_var, 0.8, median_var, 1.2, col = "red", lwd = 4)
segments(q3, 0.8, q3, 1.2, col = "blue", lwd = 4)
# Aggiungi i punti per il minimo e il massimo
points(min_var, 1, pch = 19, col = "black",
cex = 1.5) # Punto per il minimo
points(max_var, 1, pch = 17, col = "black",
cex = 1.5) # Punto per il massimo
legend(x = x_legend, y = y_legend,
legend = c(paste("Q1 =", q1),
paste("Median =", round(median_var, 2)),
paste("Q3 =", q3),
paste("Min =", round(min_var, 2)),
paste("Max =", round(max_var, 2))),
col = c("green", "red", "blue", "black", "black"),
pch = c(NA, NA, NA, 19, 17),
lwd = c(2, 2, 2, NA, NA),
title = "Legend",
cex = cex,
xpd = TRUE, # xpd ed inset mi permettono di estendere la legenda al di fuori del grafico
)
}Funzione boxplot condizionato:
boxplot_condizionato = function(var1, var2) {
# Descrizione
# La funzione boxplot_condizionato consente di creare due grafici boxplot affiancati
# per visualizzare la distribuzione di una variabile continua (var1).
# Il primo grafico rappresenta il boxplot semplice di var1, mentre il secondo rappresenta
# il boxplot di var1 condizionato a una variabile categoriale o fattoriale (var2).
# Sintassi:
# boxplot_condizionato(var1, var2)
# Argomenti:
# - var1: Una variabile numerica continua da rappresentare come boxplot.
# Può essere specificata direttamente o con un riferimento a un dataframe (es. data$var1).
# - var2: Una variabile categoriale o fattoriale utilizzata per condizionare il boxplot di var1.
# Anch'essa può essere specificata direttamente o con un riferimento a un dataframe (es. data$var2).
# Dettagli:
# 1. La funzione estrae il nome leggibile delle variabili utilizzando deparse(substitute()).
# Se la variabile è specificata con una sintassi come data$variabile, il nome viene ridotto a variabile.
# 2. Vengono creati due boxplot:
# - Il primo rappresenta solo la distribuzione di var1.
# - Il secondo condiziona var1 ai livelli di var2.
# 3. I grafici vengono disposti affiancati (par(mfrow = c(1, 2))).
# 4. I colori predefiniti sono:
# - Rosa per il primo boxplot.
# - Azzurro chiaro per il secondo.
# Esempio di utilizzo:
# data <- data.frame(
# Età = c(25, 30, 22, 28, 35, 40, 45, 50),
# Genere = factor(c("M", "F", "M", "F", "M", "F", "M", "F"))
# )
# boxplot_condizionato(data$Età, data$Genere)
# Output atteso:
# 1. Primo grafico: Boxplot della variabile var1.
# 2. Secondo grafico: Boxplot della variabile var1 condizionato a var2.
# Note:
# - var2 dovrebbe essere categoriale o fattoriale; altrimenti, potrebbero verificarsi errori.
# - La funzione è progettata per variabili in dataframe, ma accetta anche vettori semplici.
# Requisiti: R base; nessuna dipendenza da pacchetti esterni.
par(mfrow = c(1,2))
# Ottieni il nome della variabile target
var_name_1 <- deparse(substitute(var1))
# Se il nome della variabile è del tipo "data$Anni.madre", estrai solo "Anni.madre"
if (grepl("\\$", var_name_1)) {
var_name_1 <- sub(".*\\$", "", var_name_1)
}
# Crea il boxplot per var1
boxplot(var1,
main = paste("Boxplot of", var_name_1),
xlab = var_name_1,
col = "pink",
border = "black",
horizontal = FALSE # Boxplot orizzontale
)
# Ottieni il nome della variabile 2
var_name_2 <- deparse(substitute(var2))
# Se il nome della variabile è del tipo "data$Anni.madre", estrai solo "Anni.madre"
if (grepl("\\$", var_name_2)) {
var_name_2 <- sub(".*\\$", "", var_name_2)
}
# Crea il boxplot per var2
boxplot(var1 ~ var2,
main = paste("Boxplot of",
var_name_1, "conditioned",
var_name_2),
xlab = var_name_2,
col = "lightcyan",
border = "black",
horizontal = FALSE # Boxplot orizzontale
)
}Anni.madre
Applichiamo le funzioni sulla variabile Anni.madre che rappresentà l’eta della madre ed è una variabile quantitativa discreta.
# visualizziamo l'istogramma
histogram(data$Anni.madre, 40, 0.08,
cex = 0.9,
ylim_max = 0.08)
# visualizziamo il boxplot
# Imposta una finestra grafica più grande
par(mar = c(5, 1, 4, 10)) # margini: (basso, sinistra, alto, destra)
boxplot_1(data$Anni.madre, 50, 1.5)
Dall’analisi visiva Anni.madre ci sembra distribuita normalmente e la mediana coincide quasi con la media.
Nel boxplot invece vediamo come sono presenti dei valori anomali soprattutto i più bassi, infatti il minimo di tale variabile è 0 ma risulta impossibile che una persona con 0 anni o comunque in età infantile possa avere un bambino.
Andiamo a visualizzare queste righe:
data[data$Anni.madre < 10, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1152 1 1 0 41 3250 490 350
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1152 Nat osp2 F
## 1380 Nat osp3 MDato che abbiamo solo 2 righe con i valori di Anni.madre che risultano per certo impossibili, andremo ad eliminare direttamente queste righe.
Creiamo il nuovo dataset filtrando tutte le righe con Anni.madre maggiore di 10:
data_2 = data[data$Anni.madre > 10, ]
row.names(data_2) = NULL # resetta gli indici e li fa ripartire da 1
n = nrow(data_2)
attach(data_2)Vediamo anche gli indici di asimmetria e curtosi.
Se questi sono vicini allo zero significa che la distribuzione può essere approssimata ad una normale.
moments::skewness(data_2$Anni.madre)## [1] 0.1510624moments::kurtosis(data_2$Anni.madre) - 3## [1] -0.1056061L’asimmetria e la curtosi sono vicine allo zero ad indicare che la distribuzione è simile ad una normale.
.
N_gravidanze
Indica il numero di gravidanze precedenti.
Questa è una variabile quantitativa discreta quindi sarà visualizzata meglio con un barplot e la relativa distribuzione di frequenza.
fi_gravidanze = table(data_2$N.gravidanze)
fr_gravidanze = round(table(data_2$N.gravidanze)/n, 3)
f_cum_gravidanze = cumsum(fi_gravidanze)
fr_cum_gravidanze = cumsum(fr_gravidanze)
f = cbind(fi_gravidanze, fr_gravidanze, f_cum_gravidanze, fr_cum_gravidanze)
data.frame(f)## fi_gravidanze fr_gravidanze f_cum_gravidanze fr_cum_gravidanze
## 0 1095 0.438 1095 0.438
## 1 817 0.327 1912 0.765
## 2 340 0.136 2252 0.901
## 3 150 0.060 2402 0.961
## 4 48 0.019 2450 0.980
## 5 21 0.008 2471 0.988
## 6 11 0.004 2482 0.992
## 7 1 0.000 2483 0.992
## 8 8 0.003 2491 0.995
## 9 2 0.001 2493 0.996
## 10 3 0.001 2496 0.997
## 11 1 0.000 2497 0.997
## 12 1 0.000 2498 0.997Grafichiamo le frequenze con un barplot ed un boxplot:
barplot(fr_gravidanze,
col = "skyblue",
main = "Distribuzione delle N.gravidanze",
xlab = "N.gravidanze",
ylab = "Frequenze relative",
ylim = c(0, 0.5),
border = "black") 
# creiamo il boxplot
par(mar = c(5.1, 2.1, 4.1, 8)) # c(5.1, 4.1, 4.1, 2.1), (basso, sinistra, alto, destra)
boxplot_1(data_2$N.gravidanze, 13, 1.5, 0.9)
La variabile N.gravidanze risulta sbilanciata in quanto la maggior parte dei dati ha 0,1,2 gravidanze, questo si riflette sul boxplot che indetifica i valori maggiori di 3 come outlier.
Dalla frequenza relativa comulata vediamo che circa il 90% dei nostri dati presenta N.gravidanze <= 2.
boxplot(data_2$Peso ~ data_2$N.gravidanze,
main = "Boxplot",
xlab = "N gravidanze",
ylab = "Peso alla nascita",
col = "lightcyan")
Essendoci pochissimi dati sopra le 1/2 gravidanze già avute, per semplificare e comprendere meglio la relazione tra il peso e le gravidanze, un’idea è raggruppare la variabile in due gruppi:
primigravide (1 gravidanza) vs. multigravide (2 o più gravidanze) ed effettuare un test.t.
# Creazione di una nuova variabile categoriale
data_2$Gravidanza_tipo <- ifelse(data_2$N.gravidanze <= 1, "Primigravida", "Multigravida")
# Converti in fattore (opzionale, ma utile per analisi successive)
data_2$Gravidanza_tipo <- factor(data_2$Gravidanza_tipo, levels = c("Primigravida", "Multigravida"))fi_gravidanze2 = table(data_2$Gravidanza_tipo)
fr_gravidanze2 = round(fi_gravidanze2/n, 2)
f_gravidanze2 = cbind(fi_gravidanze2, fr_gravidanze2)
data.frame(f_gravidanze2)## fi_gravidanze2 fr_gravidanze2
## Primigravida 1912 0.77
## Multigravida 586 0.23boxplot_condizionato(data_2$Peso, data_2$Gravidanza_tipo)
wilcox.test(Peso ~ Gravidanza_tipo, data = data_2)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Gravidanza_tipo
## W = 538433, p-value = 0.1539
## alternative hypothesis: true location shift is not equal to 0Concludiamo che non ci sono differenze significative in media tra i due gruppi.
.
Fumatrici
Fumatrici è una variabile binaria qualitativa che corrisponde a 0 se la madre non fuma ed 1 se fuma.
fr_fumatrici = table(data_2$Fumatrici)/n
round(fr_fumatrici, 3)##
## 0 1
## 0.958 0.042barplot(fr_fumatrici,
col = "skyblue",
main = "Distribuzione di Fumatrici",
xlab = "Fumatrici (0 = No, 1 = SI)",
ylab = "Frequenze relative",
ylim = c(0, 1),
border = "black") 
Vediamo come c’è una chiara e netta distinzione ad indicare che il dataset è fortemente sbilanciato verso le non fumatrici.
Se la variabile “Fumatrici” è codificata come:
1 = fumatrici
0 = non fumatrici
ma più 90% delle osservazioni corrispondono a non fumatrici, il modello potrebbe avere difficoltà a rilevare correttamente la relazione tra “essere fumatrici” e la variabile dipendente (peso), portando a stime inaffidabili.
In un modello di regressione il coefficiente di fumatrici potrebbe risultare non significativo a causa della scarsità di dati nel gruppo minoritario. Anche se esiste un effetto reale, il modello potrebbe non rilevarlo correttamente a causa della bassa potenza statistica.
Andiamo comunque a graficare con un boxplot la relazione tra il peso e le Fumatrici, ed effettuaimo un test t per verificare la differenza in media tra i due gruppi:
boxplot_condizionato(data_2$Peso, data_2$Fumatrici)
wilcox.test(Peso ~ Fumatrici, data = data_2)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Fumatrici
## W = 138069, p-value = 0.05928
## alternative hypothesis: true location shift is not equal to 0P-value > 0.05 concludiamo che per i dati che abbiamo non c’è una differenza significativa in media tra i due gruppi.
.
Gestazione
Indica il numero di settimane di gestazione del neonato.
Viene trattata come quantitativa discreta.
# creiamo l'istogramma
histogram(data_2$Gestazione, 25, 0.5, 1, 0.5)
# creiamo il boxplot
par(mar = c(5.1, 2.1, 4.1, 8)) # c(basso, sinistra, alto, destra)
# default (5.1, 4.1, 4.1, 2.1)
boxplot_1(data_2$Gestazione,
x_legend = 45,
y_legend = 1.5,
cex = 0.8)
è risaputo che la gestazione è di circa 9 mesi che equivalgono a circa 39 settimane.
Infatti la nostra mediana corrisponde a 39 settimane.
moments::skewness(data_2$Gestazione)## [1] -2.065131moments::kurtosis(data_2$Gestazione) - 3## [1] 8.255516A giudicare dall’istogramma e dagli indici di fisher e curtosi, vediamo che la distribuzione non è normale in quanto presenta una coda più allungata a sinistra e dei dati molto più concentrati intorno alla mediana.
L’asimmetria è data dai valori outlier trovati dal boxplot che risultano minori di 35 settimane di gestazione.
outliers_gestazione=data_2[data_2$Gestazione < 35, ]
nrow(outliers_gestazione)## [1] 67Ci sono 67 righe che hanno un valore di gestazione minore di 35.
boxplot(data_2$Peso ~ data_2$Gestazione,
main = "Boxplot",
xlab = "Settimane di Gestazione",
ylab = "Peso alla nascita",
col = "lightcyan")
Da un analisi condizionata tra il peso e le settimane di gestazione vediamo come ci sia un trend crescente tra le due variabili.
Questo ci suggerisce che Gestazione sarà una variabile significativa per prevedere il peso dei neonati.
Tuttavia essendoci solo 67 valori al di sotto di 35 settimane, su circa 2500 unità statistiche, probabilmente il modello farà fatica a catturare la relazione del peso con le settimane di gestazione più basse.
La relazione tra Gestazione e peso potrebbe essere anche non lineare in quanto tra 40, 41 e 42 settimane di gestazione il peso sembra essere abbastanza costante.
.
PESO
Rappresenta il peso dei neonati alla nascita ed è una variabile quantitativa continua espressa in grammi.
Questa è la nostra variabile target Y.
histogram(variable = data_2$Peso,
ylim_max = 0.0009,
y_legend = 0.0009,
x_legend = 700,
cex = 0.9)
par(mar = c(5.1, 2.1, 4.1, 9))
boxplot_1(variable = data_2$Peso,
x_legend = 5200,
cex = 0.9)
# Q-Q plot
qqnorm(data_2$Peso, main = "Q-Q Plot di Peso", xlab = "Quantili Normali", ylab = "Quantili dei Dati")
qqline(data_2$Peso, col = "red", lwd = 2)
moments::skewness(data_2$Peso)## [1] -0.6474036moments::kurtosis(data_2$Peso) - 3## [1] 2.028753shapiro.test(data_2$Peso)##
## Shapiro-Wilk normality test
##
## data: data_2$Peso
## W = 0.97068, p-value < 2.2e-16nrow(data_2[data_2$Peso < 2000, ])## [1] 50La distribuione è asimmetrica negativa e lo shapiro test ci suggerisce che non sia normale, anche il qq plot per i primi valori non si distribuiscono lungo la retta, probabilmente condizionati dagli outliers.
Se la variabile risposta non è normale, i residui del modello potrebbero non essere distribuiti normalmente, il che violerebbe una delle assunzioni chiave della regressione lineare.
La non normalità dei residui influisce sulla validità dei test statistici e sulle stime dei coefficienti, aumentando il rischio di falsi positivi o falsi negativi.
Per affrontare questo problema, possiamo:
Verificare i residui e applicare trasformazioni alla variabile target.
Utilizzare modelli robusti o modelli alternativi (come i GLM)
Dal boxplot condizionato peso-gestazione e dal boxplot solo del peso notiamo che i valori anomali minori di 35 di settimane di gestazione si trovano principalmente sotto i 2000 grammi di peso.
Lo stesso boxplot del peso indica i valori sotto 2000g come outilers, tuttavia abbiamo visto dal boxplot condizionato peso-gestazione che ci potrebbe essere un legame tra le due variabili.
potrebbe essere che il boxplot del peso classifichi sotto i 2000g il peso dei neonati come outliers poichè ci sono poche madri (unità statistiche) che hanno una gestazione sotto le 35 settimane circa.
Ma questi valori potrebbero risultare molto significativi se avessimo più unità statistiche di madri che hanno una gestazione inferiore a 35, poichè ci permetterebbero di saggiare con più certezza questa relazione crescente tra peso e gestazione.
.
Lunghezza
Rappresenta la lunghezza del neonato ed è una variabile quantitativa continua.
histogram(data_2$Lunghezza,
ylim_max = 0.025,
x_legend = 300,
y_legend = 0.025,
cex = 0.9)
par(mar = c(5.1, 2.1, 4.1, 9))
boxplot_1(data_2$Lunghezza,
x_legend = 580,
cex = 0.9)
Notiamo subito che la distribuzione della lunghezza è asimmetrica negativa.
La rapportiamo al peso con uno scatterplot:
plot(data_2$Lunghezza,
data_2$Peso)
Anche qui c’è un trend crescente abbastanza lineare e i punti sotto a 2000 di peso sono minori.
Ci aspetteremo che la Lunghezza sia una delle variabili significative per prevedere il peso dei neonati.
.
Cranio
Diametro carniale (quantitativa continua)
histogram(data$Cranio,
ylim_max = 0.035,
x_legend = 220,
y_legend = 0.035,
cex = 0.9)
#boxplot
par(mar = c(5.1, 2.1, 4.1, 9))
boxplot_1(data$Cranio,
x_legend = 400,
cex = 0.9)
#scatterplot
plot(data_2$Cranio,
data_2$Peso,
main = "Relazione tra Cranio e Peso",
xlab = "Dimensione del Cranio",
ylab = "Peso in g"
)
Anche qui vediamo come ci sia un trend crescente tra peso e grandezza del cranio.
Cranio sarà una variabile significativa per prevedere il peso dei neonati
.
Tipo parto
Variabile qualitativa binaria che rappresenta parto naturale (nat) o cesareo (ces).
Per tale variabile creeremo una distribuzione di frequenza e tracceremo un barplot.
fi_tipo_parto = table(data_2$Tipo.parto)
fr_tipo_parto = round(fi_tipo_parto/n, 2)
f_tipo_parto = cbind(fi_tipo_parto, fr_tipo_parto)
data.frame(f_tipo_parto)## fi_tipo_parto fr_tipo_parto
## Ces 728 0.29
## Nat 1770 0.71barplot(fr_tipo_parto,
col = "skyblue",
main = "Distribuzione del tipo parto",
xlab = "tipo parto (Ces, Nat)",
ylab = "Frequenze relative",
ylim = c(0, 1),
border = "black") 
boxplot_condizionato(data_2$Peso, data_2$Tipo.parto)
Tipo.parto risulta sbilanciata a favore dei parti naturali (71%) contro i cesari (29%).
Graficamente la mediana per peso e i due gruppi di Tipo.parto sono identiche.
Applichiamo un test un test t per saggiare l’ipotesi di differenza tra i gruppi.
car::leveneTest(Peso ~ Tipo.parto, data = data_2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 4.5678 0.03267 *
## 2496
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nel leneve test p< alpha l’ipotesi nulla viene rifiutata e concludiamo che le varianze non sono omogenee.
Per questo usiamo un test t non parametrico il Wilcoxon test:
wilcox.test(Peso ~ Tipo.parto, data = data_2)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Tipo.parto
## W = 633851, p-value = 0.5244
## alternative hypothesis: true location shift is not equal to 0p-value > 0.05, non c’è evidenza sufficiente per rifiutare l’ipotesi nulla (ovvero, non ci sono differenze significative tra i gruppi).
Ci aspetteremo quindi che Tipo.parto non sia rilevante al fine di predirre il peso.
.
Ospedale
Variabile Qualitativa che può assumere 3 modalità.
fi_ospedale = table(data_2$Ospedale)
fr_ospedale = round(fi_ospedale/n, 2)
f_ospedale = cbind(fi_ospedale, fr_ospedale)
data.frame(f_ospedale)## fi_ospedale fr_ospedale
## osp1 816 0.33
## osp2 848 0.34
## osp3 834 0.33barplot(fr_ospedale,
col = "skyblue",
main = "Distribuzione di Ospedale",
xlab = "Ospedale (osp1, osp2, osp3)",
ylab = "Frequenze relative",
ylim = c(0, 1),
border = "black") 
boxplot_condizionato(data_2$Peso, data_2$Ospedale)
Ospedale è una variabile bilanciata graficamente non sembraci essere differenze in media di peso tra i gruppi.
Saggiamo questa ipotesi anche con un Anova test.
ma per prima cosa saggiamo l’ipotesi di varianze omogenee.
car::leveneTest(Peso ~ Ospedale,
data = data_2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5912 0.5537
## 2495Per il test di levante p > alpha, concludiamo che le varianze sono omogenee.
Dato che la variabile Ospedale presenta 3 modalità utilizzeremo un Anova per confermare la non differenza significativa tra le medie dei gruppi.
anova = aov(Peso ~ Ospedale, data = data_2)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Ospedale 2 948538 474269 1.72 0.179
## Residuals 2495 687888603 275707Vediamo con un test di shapiro se i residui dell’anova si distribuiscono normalmente.
shapiro.test(residuals(anova))##
## Shapiro-Wilk normality test
##
## data: residuals(anova)
## W = 0.97071, p-value < 2.2e-16Il p-value è molto piccolo, questo ci suggerisce che i residui non sono distribuiti normalmente e quindi l’Anova test potrebbe non essere affidabile.
Effettueremo al suo posto un test non parametrico di kruskal-wallis.
kruskal.test(Peso ~ Ospedale, data = data_2)##
## Kruskal-Wallis rank sum test
##
## data: Peso by Ospedale
## Kruskal-Wallis chi-squared = 3.0609, df = 2, p-value = 0.2164Il p-value > 0.05 quindi concludiamo che non c’è una differenza significativa tra le medie dei due gruppi.
La variabile Ospedale non dovrebbe essere rilevante per prevedere il peso dei neonati.
.
Sesso
Sesso rappresenta il sesso del bambino ed è una variabile qualitativa binaria:
fi_sesso = table(data_2$Sesso)
fr_sesso = round(fi_sesso/n, 2)
f_sesso = cbind(fi_sesso, fr_sesso)
data.frame(f_sesso)## fi_sesso fr_sesso
## F 1255 0.5
## M 1243 0.5barplot(fr_sesso,
col = "skyblue",
main = "Distribuzione di Sesso",
xlab = "Sesso (M, F)",
ylab = "Frequenze relative",
ylim = c(0, 1),
border = "black") 
La variabile sesso risulta bilanciata.
boxplot_condizionato(data_2$Peso, data_2$Sesso)
Dal boxplot condizionato sembra esserci una differenza tra i gruppi delle due distribuzioni.
car::leveneTest(Peso ~ Sesso,
data = data_2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.8222 0.3646
## 2496p-value > alpha, le varianze sono omogenee.
Possiamo applicare il test t:
t.test(data_2$Peso ~ data_2$Sesso)##
## Welch Two Sample t-test
##
## data: data_2$Peso by data_2$Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496Per il test t p-value < alpha 0.05 quindi concludiamo che esiste una differenza significativa tra le medie dei due gruppi maschio e femmina.
Questo ci suggerisce che il Sesso sarà una delle variabili chiave per determinare il peso dei neonati.
.
Correlazioni variabili quantitative
# Caricare il pacchetto
library(corrplot)
numerical_data <- Filter(is.numeric, data_2)
numerical_data <- numerical_data[, !names(numerical_data) %in% "Fumatrici"]
cor_matrix <- cor(numerical_data)
# Calcola la matrice di correlazione usando il metodo di Spearman
cor_matrix_2 <- cor(numerical_data, method = "spearman")
par(mfrow = c(1, 2))
# Grafico della matrice di correlazione
# Grafico della matrice di correlazione con numeri e colori
corrplot(cor_matrix,
method = "color", # Quadrati colorati
type = "lower", # Mostra solo la parte superiore
addCoef.col = "black", # Aggiunge i numeri delle correlazioni in nero
tl.col = "black", # Colore delle etichette
tl.srt = 45, # Rotazione delle etichette
number.cex = 0.7) # Dimensione dei numeri
# Titolo per il primo grafico
title(main = "Correlazione Lineare", col.main = "blue", font.main = 2, cex.main = 1.2)
# Grafico della matrice di correlazione
corrplot(cor_matrix_2,
method = "color",
type = "lower",
addCoef.col = "black",
tl.col = "black",
tl.srt = 45,
number.cex = 0.7)
# Titolo per il secondo grafico
title(main = "Correlazione di Spearman", col.main = "red", font.main = 2, cex.main = 1.2)
.
Creazione del Modello di Regressione
Verrà sviluppato un modello di regressione lineare multipla con la procedura stepwise mista.
Questa procedura trova il miglior modello in automatico in base ad un criterio di selezione.
Come criterio useremo il BIC (Criterio di informazione bayesiano).
Il BIC bilancia la bontà di adattamento del modello (quanto bene il modello si adatta ai dati) e la complessità del modello (per evitare overfitting), il modello con il minor BIC sarà classificato come migliore.
Partiremo dal modello di regressione formato da tutte le variabili e applicheremo la procedura stepwise mista.
Dalla matrice di correlazione, dai boxplot condizionati e test abbiamo assunto che le variabili più rilevanti per determinare il peso dei neonati sono:
Lunghezza, cranio, gestazione (per le variabili quantitative)
Sesso (per le variabili qualitative)
Ci aspetteremo quindi che il modello migliore conservi queste variabili come esplicative.
Andiamo a costruire il primo modello di regressione con tutte le variabili:
mod_1 = lm(Peso ~ . -N.gravidanze, data = data_2)
stepwise.mod = MASS::stepAIC(mod_1,
direction = "both",
k = log(n))## Start: AIC=28114.74
## Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Gravidanza_tipo) -
## N.gravidanze
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 25688 186479901 28107
## - Ospedale 2 661287 187115500 28108
## - Fumatrici 1 92472 186546686 28108
## - Tipo.parto 1 466295 186920508 28113
## <none> 186454214 28115
## - Gravidanza_tipo 1 735224 187189438 28117
## - Sesso 1 3589476 190043690 28155
## - Gestazione 1 5418911 191873125 28179
## - Cranio 1 45634852 232089066 28654
## - Lunghezza 1 87826778 274280992 29071
##
## Step: AIC=28107.26
## Peso ~ Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + Gravidanza_tipo
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 665739 187145640 28101
## - Fumatrici 1 92791 186572692 28101
## - Tipo.parto 1 467425 186947327 28106
## <none> 186479901 28107
## - Gravidanza_tipo 1 931085 187410986 28112
## + Anni.madre 1 25688 186454214 28115
## - Sesso 1 3596273 190076174 28147
## - Gestazione 1 5405527 191885428 28171
## - Cranio 1 45972153 232452055 28650
## - Lunghezza 1 87804374 274284275 29063
##
## Step: AIC=28100.52
## Peso ~ Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso + Gravidanza_tipo
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 101914 187247555 28094
## - Tipo.parto 1 491756 187637397 28099
## <none> 187145640 28101
## - Gravidanza_tipo 1 1001742 188147382 28106
## + Ospedale 2 665739 186479901 28107
## + Anni.madre 1 30140 187115500 28108
## - Sesso 1 3637160 190782801 28141
## - Gestazione 1 5507430 192653070 28165
## - Cranio 1 46099531 233245172 28643
## - Lunghezza 1 87504213 274649853 29051
##
## Step: AIC=28094.05
## Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso +
## Gravidanza_tipo
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 483461 187731016 28093
## <none> 187247555 28094
## - Gravidanza_tipo 1 975746 188223301 28099
## + Fumatrici 1 101914 187145640 28101
## + Ospedale 2 674863 186572692 28101
## + Anni.madre 1 30556 187216999 28102
## - Sesso 1 3620162 190867717 28134
## - Gestazione 1 5437479 192685034 28158
## - Cranio 1 46135760 233383315 28636
## - Lunghezza 1 87936166 275183721 29048
##
## Step: AIC=28092.67
## Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Gravidanza_tipo
##
## Df Sum of Sq RSS AIC
## <none> 187731016 28093
## + Tipo.parto 1 483461 187247555 28094
## - Gravidanza_tipo 1 932091 188663107 28097
## + Ospedale 2 698769 187032247 28099
## + Fumatrici 1 93620 187637397 28099
## + Anni.madre 1 31896 187699120 28100
## - Sesso 1 3626871 191357887 28133
## - Gestazione 1 5458199 193189215 28156
## - Cranio 1 46492096 234223112 28638
## - Lunghezza 1 87503142 275234158 29041summary(stepwise.mod)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1171.65 -182.97 -15.48 163.32 2647.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6675.8431 135.4614 -49.282 < 2e-16 ***
## Gestazione 32.2513 3.7889 8.512 < 2e-16 ***
## Lunghezza 10.2362 0.3003 34.081 < 2e-16 ***
## Cranio 10.5568 0.4250 24.843 < 2e-16 ***
## SessoM 77.7251 11.2018 6.939 5.03e-12 ***
## Gravidanza_tipoMultigravida 45.8601 13.0377 3.518 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7275, Adjusted R-squared: 0.7269
## F-statistic: 1330 on 5 and 2492 DF, p-value: < 2.2e-16best_stepwise_mod = lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Gravidanza_tipo, data = data_2)
summary(best_stepwise_mod)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1171.65 -182.97 -15.48 163.32 2647.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6675.8431 135.4614 -49.282 < 2e-16 ***
## Gestazione 32.2513 3.7889 8.512 < 2e-16 ***
## Lunghezza 10.2362 0.3003 34.081 < 2e-16 ***
## Cranio 10.5568 0.4250 24.843 < 2e-16 ***
## SessoM 77.7251 11.2018 6.939 5.03e-12 ***
## Gravidanza_tipoMultigravida 45.8601 13.0377 3.518 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7275, Adjusted R-squared: 0.7269
## F-statistic: 1330 on 5 and 2492 DF, p-value: < 2.2e-16Tutti i coefficienti sono significativi e L’R2 aggiustato ci sta ad indicare che il modello spiega circa il 72% della variabilità del peso dei neonati.
Nel modello è stata inserita anche la variabile Gravidanza_tipo creata da noi dal raggruppamento della variabile N.gravidanze.
1. Grafico degli osservati vs predetti (per valutare la bontà del modello)
Questo grafico è utile per vedere come i valori predetti dal modello si confrontano con i valori osservati. Idealmente, i punti dovrebbero distribuire attorno alla linea di identità, che indica che i valori osservati e predetti sono uguali.
predetti_originali <- fitted(best_stepwise_mod)
plot(data_2$Peso, predetti_originali,
main = "Osservati vs Predetti",
xlab = "Osservati",
ylab = "Predetti",
pch = 20, col = "blue")
# Aggiungi la linea ideale (linea di identità)
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2) 
Notiamo che sembrerebbe esserci un pattern di non linearità in questo scatterplot e non tutti i punti si distriuiscono intorno alla linea di identità.
Dato che abbiamo ipotizzato esserci un effetto non lineare tra gestazione e peso, inseriamo questo effetto nel modello come funzione radice quadrata di gestazione.
mod_3 = lm(Peso ~ sqrt(Gestazione) + Lunghezza + Cranio + Sesso + Gravidanza_tipo, data = data_2)
summary(mod_3)##
## Call:
## lm(formula = Peso ~ sqrt(Gestazione) + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1173.40 -183.59 -15.48 163.06 2646.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7839.8149 231.5767 -33.854 < 2e-16 ***
## sqrt(Gestazione) 388.1491 46.2400 8.394 < 2e-16 ***
## Lunghezza 10.2380 0.3013 33.981 < 2e-16 ***
## Cranio 10.5495 0.4254 24.801 < 2e-16 ***
## SessoM 78.0797 11.2055 6.968 4.1e-12 ***
## Gravidanza_tipoMultigravida 45.7317 13.0425 3.506 0.000462 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared: 0.7273, Adjusted R-squared: 0.7267
## F-statistic: 1329 on 5 and 2492 DF, p-value: < 2.2e-16R2 è rimasto pressocchè identico quindi non aggiungeremo l’effetto non lineare al nostro modello.
Andiamo a vedere possibili collinearità tra le variabili del modello:
car::vif(best_stepwise_mod)## Gestazione Lunghezza Cranio Sesso Gravidanza_tipo
## 1.662111 2.072705 1.615678 1.040190 1.012058Sono tutti minori di 5 quindi non dovremmo avere problemi di multicollinearità tra le variabili.
.
Analisi dei residui
Ora andiamo ad analizzare i residui e possibili valori outliers e di leverage nel modello.
I valori di leva (leverage) sono dati che si trovano lontano dalla media delle variabili indipendenti in un modello di regressione, e possono influenzare in modo significativo i parametri stimati del modello.
Un valore di leva non è necessariamente un valore anomalo (cioè non è un valore estremo della variabile dipendente), ma può avere un grande impatto sui coefficienti della regressione, alterando la stima dei parametri.
Le assunzioni dei residui verranno testate con:
Test di omoschedasticità Breusch-Pagan
Scopo: Il test di Breusch-Pagan serve per verificare se vi è eteroscedasticità nei residui di un modello di regressione. L’eteroscedasticità si verifica quando la varianza dei residui non è costante lungo i valori predetti o delle variabili indipendenti.
Test di corellazione Darwin-Watson
Scopo: Il test di Durbin-Watson viene utilizzato per verificare la autocorrelazione dei residui. L’autocorrelazione nei residui è problematica, poiché indica che c’è una relazione tra i residui adiacenti nel modello di regressione, il che può invalidare le inferenze statistiche basate sui residui.
Test di normalità Shapiro-Wilk
Scopo: Il test di Shapiro-Wilk viene utilizzato per verificare la normalità dei residui. I residui dovrebbero essere normalmente distribuiti affinché le inferenze del modello siano valide, come nel caso dell’analisi di regressione lineare.
par(mfrow = c(2,2),
mar = c(2, 4.1, 2, 2.1))
plot(best_stepwise_mod)
Notiamo subito nel primo grafico in alto a destra che i residui non sembrano seguire un pattern evidente e sembrano distribuirsi attorno la media 0 eccetto per alcuni punti.
Anche per il q-q plot sembrano distribuirsi lungo la linea della bisettrice eccetto per i valori estremi.
Inoltre, sembrerebbe esserci un valore molto vicino alla soglia di leverage che potrebbe causare problemi.
Quello che farò è individuare possbili outliers e valori di leverage e valutare se togliere quei dati e riaddestrare il modello per ottenere risultati migliori.
Troviamo gli outliers:
plot(rstudent(best_stepwise_mod))
abline(h = c(-2, 2), col = "red", lty = 2, lwd=2)
Generalmente valori sopra o sotto le soglie di 2 e -2 sono indicati come outliers.
Svolgiamo un test t per individuare gli outlier:
car::outlierTest(best_stepwise_mod)## rstudent unadjusted p-value Bonferroni p
## 1549 10.089945 1.7189e-23 4.2938e-20
## 155 5.020133 5.5284e-07 1.3810e-03
## 1305 4.823687 1.4943e-06 3.7328e-03
## 1397 -4.290977 1.8468e-05 4.6133e-02La funzione outlierTest è utilizzata per indicare outlier estremi in un modello di regressione.
L’output della funzione si basa sui residui studentizzati e calcola un valore t per ogni osservazione, insieme al relativo p-value.
Ci identifica il numero di riga del possibile outlier.
Andiamo a visualizzare queste righe.
# Indici delle righe da selezionare
indici_outlier <- c(1549, 155, 1305, 1397)
# Filtrare il dataset per selezionare solo le righe con gli indici specificati
data_2[indici_outlier, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1549 35 1 0 38 4370 315 374
## 155 30 0 0 36 3610 410 330
## 1305 23 0 0 41 4900 510 352
## 1397 42 2 0 38 2560 525 349
## Tipo.parto Ospedale Sesso Gravidanza_tipo
## 1549 Nat osp3 F Primigravida
## 155 Nat osp1 M Primigravida
## 1305 Nat osp2 F Primigravida
## 1397 Ces osp2 M MultigravidaIndividuiamo i valori di leverage:
lev = hatvalues(best_stepwise_mod)
plot(lev)
p = sum(lev)
n = length(lev)
soglia = 2*p/n
abline(h=soglia, col=2, lwd=2)
Cook distance:
cook = cooks.distance(best_stepwise_mod)
cook[cook > 0.5]## 1549
## 0.840727L’unico valore che si trova sopra la soglia di avveritimento della distanza di Cook di 0.05 è quello ad indice 1549 che è classificato anche come outliers come abbiamo visto in precedenza.
Proverò a creare lo stesso modello togliendo i 4 valori outlier che abbiamo trovato incluso quello di leverage:
data_3 = data_2[- indici_outlier, ]
row.names(data_3) = NULL
n = nrow(data_3)mod_without_outliers <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Gravidanza_tipo, data = data_3)
summary(mod_without_outliers)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -968.05 -179.43 -13.78 162.28 1146.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6689.8531 131.0837 -51.035 < 2e-16 ***
## Gestazione 28.5942 3.6752 7.780 1.05e-14 ***
## Lunghezza 11.0438 0.2981 37.051 < 2e-16 ***
## Cranio 9.8323 0.4156 23.658 < 2e-16 ***
## SessoM 78.0066 10.8395 7.197 8.14e-13 ***
## Gravidanza_tipoMultigravida 52.4349 12.6178 4.156 3.35e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.4 on 2488 degrees of freedom
## Multiple R-squared: 0.744, Adjusted R-squared: 0.7435
## F-statistic: 1446 on 5 and 2488 DF, p-value: < 2.2e-16Quest’ultimo modello risulta migliore in quanto si sono abbassati gli errori standard, l’R2 è aumentato leggermente e per la variabile Gravidanza_tipo ora il coefficiente risulta più significativo.
.
verifichiamo i test di eteroschedasticità, correlazione e normalità:
library(lmtest)
bptest(mod_without_outliers)##
## studentized Breusch-Pagan test
##
## data: mod_without_outliers
## BP = 8.2522, df = 5, p-value = 0.1429dwtest(mod_without_outliers)##
## Durbin-Watson test
##
## data: mod_without_outliers
## DW = 1.9551, p-value = 0.1308
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(mod_without_outliers$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_without_outliers$residuals
## W = 0.99231, p-value = 3.18e-10Nel breusch-Pagan test il p-value è alto (> 0.05): Questo indica che l’eteroscedasticità non è presente. Ciò significa che la varianza dei residui è costante e va bene per il modello.
Nel Durbin-Watson test il p-value è alto (> 0.05): Questo significa che non c’è evidenza di autocorrelazione nei residui, il che è un buon segno per l’affidabilità del modello.
Nello shapiro-test il p-value è basso (< 0.05): Questo indica che i residui non seguono una distribuzione normale, il che potrebbe essere problematico. Dobbiamo tenere conto però che lo shapiro test è molto influenzato dalle code della distribuzione. In tal caso comunque, potremmo considerare trasformazioni dei dati o metodi robusti per trattare la non-normalità.
predetti_originali <- fitted(mod_without_outliers)
plot(data_3$Peso, predetti_originali,
main = "Osservati vs Predetti",
xlab = "Osservati",
ylab = "Predetti",
pch = 20, col = "blue")
# Aggiungi la linea ideale (linea di identità)
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2)
par(mfrow = c(2,2),
mar = c(2, 4.1, 2, 2.1))
plot(mod_without_outliers)
Andiamo ad applicare una trasformazione logaritmica alla variabile target peso poi confrontiamo il modello logaritmico con quest’ultimo precedente:
# Applicare la trasformazione logaritmica alla variabile Peso
mod_log <- lm(log(Peso) ~ Gestazione + Lunghezza + Cranio + Sesso + Gravidanza_tipo, data = data_3)
# Visualizzare il riepilogo del modello
summary(mod_log)##
## Call:
## lm(formula = log(Peso) ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40966 -0.05335 0.00069 0.05482 0.39047
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.396e+00 4.133e-02 106.366 < 2e-16 ***
## Gestazione 1.739e-02 1.159e-03 15.010 < 2e-16 ***
## Lunghezza 3.788e-03 9.398e-05 40.313 < 2e-16 ***
## Cranio 3.295e-03 1.310e-04 25.147 < 2e-16 ***
## SessoM 1.788e-02 3.418e-03 5.231 1.82e-07 ***
## Gravidanza_tipoMultigravida 1.666e-02 3.978e-03 4.188 2.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08366 on 2488 degrees of freedom
## Multiple R-squared: 0.7917, Adjusted R-squared: 0.7912
## F-statistic: 1891 on 5 and 2488 DF, p-value: < 2.2e-16Notiamo subito come l’R2 sia aumentato.
par(mfrow = c(2,2),
mar = c(2, 4.1, 2, 2.1))
plot(mod_log)
bptest(mod_without_outliers)##
## studentized Breusch-Pagan test
##
## data: mod_without_outliers
## BP = 8.2522, df = 5, p-value = 0.1429dwtest(mod_without_outliers)##
## Durbin-Watson test
##
## data: mod_without_outliers
## DW = 1.9551, p-value = 0.1308
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(mod_without_outliers$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_without_outliers$residuals
## W = 0.99231, p-value = 3.18e-10# Ottieni i residui del modello di regressione
residui <- residuals(mod_log)
# Crea l'istogramma dei residui
hist(residui,
main = "Istogramma dei residui",
xlab = "Residui",
col = "lightblue",
border = "black",
breaks = 20)
# Aggiungi una linea verticale sulla media dei residui
abline(v = mean(residui),
col = "blue", # Colore della linea
lwd = 2, # Spessore della linea
lty = 2) # Stile della linea
# Aggiungi la curva di densità
lines(density(residui), col = "red", lwd = 2)
Nel grafico qui sopra non riesco a capire perchè non mi tracca la curva di densità correttamente, attendo feedback per come posso modificare il codice.
.
Nonostante lo shapiro test sia ancora basso il grafico di distribuzione dei residui più gli altri relativi grafici ci suggeriscono che è possibile approssimare tale distribuzione ad una normale.
# Calcola i valori predetti nel dominio originale (inverso della trasformazione logaritmica)
predetti_originali <- exp(fitted(mod_log))
# Grafico degli osservati vs predetti (sui dati originali)
plot(data_3$Peso, predetti_originali,
main = "Osservati vs Predetti",
xlab = "Osservati",
ylab = "Predetti",
pch = 20, col = "blue")
# Aggiungi la linea ideale (linea di identità)
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2) # Linea di identità
I punti dei valori predetti ed osservati si distribuiscono meglio attorno alla linea di identità, che indica che i valori osservati e predetti sono uguali.
.
Andiamo a confrontare gli AIC i BIC dei modelli che abbiamo creato:
AIC(best_stepwise_mod,
mod_3,
mod_without_outliers,
mod_log)## df AIC
## best_stepwise_mod 7 35148.750
## mod_3 7 35150.690
## mod_without_outliers 7 34923.995
## mod_log 7 -5289.398BIC(best_stepwise_mod,
mod_3,
mod_without_outliers,
mod_log)## df BIC
## best_stepwise_mod 7 35189.513
## mod_3 7 35191.453
## mod_without_outliers 7 34964.747
## mod_log 7 -5248.646Il grafico dei valori predetti ed osservati, l’R2 ed il BIC più basso ci indica che il modello logaritmico è il migliore tra quelli testati.
summary(mod_log)##
## Call:
## lm(formula = log(Peso) ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40966 -0.05335 0.00069 0.05482 0.39047
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.396e+00 4.133e-02 106.366 < 2e-16 ***
## Gestazione 1.739e-02 1.159e-03 15.010 < 2e-16 ***
## Lunghezza 3.788e-03 9.398e-05 40.313 < 2e-16 ***
## Cranio 3.295e-03 1.310e-04 25.147 < 2e-16 ***
## SessoM 1.788e-02 3.418e-03 5.231 1.82e-07 ***
## Gravidanza_tipoMultigravida 1.666e-02 3.978e-03 4.188 2.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08366 on 2488 degrees of freedom
## Multiple R-squared: 0.7917, Adjusted R-squared: 0.7912
## F-statistic: 1891 on 5 and 2488 DF, p-value: < 2.2e-16La trasformazione logaritmica cambia l’interpretazione dei coefficienti da unità assolute a percentuali.
Variabili continue (Gestazione, Lunghezza, Cranio):
I coefficienti di queste variabili rappresentano la variazione percentuale del peso del neonato associata a un incremento di 1 unità nella variabile di riferimento:
Gestazione: Ogni settimana di gestazione in più è associata a un aumento del peso del 1.73%.
Lunghezza: Ogni centimetro in più di lunghezza è associato a un aumento del peso del 0.37%.
Cranio: Ogni unità in più nella misura del cranio è associata a un aumento del peso del 0.32%.
Variabili categoriche (Sesso, Gravidanza_tipo):
I coefficienti di queste variabili indicano la variazione percentuale del peso rispetto alla categoria di riferimento:
Sesso (Maschio): Essere maschio (rispetto al riferimento “femmina”) è associato a un peso 1.78% più alto.
Gravidanza_tipo (Multigravida): Avere avuto già più di una gravidanza è associato ad un peso del 1.66% più alto rispetto ad averne avuta 1 o nessuna (Primigravida).
Per avere un modello più semplice e di maggiore interpretazione andremo a testare anche il modello senza trasformazione e senza la variabile Gravidanza_tipo.
Inoltre, dato che Lunghezza e Cranio sono variabili abbastanza correlate tra loro proveremo a togliere anche una delle due.
mod_few_var = lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, data = data_3)
summary(mod_few_var)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -947.12 -181.01 -15.55 163.23 1130.97
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6662.0035 131.3395 -50.724 < 2e-16 ***
## Gestazione 27.5829 3.6791 7.497 9.01e-14 ***
## Lunghezza 10.9964 0.2988 36.799 < 2e-16 ***
## Cranio 9.9691 0.4156 23.985 < 2e-16 ***
## SessoM 79.5675 10.8683 7.321 3.30e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.2 on 2489 degrees of freedom
## Multiple R-squared: 0.7423, Adjusted R-squared: 0.7418
## F-statistic: 1792 on 4 and 2489 DF, p-value: < 2.2e-16mod_few_var_2 = lm(formula = Peso ~ Gestazione + Lunghezza + Sesso, data = data_3)
summary(mod_few_var_2)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Sesso, data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1153.80 -190.61 -22.48 183.29 1449.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5360.9312 132.6922 -40.401 < 2e-16 ***
## Gestazione 38.6165 4.0493 9.537 < 2e-16 ***
## Lunghezza 14.3380 0.2933 48.892 < 2e-16 ***
## SessoM 89.7904 12.0473 7.453 1.25e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 295.3 on 2490 degrees of freedom
## Multiple R-squared: 0.6827, Adjusted R-squared: 0.6823
## F-statistic: 1786 on 3 and 2490 DF, p-value: < 2.2e-16summary(best_stepwise_mod)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gravidanza_tipo, data = data_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1171.65 -182.97 -15.48 163.32 2647.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6675.8431 135.4614 -49.282 < 2e-16 ***
## Gestazione 32.2513 3.7889 8.512 < 2e-16 ***
## Lunghezza 10.2362 0.3003 34.081 < 2e-16 ***
## Cranio 10.5568 0.4250 24.843 < 2e-16 ***
## SessoM 77.7251 11.2018 6.939 5.03e-12 ***
## Gravidanza_tipoMultigravida 45.8601 13.0377 3.518 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7275, Adjusted R-squared: 0.7269
## F-statistic: 1330 on 5 and 2492 DF, p-value: < 2.2e-16AIC(best_stepwise_mod,
mod_3,
mod_without_outliers,
mod_log,
mod_few_var,
mod_few_var_2)## df AIC
## best_stepwise_mod 7 35148.750
## mod_3 7 35150.690
## mod_without_outliers 7 34923.995
## mod_log 7 -5289.398
## mod_few_var 6 34939.246
## mod_few_var_2 5 35455.821BIC(best_stepwise_mod,
mod_3,
mod_without_outliers,
mod_log,
mod_few_var,
mod_few_var_2)## df BIC
## best_stepwise_mod 7 35189.513
## mod_3 7 35191.453
## mod_without_outliers 7 34964.747
## mod_log 7 -5248.646
## mod_few_var 6 34974.176
## mod_few_var_2 5 35484.930par(mfrow = c(2,2),
mar = c(2, 4.1, 2, 2.1))
plot(mod_few_var)
# Ottieni i residui del modello di regressione
residui <- residuals(mod_few_var)
# Crea l'istogramma dei residui
hist(residui,
main = "Istogramma dei residui",
xlab = "Residui",
col = "lightblue",
border = "black",
breaks = 20)
# Aggiungi una linea verticale sulla media dei residui
abline(v = mean(residui),
col = "blue", # Colore della linea
lwd = 2, # Spessore della linea
lty = 2) # Stile della linea
# Aggiungi la curva di densità
lines(density(residui), col = "red", lwd = 2)
# Calcola i valori predetti nel dominio originale (inverso della trasformazione logaritmica)
predetti_originali <- fitted(mod_few_var)
# Grafico degli osservati vs predetti (sui dati originali)
plot(data_3$Peso, predetti_originali,
main = "Osservati vs Predetti",
xlab = "Osservati",
ylab = "Predetti",
pch = 20, col = "blue")
# Aggiungi la linea ideale (linea di identità)
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2) # Linea di identità
Calcoliamo RMSE
Il Root Mean Squared Error (MSE) è una metrica utilizzata per valutare la qualità di un modello di regressione.
residui <- data_3$Peso - predict(mod_few_var, newdata = data_3)
# Calcolo del Mean Squared Error (MSE)
mse <- mean(residui^2)
# Calcolo del Root Mean Squared Error (RMSE)
rmse <- sqrt(mse)
# Stampa il risultato
cat("Il Root Mean Squared Error (RMSE) del modello è:", rmse, "\n")## Il Root Mean Squared Error (RMSE) del modello è: 265.9521In media il modello sbaglia di circa 266 g.
Useremo il nostro modello “mod_few_var” per fare le nostre predizioni poichè, dal BIC, R2 e gli altri dati trovati, ha un rapporto efficienza/interpretazione migliore rispetto agli altri.
mod_few_var = lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, data = data_3)
summary(mod_few_var)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -947.12 -181.01 -15.55 163.23 1130.97
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6662.0035 131.3395 -50.724 < 2e-16 ***
## Gestazione 27.5829 3.6791 7.497 9.01e-14 ***
## Lunghezza 10.9964 0.2988 36.799 < 2e-16 ***
## Cranio 9.9691 0.4156 23.985 < 2e-16 ***
## SessoM 79.5675 10.8683 7.321 3.30e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.2 on 2489 degrees of freedom
## Multiple R-squared: 0.7423, Adjusted R-squared: 0.7418
## F-statistic: 1792 on 4 and 2489 DF, p-value: < 2.2e-16Interpretazione
Tenendo costanti le altre variabili tranne quella di interesse:
Per ogni settimana in più di gestazione il peso aumenta in media di 27.58g.
Per ogni centimetro in più di lunghezza, il peso aumenta in media di circa 11g
Per ogni centimetro in più del cranio il peso aumenta in media di circa 10g
l coefficiente associato alla variabile Sesso (maschio) signifca che i maschi a parità di tutte le altre variabili pesano in media circa 80g in più rispetto alle femmine. L’interpretazione del coefficiente è valida solo a parità di tutte le altre variabili incluse nel modello. Questo è un aspetto cruciale da ricordare, poiché l’effetto di Sesso non è “isolato”, ma è condizionato dai valori di Gestazione, Lunghezza, e Cranio.
Il modello spiega circa il 74% della variabilità del peso dei neonati.
mod_few_var_3 = lm(formula = Peso ~ Gestazione + Sesso + Lunghezza, data = data_3)
summary(mod_few_var_3)##
## Call:
## lm(formula = Peso ~ Gestazione + Sesso + Lunghezza, data = data_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1153.80 -190.61 -22.48 183.29 1449.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5360.9312 132.6922 -40.401 < 2e-16 ***
## Gestazione 38.6165 4.0493 9.537 < 2e-16 ***
## SessoM 89.7904 12.0473 7.453 1.25e-13 ***
## Lunghezza 14.3380 0.2933 48.892 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 295.3 on 2490 degrees of freedom
## Multiple R-squared: 0.6827, Adjusted R-squared: 0.6823
## F-statistic: 1786 on 3 and 2490 DF, p-value: < 2.2e-16Calcolare il peso stimato
# Creiamo un nuovo dataset con i valori delle variabili da utilizzare per la previsione
new_data <- data.frame(
Gestazione = 39, # Settimana di gestazione
Lunghezza = 480, # Lunghezza
Cranio = 330, # Circonferenza cranio
Sesso = "F" # Sesso (nome esattamente come nella tua variabile originale)
)
# Calcoliamo il peso stimato utilizzando la funzione predict
peso_stimato <- predict(mod_few_var, newdata = new_data)
# Stampiamo il risultato
cat("Il peso stimato della neonata è:", peso_stimato, "unità\n")## Il peso stimato della neonata è: 2981.83 unità.
Visualizzazioni
Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del sesso sul peso previsto.
# Carichiamo i pacchetti necessari
library(ggplot2)
library(dplyr)
# Creiamo un dataset predittivo con varie combinazioni di Gestazione e Sesso
gestazione_values <- seq(30, 42, by = 1) # Valori da 30 a 42 settimane
sesso_values <- c("F", "M") # Categorie di Sesso
# Generiamo un nuovo dataset con tutte le combinazioni
new_data <- expand.grid(
Gestazione = gestazione_values,
Lunghezza = mean(data_3$Lunghezza, na.rm = TRUE), # Lunghezza media
Cranio = mean(data_3$Cranio, na.rm = TRUE), # Cranio medio
Sesso = sesso_values
)
# Prediciamo il peso con il modello
new_data$Peso_Previsto <- predict(mod_few_var, newdata = new_data)
# Visualizzazione 1: Relazione tra Gestazione, Sesso e Peso Previsto
ggplot(new_data, aes(x = Gestazione, y = Peso_Previsto, color = Sesso)) +
geom_line(size = 1.2) +
labs(
title = "Impatto della Gestazione e del Sesso sul Peso Previsto",
x = "Settimane di Gestazione",
y = "Peso Previsto",
color = "Sesso"
) +
theme_minimal()
# Visualizzazione 2: Boxplot per l'effetto del Sesso sul Peso (mediante dati osservati)
ggplot(data_3, aes(x = Sesso, y = Peso, fill = Sesso)) +
geom_boxplot() +
labs(
title = "Distribuzione del Peso per Sesso",
x = "Sesso",
y = "Peso"
) +
theme_minimal() +
scale_fill_manual(values = c("F" = "pink", "M" = "lightcyan"))
Graficamente vediamo come con l’aumentare delle settimane di gestazione il peso previsto tende ad aumentare.
è chiara la differenza tra maschi e femmine dove per le femmine c’è una retta più bassa ad indicare che il peso previsto delle femmine è inferiore rispetto a quello dei maschi a parità di tutte le altre variabili.