Modello Statistico per la Previsione del Peso Neonatale
Carlo Calì
2025-03-08
Libraries
library(dplyr)##
## Caricamento pacchetto: 'dplyr'## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(summarytools)
library(readr)
library(moments)
library(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.4.3## Caricamento del pacchetto richiesto: carData## Warning: il pacchetto 'carData' è stato creato con R versione 4.4.3##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recodelibrary(corrplot)## Warning: il pacchetto 'corrplot' è stato creato con R versione 4.4.3## corrplot 0.95 loadedlibrary(MASS)## Warning: il pacchetto 'MASS' è stato creato con R versione 4.4.3##
## Caricamento pacchetto: 'MASS'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## selectLoad Data
Neonati_df <- read_csv("neonati.csv", show_col_types = FALSE)0) Primo sguardo ai dati
Neonati_df## # A tibble: 2,500 × 10
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 26 0 0 42 3380 490 325
## 2 21 2 0 39 3150 490 345
## 3 34 3 0 38 3640 500 375
## 4 28 1 0 41 3690 515 365
## 5 20 0 0 38 3700 480 335
## 6 32 0 0 40 3200 495 340
## 7 26 1 0 39 3100 480 345
## 8 25 0 0 40 3580 510 349
## 9 22 1 0 40 3670 500 335
## 10 23 0 0 41 3700 510 362
## # ℹ 2,490 more rows
## # ℹ 3 more variables: Tipo.parto <chr>, Ospedale <chr>, Sesso <chr>Verifico che non ci siano valori mancanti e duplicati
colSums(is.na(Neonati_df))## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0sum(duplicated(Neonati_df))## [1] 0Non sono presenti né valori mancanti né duplicati all’interno del dataset
1) Analisi Descrittiva
# Suddivisione delle variabili in base alla loro tipologia
bin_var <- c("Fumatrici", "Tipo.parto", "Sesso")
cat_var <- c("Ospedale")
num_var <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")# Statistica descrittiva delle variabili binarie
for (temp_var in bin_var){
print(temp_var)
ni <- table(Neonati_df[[temp_var]])
fi <- round(100*prop.table(ni), 1)
print(cbind(ni, fi))
barplot(table(Neonati_df[[temp_var]]), col = c("red", "blue"), main = paste("Distribuzione", temp_var))
}## [1] "Fumatrici"
## ni fi
## 0 2396 95.8
## 1 104 4.2
## [1] "Tipo.parto"
## ni fi
## Ces 728 29.1
## Nat 1772 70.9
## [1] "Sesso"
## ni fi
## F 1256 50.2
## M 1244 49.8
è preponderante la presenza di madri non fumatrici rispetto a quelle fumatrici. Quasi il 30% dei parti presenti all’interno del dataset è cesario. Il sesso dei nascituri è bilanciato, con valori prossimi al 50% per entrambi i valori
# Statistica descrittiva delle variabili categoriali
for (temp_var in cat_var){
print(temp_var)
ni <- table(Neonati_df[[temp_var]])
fi <- round(100*prop.table(ni), 1)
print(cbind(ni, fi))
barplot(table(Neonati_df[[temp_var]]), col = c("red", "blue", "green"), main = paste("Distribuzione", temp_var))
}## [1] "Ospedale"
## ni fi
## osp1 816 32.6
## osp2 849 34.0
## osp3 835 33.4
I dati resgistrati provengono in modo bilanciato da tutte e tre gli ospedali presi in esame.
# Statistica descrittiva delle variabili quantitative
num_var_stats <- data.frame(
Variable = num_var,
Min = round(sapply(Neonati_df[num_var], min, na.rm = TRUE), 1),
Q1 = round(sapply(Neonati_df[num_var], function(x) quantile(x, probs = 0.25, na.rm = TRUE)), 1),
Q3 = round(sapply(Neonati_df[num_var], function(x) quantile(x, probs = 0.75, na.rm = TRUE)), 1),
Max = round(sapply(Neonati_df[num_var], max, na.rm = TRUE), 1),
Range = round(sapply(Neonati_df[num_var], function(x) max(x, na.rm = TRUE) - min(x, na.rm = TRUE)), 1),
IQR = round(sapply(Neonati_df[num_var], IQR, na.rm = TRUE), 1),
Mean = round(sapply(Neonati_df[num_var], mean, na.rm = TRUE), 2),
Median = round(sapply(Neonati_df[num_var], median, na.rm = TRUE), 2),
Variance = round(sapply(Neonati_df[num_var], var, na.rm = TRUE), 3),
Std_Dev = round(sapply(Neonati_df[num_var], sd, na.rm = TRUE), 3),
RSD = round(sapply(Neonati_df[num_var], function(x) (sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)) * 100), 3),
Skewness = round(sapply(Neonati_df[num_var], skewness, na.rm = TRUE), 3),
Kurtosis = round(sapply(Neonati_df[num_var], kurtosis, na.rm = TRUE), 3)
)
print(num_var_stats)## Variable Min Q1 Q3 Max Range IQR Mean Median
## Anni.madre Anni.madre 0 25 32 46 46 7 28.16 28
## N.gravidanze N.gravidanze 0 0 1 12 12 1 0.98 1
## Gestazione Gestazione 25 38 40 43 18 2 38.98 39
## Peso Peso 830 2990 3620 4930 4100 630 3284.08 3300
## Lunghezza Lunghezza 310 480 510 565 255 30 494.69 500
## Cranio Cranio 235 330 350 390 155 20 340.03 340
## Variance Std_Dev RSD Skewness Kurtosis
## Anni.madre 27.811 5.274 18.725 0.043 3.380
## N.gravidanze 1.640 1.281 130.512 2.514 13.989
## Gestazione 3.492 1.869 4.794 -2.065 11.258
## Peso 275665.683 525.039 15.987 -0.647 5.032
## Lunghezza 692.671 26.319 5.320 -1.515 9.487
## Cranio 269.791 16.425 4.831 -0.785 5.946for (temp_var in num_var){
print(temp_var)
# Recupero dei valori di Q1 e IQR dalla tabella `num_var_stats`
Q1 <- num_var_stats[num_var_stats$Variable == temp_var, "Q1"]
Q3 <- num_var_stats[num_var_stats$Variable == temp_var, "Q3"]
IQR <- num_var_stats[num_var_stats$Variable == temp_var, "IQR"]
# Calcolo del numero di outlier
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
N_Outliers <- sum(Neonati_df[[temp_var]] < lower_bound | Neonati_df[[temp_var]] > upper_bound, na.rm = TRUE)
print(paste("Numero di Outlier per", temp_var, ":", N_Outliers))
boxplot(Neonati_df[[temp_var]], col = "orange", main = paste("Distribuzione", temp_var))
}## [1] "Anni.madre"
## [1] "Numero di Outlier per Anni.madre : 13"
## [1] "N.gravidanze"
## [1] "Numero di Outlier per N.gravidanze : 246"
## [1] "Gestazione"
## [1] "Numero di Outlier per Gestazione : 67"
## [1] "Peso"
## [1] "Numero di Outlier per Peso : 69"
## [1] "Lunghezza"
## [1] "Numero di Outlier per Lunghezza : 59"
## [1] "Cranio"
## [1] "Numero di Outlier per Cranio : 48"
Ad esclusione dell’età della madre, tutte le variabili quantitative presentano distribuzioni non normali e asimmetriche. Tuttavia bisogna far presente che dal grafico della variabile Anni madre ci sono due valori insolitamente bassi, di cui uno di essi è 0. Questo potrebbe essere dato da qualche errore nella trasposizione dei dati. Si può optare o di rimuovere queste osservazioni perdendo però anche le informazioni sulle altre variabili o di sostituiire il valore errato con la media della distribuzione, questo ci permette di mantenere le altre informazioni e non dovrebbe comportare
# Ispezione delle osservazioni anomale
Neonati_df[Neonati_df$Anni.madre<10, ]## # A tibble: 2 × 10
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 1 0 41 3250 490 350 Nat
## 2 0 0 0 39 3060 490 330 Nat
## # ℹ 2 more variables: Ospedale <chr>, Sesso <chr>Le due osservazioni anomale non presentano altre stranezze negli altri parametri, er tale motivo si opterà per la sostituizione dei valori anomali con il valore medio.
Neonati_df$Anni.madre[Neonati_df$Anni.madre < 10] <- 28Adesso si può ricalcolare e graficare nuovamente le variabili quantitiative.
# Statistica descrittiva delle variabili quantitative
num_var_stats <- data.frame(
Variable = num_var,
Min = round(sapply(Neonati_df[num_var], min, na.rm = TRUE), 1),
Q1 = round(sapply(Neonati_df[num_var], function(x) quantile(x, probs = 0.25, na.rm = TRUE)), 1),
Q3 = round(sapply(Neonati_df[num_var], function(x) quantile(x, probs = 0.75, na.rm = TRUE)), 1),
Max = round(sapply(Neonati_df[num_var], max, na.rm = TRUE), 1),
Range = round(sapply(Neonati_df[num_var], function(x) max(x, na.rm = TRUE) - min(x, na.rm = TRUE)), 1),
IQR = round(sapply(Neonati_df[num_var], IQR, na.rm = TRUE), 1),
Mean = round(sapply(Neonati_df[num_var], mean, na.rm = TRUE), 2),
Median = round(sapply(Neonati_df[num_var], median, na.rm = TRUE), 2),
Variance = round(sapply(Neonati_df[num_var], var, na.rm = TRUE), 3),
Std_Dev = round(sapply(Neonati_df[num_var], sd, na.rm = TRUE), 3),
RSD = round(sapply(Neonati_df[num_var], function(x) (sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)) * 100), 3),
Skewness = round(sapply(Neonati_df[num_var], skewness, na.rm = TRUE), 3),
Kurtosis = round(sapply(Neonati_df[num_var], kurtosis, na.rm = TRUE), 3)
)
print(num_var_stats)## Variable Min Q1 Q3 Max Range IQR Mean Median
## Anni.madre Anni.madre 13 25 32 46 33 7 28.19 28
## N.gravidanze N.gravidanze 0 0 1 12 12 1 0.98 1
## Gestazione Gestazione 25 38 40 43 18 2 38.98 39
## Peso Peso 830 2990 3620 4930 4100 630 3284.08 3300
## Lunghezza Lunghezza 310 480 510 565 255 30 494.69 500
## Cranio Cranio 235 330 350 390 155 20 340.03 340
## Variance Std_Dev RSD Skewness Kurtosis
## Anni.madre 27.197 5.215 18.503 0.151 2.897
## N.gravidanze 1.640 1.281 130.512 2.514 13.989
## Gestazione 3.492 1.869 4.794 -2.065 11.258
## Peso 275665.683 525.039 15.987 -0.647 5.032
## Lunghezza 692.671 26.319 5.320 -1.515 9.487
## Cranio 269.791 16.425 4.831 -0.785 5.946for (temp_var in num_var){
print(temp_var)
# Recupero dei valori di Q1 e IQR dalla tabella `num_var_stats`
Q1 <- num_var_stats[num_var_stats$Variable == temp_var, "Q1"]
Q3 <- num_var_stats[num_var_stats$Variable == temp_var, "Q3"]
IQR <- num_var_stats[num_var_stats$Variable == temp_var, "IQR"]
# Calcolo del numero di outlier
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
N_Outliers <- sum(Neonati_df[[temp_var]] < lower_bound | Neonati_df[[temp_var]] > upper_bound, na.rm = TRUE)
print(paste("Numero di Outlier per", temp_var, ":", N_Outliers))
boxplot(Neonati_df[[temp_var]], col = "orange", main = paste("Distribuzione", temp_var))
}## [1] "Anni.madre"
## [1] "Numero di Outlier per Anni.madre : 11"
## [1] "N.gravidanze"
## [1] "Numero di Outlier per N.gravidanze : 246"
## [1] "Gestazione"
## [1] "Numero di Outlier per Gestazione : 67"
## [1] "Peso"
## [1] "Numero di Outlier per Peso : 69"
## [1] "Lunghezza"
## [1] "Numero di Outlier per Lunghezza : 59"
## [1] "Cranio"
## [1] "Numero di Outlier per Cranio : 48"
Valutiamo adesso tramite test di ipotesi se ci sono ospedali in cui si fanno un maggior numero di parti cesari. Dato che il confronto verte su due tipologie di variabili categoriali verrà eseguito un Chi-square test.
preg_hosp_table <- table(Neonati_df$Tipo.parto, Neonati_df$Ospedale)
print(preg_hosp_table)##
## osp1 osp2 osp3
## Ces 242 254 232
## Nat 574 595 603print(prop.table(preg_hosp_table))##
## osp1 osp2 osp3
## Ces 0.0968 0.1016 0.0928
## Nat 0.2296 0.2380 0.2412chisq.test(preg_hosp_table)##
## Pearson's Chi-squared test
##
## data: preg_hosp_table
## X-squared = 1.0972, df = 2, p-value = 0.5778Dato che dal test il p-value è maggiore di 0.05 non siamo in grado di rifiutare l’ipotesi nulla e quindi accettiamo che non ci sono prove evidenti che un ospedale faccia più parti cesari rispetto ad altri.
Si vuole anche valutare se i dati del peso e della lunghezza raccolti dal campione sono significativamente uguali a quelli della popolazione. Per fare ciò si svolgerà un one-sample t-test per ciascuna caratteristica in cui impongo come valore di confronto per la media della popolazione il valore riportato dal sito dell’ospedale Bambino Gesù.
# One-Sample t-Test per il peso
t.test(Neonati_df$Peso, mu=3300)##
## One Sample t-test
##
## data: Neonati_df$Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081# One-Sample t-Test per la lunghezza
t.test(Neonati_df$Lunghezza, mu=500)##
## One Sample t-test
##
## data: Neonati_df$Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692Dai dati ottenuti si evince che nel caso del peso il p-value è maggiore di 0.05 e quindi non possiamo rifiutare l’ipotesi nulla, ciò significa che il nostro campione riporta una media del peso dei neonati simile a quella della popolazione, ciò è anche confermato dal fatto che all’interno dell’intervallo di confidenza ricade il valore medio atteso. Viceversa per quanto riguarda la lunghezza è stato ottenuto un p-value prossimo allo 0, ciò significa che dobbiamo rifiutare l’ipotesi di uguaglianza e accettare che il nostro campione non riporta una media delle lunghezze dei neonati significativamente uguale a quella della popolazione, ciò è avvalorato dal fatto che l’intervallo di confidenza non contiene il valore medio della popolazione.
Infine come ultima valutazione di confronto tra i dati raccolti all’interno del dataset si vuole verificare se ci sono significative differenze tra i sessi per le misure antropometriche(Peso, Lugnhezza, Dimensione Cranio)
Prima di eseguire il test di confronto bisogna valutare normalità di ogni singola distribuzione tramite il test di Shapiro-Wilk, se ciò non fosse rispettata i gruppi verrano confrontati tramite Wilcoxon test. Viceversa, se i gruppi hanno una distribuzione normale si passerà al test di Levene per verificare che l’uguagluanzia tra le varianze. Qualora le varianze tra i gruppi sono uguali per la medesima proprietà allora si potrà eseguire il confronto tramite il t-test Student, altrimenti si eseguirà il t-test di Welch.
shapiro.test(Neonati_df$Peso[Neonati_df$Sesso == "M"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Peso[Neonati_df$Sesso == "M"]
## W = 0.96647, p-value = 2.321e-16shapiro.test(Neonati_df$Peso[Neonati_df$Sesso == "F"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Peso[Neonati_df$Sesso == "F"]
## W = 0.96285, p-value < 2.2e-16Le distribuzioni del peso per entrambi i sessi non sono normali quindi si eseguirà un Wilcoxon Test
wilcox.test(Peso ~ Sesso, data = Neonati_df)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Sesso
## W = 538641, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Il test per il confronto dei pesi tra i sessi ha un p-value molto più piccolo di 0.05. Per tale ragione dobbiamo rifiutare l’ipotesi nulla di uguaglianza e affermare che c’è una differenza significativa nel peso dei neonati tra maschi e femmine.
median(Neonati_df$Peso[Neonati_df$Sesso == "M"])## [1] 3430median(Neonati_df$Peso[Neonati_df$Sesso == "F"])## [1] 3160Procediamo in modo analogo con le altre caratteristiche.
shapiro.test(Neonati_df$Lunghezza[Neonati_df$Sesso == "M"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Lunghezza[Neonati_df$Sesso == "M"]
## W = 0.92028, p-value < 2.2e-16shapiro.test(Neonati_df$Lunghezza[Neonati_df$Sesso == "F"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Lunghezza[Neonati_df$Sesso == "F"]
## W = 0.89953, p-value < 2.2e-16Anche in questo caso le distribuzioni delle lunghezze per entrambi i sessi non sono normali quindi si eseguirà un Wilcoxon Test
wilcox.test(Lunghezza ~ Sesso, data = Neonati_df)##
## Wilcoxon rank sum test with continuity correction
##
## data: Lunghezza by Sesso
## W = 594455, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Il test per il confronto delle lunghezze tra i sessi ha un p-value molto più piccolo di 0.05. Per tale ragione dobbiamo rifiutare l’ipotesi nulla di uguaglianza e affermare che c’è una differenza significativa nelle lunghezze dei neonati tra maschi e femmine.
median(Neonati_df$Lunghezza[Neonati_df$Sesso == "M"])## [1] 500median(Neonati_df$Lunghezza[Neonati_df$Sesso == "F"])## [1] 490shapiro.test(Neonati_df$Cranio[Neonati_df$Sesso == "M"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Cranio[Neonati_df$Sesso == "M"]
## W = 0.97046, p-value = 3.006e-15shapiro.test(Neonati_df$Cranio[Neonati_df$Sesso == "F"])##
## Shapiro-Wilk normality test
##
## data: Neonati_df$Cranio[Neonati_df$Sesso == "F"]
## W = 0.95543, p-value < 2.2e-16Anche in quest’ultimo caso le distribuzioni delle dimensioni del cranio per entrambi i sessi non sono normali quindi si eseguirà un Wilcoxon Test
wilcox.test(Cranio ~ Sesso, data = Neonati_df)##
## Wilcoxon rank sum test with continuity correction
##
## data: Cranio by Sesso
## W = 641638, p-value = 9.633e-15
## alternative hypothesis: true location shift is not equal to 0Il test per il confronto delle dimensioni del cranio tra i sessi ha un p-value molto più piccolo di 0.05. Per tale ragione dobbiamo rifiutare l’ipotesi nulla di uguaglianza e affermare che c’è una differenza significativa nelle dimensioni dei neonati tra maschi e femmine.
median(Neonati_df$Cranio[Neonati_df$Sesso == "M"])## [1] 343median(Neonati_df$Cranio[Neonati_df$Sesso == "F"])## [1] 340Alla luce di questi test si evince che per tutte e tre le misure antropometriche ci sono delle differenze significative tra maschi e femmine. In generale i maschi tendono ad avere valori superiori alle femmine.
2) Creazione del modello
Si vuole sviluppare un modello di regressione lineare multipla per prevedere il peso del neonato in base alle variabili rilevanti.
2.1) Selezione delle varaibili
Innanzitutto bisogna valutare e quantificare l’impatto di ogni variabile sul fattore peso del neonato.
Come primo approcio si può calcolare la matrice di correlazione tra le variabili quantitative. In questa circostanza ci è utile codificare le variabili categoriali binarie in valori numerici 0 e 1, senza però tralascaire il fatto che per il calcolo della correlazioni tra valori continui e binari non è un approccio ottimale. Così facendo l’unica varaibile esclusa sarà l’ospedale presso cui è stata eseguita l’osservazione, che verrà valutata per completezza a parte.
#Conversione dei valori binari in numerici
Neonati_df$Tipo.parto_num <- ifelse(Neonati_df$Tipo.parto == "Ces", 1, 0)
Neonati_df$Sesso_num <- ifelse(Neonati_df$Sesso == "M", 1, 0)
#Matrice di correlazione delle varaibili
num_fact <- Neonati_df[, c("Anni.madre", "N.gravidanze", "Fumatrici",
"Gestazione", "Lunghezza", "Cranio", "Tipo.parto_num",
"Sesso_num", "Peso")]
corr_matr <- cor(num_fact, use = "complete.obs")
#Heatmap della matrice di correlazione
corrplot(corr_matr, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", tl.srt = 45)
#Scatterplot accoppiati e linea di interpolazione
pairs(num_fact,
panel = panel.smooth,
main = "Matrice di Scatterplot")
Da questo primo test si evince che le varaibili continue Anni della
madre e Numero delle gravidanze sono fattori ininfluenti per il peso del
neonato. Viceversa la lunghezza, la dimensione del cranio e la durata
della gestazione sono altamente correlate con il peso in modo positivo.
Le variabili binarie Fumatrice e Tipo di parto non sembrano avere
nessuna correlazione, ma sarà bene eseguire un’analisi dedicata per
esse. Infine anche il sesso sembrerebbe avere ha un’impatto positivo sul
peso del neonato seppur in maniera decisamente inferiore rispetto alle
altre tre variabili. In questa analisi ci sono anche altre informazioni
da non trascurare, infatti sono evidenti delle correlazioni positive tra
le stesse variabili con cui costruire il modello (Gestazione, Lunghezza
e Cranio). Questo potrebbe tornare utile per ottimizzare il modello
sfruttando le interazioni tra le variabili.
Come ultima variabile è rimasta da valutare l’impatto degll’ospedale, che è una varibiale categoriale. Ciò verrà fatto sia numericamente tramite un test anova sulla sola relazione lineare Ospedale-Peso e sia visivamente tramite una rappresentazione boxplot del peso in relazione ai diversi ospedali in esame.
# Relazione numerica Peso-Ospedale
anova(lm(Peso ~ Ospedale, data = Neonati_df))## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Ospedale 2 936237 468118 1.6991 0.1831
## Residuals 2497 687952305 275512#Confronto visivo Peso-Ospedale
boxplot(Peso ~ Ospedale, data = Neonati_df, main="Peso vs Ospedale")
Tramite il test anova si evince che la variabile Ospedale ha un peso
infimo sulla spiegazione dell’andamento della varaibile peso, ciò si
nota dall’alto valore della somma dei quadrati dei residui rispetto a
quello della variabile ospedale e dal fatto che il p-value del test è
maggiore di 0.05. Anche visivamente sembra evidente il poco impatto che
ha la variabile sul peso. Tutte e tre i boxplot appaiono distribuiti nel
medesimo range e con i medesimi valori medi. In sintesi l’ospedale
presso cui sono stati eseguite le osservazioni non mostra evidenze
statisticamente significanti sulla variazione del peso alla nascita.
Proseguiamo adesso a fare un controllo analogo sulle variabili binarie che hanno dato bassi valori di correlazione nella matrice.
#Confronto visivo
par(mfrow=c(1,3))
boxplot(Peso~Fumatrici, data = Neonati_df)
boxplot(Peso~Tipo.parto_num, data = Neonati_df)
boxplot(Peso~Sesso_num, data = Neonati_df)
# T-Test
t.test(Peso~Fumatrici, data = Neonati_df)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1
## 3286.153 3236.346t.test(Peso~Tipo.parto_num, data = Neonati_df)##
## Welch Two Sample t-test
##
## data: Peso by Tipo.parto_num
## t = 0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -40.54037 46.27992
## sample estimates:
## mean in group 0 mean in group 1
## 3284.916 3282.047t.test(Peso~Sesso_num, data = Neonati_df)##
## Welch Two Sample t-test
##
## data: Peso by Sesso_num
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group 0 mean in group 1
## 3161.132 3408.215Sia dai grafici che dai test numerici si conferma ciò che si è notato inizialmente. Soltanto la variabile binaria Sesso è significativa nella variazione del peso.
2.2) Modello con Regressione Lineare Multipla Completo
Partiamo a creare un modello completo sfruttando tutte le variabili numeriche.
# Regressione Lineare Multipla Completa
lin_full_mod <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici +
Gestazione + Lunghezza + Cranio +
Tipo.parto_num + Sesso_num, data = Neonati_df)
summary(lin_full_mod)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto_num + Sesso_num, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1140.23 -181.78 -14.89 160.22 2630.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.2324 141.1162 -47.530 < 2e-16 ***
## Anni.madre 0.8647 1.1475 0.754 0.4512
## N.gravidanze 11.7148 4.6712 2.508 0.0122 *
## Fumatrici -31.5829 27.5739 -1.145 0.2522
## Gestazione 32.8391 3.8219 8.592 < 2e-16 ***
## Lunghezza 10.2723 0.3010 34.129 < 2e-16 ***
## Cranio 10.4844 0.4266 24.575 < 2e-16 ***
## Tipo.parto_num -30.2562 12.0994 -2.501 0.0125 *
## Sesso_num 78.0338 11.1925 6.972 3.99e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7279, Adjusted R-squared: 0.727
## F-statistic: 832.9 on 8 and 2491 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_full_mod))^2)), "\n")## RMSE: 273.8314Da un primo sguardo alle informazioni del modello completo si nota che come le variabili maggiormente significative sono le stese evidenziate in precedenza (Gestazione, Lunghezza, Cranio e Sesso). In questa analisi risultano marginalmente signifiative anche il numero delle gravidanze e il tipo di parto, dato che mostrano un p-value inferiore a 0.05. Infine i valori di R-squared e Adjusted R-squared sono nell’intorno del 73%. Dobbiamo cercare di ottimizzare le varibili affinchè il modello abbia un valore di Adjusted R-squared migliore.
2.2.1) Rimozione Variabili tramite Stepwise con AIC
Sfruttiamo il criterio di informazione di Akaike per ricavare le variabili più significative.
lin_aic_mod <- stepAIC(lin_full_mod, direction = "both", k=2, trace = TRUE)## Start: AIC=28080.56
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto_num + Sesso_num
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 42734 187501837 28079
## - Fumatrici 1 98728 187557831 28080
## <none> 187459103 28081
## - Tipo.parto_num 1 470578 187929681 28085
## - N.gravidanze 1 473298 187932401 28085
## - Sesso_num 1 3658038 191117141 28127
## - Gestazione 1 5555847 193014950 28152
## - Cranio 1 45450123 232909226 28621
## - Lunghezza 1 87653102 275112206 29038
##
## Step: AIC=28079.13
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto_num + Sesso_num
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 99840 187601677 28079
## <none> 187501837 28079
## + Anni.madre 1 42734 187459103 28081
## - Tipo.parto_num 1 471817 187973654 28083
## - N.gravidanze 1 675718 188177555 28086
## - Sesso_num 1 3665908 191167745 28126
## - Gestazione 1 5514506 193016343 28150
## - Cranio 1 45711714 233213551 28623
## - Lunghezza 1 87642114 275143951 29036
##
## Step: AIC=28078.46
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto_num +
## Sesso_num
##
## Df Sum of Sq RSS AIC
## <none> 187601677 28079
## + Fumatrici 1 99840 187501837 28079
## + Anni.madre 1 43845 187557831 28080
## - Tipo.parto_num 1 463870 188065546 28083
## - N.gravidanze 1 651066 188252743 28085
## - Sesso_num 1 3649259 191250936 28125
## - Gestazione 1 5444109 193045786 28148
## - Cranio 1 45758101 233359778 28622
## - Lunghezza 1 88054432 275656108 29039summary(lin_aic_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto_num + Sesso_num, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.31 -181.70 -16.31 161.07 2638.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6677.2629 135.5916 -49.245 < 2e-16 ***
## N.gravidanze 12.7558 4.3366 2.941 0.0033 **
## Gestazione 32.2713 3.7941 8.506 < 2e-16 ***
## Lunghezza 10.2864 0.3007 34.207 < 2e-16 ***
## Cranio 10.5057 0.4260 24.659 < 2e-16 ***
## Tipo.parto_num -30.0342 12.0969 -2.483 0.0131 *
## Sesso_num 77.9285 11.1905 6.964 4.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1110 on 6 and 2493 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_aic_mod))^2)), "\n")## RMSE: 273.9355Il metodo stepwise AIC, come è noto, tendere a essere meno severo nel rimuovere delle variabili con l’ottica di essere più preciso ed infatti in questo caso non ha escluso tutte le varaibili che durante il controllo con la matrice di correlazione erano poco influenti al peso (N.gravidanze e Tipo.parto). Tuttavia R-Squared non ha sortito grossi miglioramenti rispetto al modello lineare completo.
2.2.2) Rimozione Variabili tramite Stepwise con BIC
n <- nrow(Neonati_df)
lin_bic_mod <- stepAIC(lin_full_mod, direction = "both", k=log(n), trace = TRUE)## Start: AIC=28132.98
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto_num + Sesso_num
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 42734 187501837 28126
## - Fumatrici 1 98728 187557831 28127
## - Tipo.parto_num 1 470578 187929681 28131
## - N.gravidanze 1 473298 187932401 28132
## <none> 187459103 28133
## - Sesso_num 1 3658038 191117141 28174
## - Gestazione 1 5555847 193014950 28198
## - Cranio 1 45450123 232909226 28668
## - Lunghezza 1 87653102 275112206 29084
##
## Step: AIC=28125.73
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto_num + Sesso_num
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 99840 187601677 28119
## - Tipo.parto_num 1 471817 187973654 28124
## <none> 187501837 28126
## - N.gravidanze 1 675718 188177555 28127
## + Anni.madre 1 42734 187459103 28133
## - Sesso_num 1 3665908 191167745 28166
## - Gestazione 1 5514506 193016343 28190
## - Cranio 1 45711714 233213551 28663
## - Lunghezza 1 87642114 275143951 29077
##
## Step: AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto_num +
## Sesso_num
##
## Df Sum of Sq RSS AIC
## - Tipo.parto_num 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 43845 187557831 28127
## - Sesso_num 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso_num
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto_num 1 463870 187601677 28119
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 45044 188020502 28125
## - Sesso_num 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(lin_bic_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso_num, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## Sesso_num 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_bic_mod))^2)), "\n")## RMSE: 274.274In questo caso, rispetto al metodo AIC, il metodo di rimozione BIC è stato leggermente più aggressivo nella rimozione delle varaibili. Ha mantenuto le variabili che avevamo riscontrato con la matrice di correlazione essere le più significative, ma anche in questo caso non ha escluso la varaibile N.gravidanze. Seppur si sia abbassato leggermente il valore di R-Squared rispetto al modello completo e quello post stepwise AIC, la variazione è trascuraible rispetto alla semplificazione del modello così ottenuta.
2.2.3) Modelli con Effetti di ordine superiore
Partendo dal modello ridotto ottenuto tramite il criterio di Bayes si valuteranno eventuali cambi di perfomance introducendo le interazioni tra le varibili di secondo ordine e possibili effetti quadratici.
lin_ord2_mod <- lm(Peso ~ (N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso_num)^2 +
I(N.gravidanze^2) + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + I(Sesso_num^2),
data = Neonati_df)
summary(lin_ord2_mod)##
## Call:
## lm(formula = Peso ~ (N.gravidanze + Gestazione + Lunghezza +
## Cranio + Sesso_num)^2 + I(N.gravidanze^2) + I(Gestazione^2) +
## I(Lunghezza^2) + I(Cranio^2) + I(Sesso_num^2), data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1193.3 -179.2 -10.0 166.2 1329.6
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.866e+02 1.254e+03 0.707 0.47962
## N.gravidanze -1.324e+02 8.494e+01 -1.559 0.11920
## Gestazione 1.506e+02 7.665e+01 1.965 0.04957 *
## Lunghezza -3.215e+00 6.021e+00 -0.534 0.59335
## Cranio -2.896e+01 1.057e+01 -2.740 0.00619 **
## Sesso_num -6.479e+01 2.750e+02 -0.236 0.81378
## I(N.gravidanze^2) -3.603e+00 1.334e+00 -2.700 0.00698 **
## I(Gestazione^2) -5.160e+00 1.633e+00 -3.160 0.00160 **
## I(Lunghezza^2) 5.755e-02 6.715e-03 8.571 < 2e-16 ***
## I(Cranio^2) 4.576e-02 1.912e-02 2.393 0.01679 *
## I(Sesso_num^2) NA NA NA NA
## N.gravidanze:Gestazione -3.716e+00 2.731e+00 -1.361 0.17365
## N.gravidanze:Lunghezza 5.472e-02 2.374e-01 0.230 0.81775
## N.gravidanze:Cranio 8.166e-01 3.566e-01 2.290 0.02210 *
## N.gravidanze:Sesso_num 1.026e+01 8.947e+00 1.147 0.25167
## Gestazione:Lunghezza -3.194e-01 1.992e-01 -1.603 0.10897
## Gestazione:Cranio 1.304e+00 2.261e-01 5.769 8.94e-09 ***
## Gestazione:Sesso_num 1.106e+01 7.688e+00 1.439 0.15036
## Lunghezza:Cranio -8.782e-02 1.699e-02 -5.167 2.57e-07 ***
## Lunghezza:Sesso_num -5.888e-01 6.147e-01 -0.958 0.33826
## Cranio:Sesso_num -3.450e-02 8.632e-01 -0.040 0.96812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.8 on 2480 degrees of freedom
## Multiple R-squared: 0.7457, Adjusted R-squared: 0.7437
## F-statistic: 382.7 on 19 and 2480 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_ord2_mod))^2)), "\n")## RMSE: 264.7266Il modello ottenuto è decisamente più complesso rispetto a quello iniziale. I fattori che sembrano avere più significatività sono quelli di ordine superiore, in particolare Lunghezza^2, Gestazione-Cranio e Lunghezza-Cranio. Il sesso del neonato sembra non avere più rilevanza in questo modello, mentre il numero di gravidanze sembra dare un effetto marginale. L’Adjusted R-squared di questo modello è aumentato (74%) come è anche diminuito RMSE rispetto ai modelli precedenti, un miglioramento non così accentuato in relazione alla complessità del modello ottenuto.
Riproviamo ad eseguire le tecniche stepwise AIC e BIC su questo nuovo modello.
lin_aic_ord2_mod <- stepAIC(lin_ord2_mod, direction = "both", k=2, trace = FALSE)
summary(lin_aic_ord2_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso_num + I(N.gravidanze^2) + I(Gestazione^2) + I(Lunghezza^2) +
## I(Cranio^2) + N.gravidanze:Cranio + Gestazione:Lunghezza +
## Gestazione:Cranio + Lunghezza:Cranio, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1194.65 -180.79 -13.32 166.51 1332.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.043e+03 1.237e+03 0.843 0.39910
## N.gravidanze -1.879e+02 7.678e+01 -2.447 0.01448 *
## Gestazione 1.323e+02 7.536e+01 1.755 0.07937 .
## Lunghezza -2.308e+00 5.938e+00 -0.389 0.69754
## Cranio -2.913e+01 1.033e+01 -2.820 0.00484 **
## Sesso_num 7.302e+01 1.090e+01 6.700 2.56e-11 ***
## I(N.gravidanze^2) -3.350e+00 1.269e+00 -2.639 0.00836 **
## I(Gestazione^2) -4.770e+00 1.614e+00 -2.955 0.00316 **
## I(Lunghezza^2) 5.550e-02 6.502e-03 8.537 < 2e-16 ***
## I(Cranio^2) 5.146e-02 1.869e-02 2.754 0.00593 **
## N.gravidanze:Cranio 6.476e-01 2.255e-01 2.872 0.00412 **
## Gestazione:Lunghezza -2.873e-01 1.965e-01 -1.462 0.14378
## Gestazione:Cranio 1.227e+00 2.129e-01 5.762 9.34e-09 ***
## Lunghezza:Cranio -8.887e-02 1.671e-02 -5.319 1.14e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.7 on 2486 degrees of freedom
## Multiple R-squared: 0.7452, Adjusted R-squared: 0.7438
## F-statistic: 559.2 on 13 and 2486 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_aic_ord2_mod))^2)), "\n")## RMSE: 264.9913n <- nrow(Neonati_df)
lin_bic_ord2_mod <- stepAIC(lin_ord2_mod, direction = "both", k=log(n), trace = FALSE)
summary(lin_bic_ord2_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso_num + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) +
## Gestazione:Cranio + Lunghezza:Cranio, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.37 -182.19 -11.71 168.33 1318.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.032e+02 1.190e+03 -0.171 0.864479
## N.gravidanze 1.483e+01 4.216e+00 3.516 0.000445 ***
## Gestazione 1.802e+02 7.431e+01 2.425 0.015395 *
## Lunghezza -7.234e+00 5.657e+00 -1.279 0.201045
## Cranio -2.058e+01 1.006e+01 -2.046 0.040877 *
## Sesso_num 7.275e+01 1.092e+01 6.662 3.30e-11 ***
## I(Gestazione^2) -6.298e+00 1.032e+00 -6.102 1.21e-09 ***
## I(Lunghezza^2) 5.048e-02 5.042e-03 10.012 < 2e-16 ***
## I(Cranio^2) 5.483e-02 1.839e-02 2.981 0.002901 **
## Gestazione:Cranio 1.014e+00 2.056e-01 4.932 8.69e-07 ***
## Lunghezza:Cranio -9.269e-02 1.579e-02 -5.870 4.94e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.5 on 2489 degrees of freedom
## Multiple R-squared: 0.7435, Adjusted R-squared: 0.7424
## F-statistic: 721.3 on 10 and 2489 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_bic_ord2_mod))^2)), "\n")## RMSE: 265.8746Anche in questo caso il modello ottenuto tramite stepwise BIC è risultato essere più snello rispetto a quello AIC e i valori delle metriche sono migliori rispetto alle controparti prive di elementi di ordine superiore.
Come ultimo step si può valutare di rimuovere manualmente la varaibile Numero di gravidanze dato che in una prima fase di selezione sembrava essere poco o nulla influente sul fattore peso ed è anche l’unica variabile che nell’ulitmo modello ottenuto non sembra avere effetti di ordine superiore.
lin_bic_ord2_red_mod <- update(lin_bic_ord2_mod,~.-N.gravidanze)
summary(lin_bic_ord2_red_mod)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso_num +
## I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + Gestazione:Cranio +
## Lunghezza:Cranio, data = Neonati_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1153.54 -183.48 -11.87 168.55 1306.70
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -225.33592 1192.90777 -0.189 0.85019
## Gestazione 177.71542 74.47582 2.386 0.01710 *
## Lunghezza -6.80501 5.66829 -1.201 0.23004
## Cranio -20.62418 10.08123 -2.046 0.04088 *
## Sesso_num 74.15643 10.93780 6.780 1.50e-11 ***
## I(Gestazione^2) -6.20920 1.03409 -6.004 2.20e-09 ***
## I(Lunghezza^2) 0.04981 0.00505 9.863 < 2e-16 ***
## I(Cranio^2) 0.05571 0.01843 3.022 0.00253 **
## Gestazione:Cranio 0.99719 0.20602 4.840 1.38e-06 ***
## Lunghezza:Cranio -0.09218 0.01583 -5.824 6.48e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.1 on 2490 degrees of freedom
## Multiple R-squared: 0.7422, Adjusted R-squared: 0.7413
## F-statistic: 796.5 on 9 and 2490 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df$Peso - predict(lin_bic_ord2_red_mod))^2)), "\n")## RMSE: 266.5342La rimozione della variabile numero di gravidanze ha comportato un leggero calo delle perfomance, tuttavia è del tutto accettabile dato che è pur sempre bene scegliere un modello parsimonioso di variabili. Per tale motivo verrà scelto come modello finale lin_bic_ord2_red_mod
#2.3) Valutazione del modello Adesso che è stato scelto il modello e valutato le metriche R^2 e RMSE si vuole porre l’attenzione sull’analisi dei residui in modo da identificare eventuali anomalie.
# Plot dei grafici di analisi del modello
par(mfrow = c(2, 2))
plot(lin_bic_ord2_red_mod)
Dai risultati ottenuti si nota che: - Dal primo grafico i residui
sembrano distribuiti in modo casuale attorno allo zero con una leggera
dispersione accentuanta per valori di output più alti, che potrebbe
indicare eteroscedasticità leggera. - Dal grafico Q-Q i punti centrali
si allineano bene con la retta tratteggiata ciò implica una buona
approssimazione alla normalità, tuttavia le code (soprattutto quella
superiore) mostrano deviazioni dalla normalità riconducibili ad alcuni
outlier. - Dal grafico Scale-Location i punti sono relativamente
dispersi in modo uniforme, ma c’è un leggero aumento della varianza
all’aumentare dei valori predetti sintomo di una potenziale lieve
eteroscedasticità. - Infine dal grafico del Leverage si nota che la
maggior parte dei punti ha basso leverage, tuttavia alcuni punti (1551,
1429, 310) sono più distanti dal gruppo e tendono ad essere vicini se
non addirittua oltre le linee di Cook’s distance.
Questi punti anomali potrebbero essere influenti e hanno bisogno di un’analisi approfondita.
#2.3.1) Analisi dei punti influenti
#Leverage
lev <- hatvalues(lin_bic_ord2_red_mod)
#Soglia del valore di Leverage
p = sum(lev)
threshold_lev = 2*p/n
# Grafico dei punti
plot(lev)
abline(h=threshold_lev, col="red")
Dal grafico dei punti in base al Leverage si nota che la soglia è
prossima allo zero, questo comporta un elevato di punti al di sopra di
essa. Tuttavia si nota come soltanto uno in particolare sia al di sopra
dello 0.5. per tale motivo é bene identificarlo.
# Identificazione del punto estremo
lev[lev>0.5]## 1551
## 0.6534458Il punto identificato è il 1551, lo stesso che era risulato evidente durante l’analisi grafica dei residui.
Passiamo adesso ad una valutazione degli outliers.
#Outliers
resid <- rstudent(lin_bic_ord2_red_mod)
#Grafico dei punti
plot(resid)
abline(h=c(-2,2), col=2)
# Test di identificazione degli outliers
outlierTest(lin_bic_ord2_red_mod)## rstudent unadjusted p-value Bonferroni p
## 1306 4.920817 9.1813e-07 0.0022953
## 155 4.520239 6.4660e-06 0.0161650
## 1694 4.373643 1.2722e-05 0.0318040
## 1399 -4.340729 1.4769e-05 0.0369210Dal grafico anche in questo risultano essere molti i punti al di fuori o in prossimità delle sogli e di ouliers, ma dal test eseguito è emerso che soltanto 4 punti sono da considerare effettivamente come punti anomali.
Come ultimo test sui punti passiamo al calcolo della distanza di Cook.
# Calcolo della Cook's distance
cooksD <- cooks.distance(lin_bic_ord2_red_mod)
#Grafico dei punti
plot(cooksD)
abline(h=c(0.5,1), col="red")
é evidente dal grafico che tutti i punti ad esclusione di uno siano
lontatni dalle soglie di criticità di Cook. Passiamo adesso ad
indentificare il punto che sta ben oltre la soglia di 1.
# Identificazione del punto estremo
cooksD[cooksD>1]## 1551
## 2.92047# Soglia di Cook's Distance
influential_point <- which(cooksD > 1)
# Dataset con sole osservazioni influenti
Neonati_influenti_df <- Neonati_df[influential_point, ]
Neonati_influenti_df## # A tibble: 1 × 12
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 35 1 0 38 4370 315 374 Nat
## # ℹ 4 more variables: Ospedale <chr>, Sesso <chr>, Tipo.parto_num <dbl>,
## # Sesso_num <dbl>Anche in questo caso il punto più influente è lo stesso visto durante l’analisi tramite Leverage, cioè il punto 1551. Quello che si nota da questa osservazione è che il peso è insolitamente alto rispetto alla lunghezza del neonato se confrontato coi valori medi e mediani riportati nella tabella dell’analisi statistica. Qeusto ci spiega perchè il modello non riesce a fittarlo in maniera idonea. Si potrebbe ipotizzare di riomuovere tale osservazione e ricalcolare le metriche del modello.
# Rimozione dei punti influenti
Neonati_df_clean <- Neonati_df[-influential_point, ]
# Riadattamento del modello
lin_mod_clean <- lm(Peso ~ Gestazione + Lunghezza + Cranio +
Sesso_num + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) +
Gestazione:Cranio + Lunghezza:Cranio, data = Neonati_df_clean)
summary(lin_mod_clean)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso_num +
## I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + Gestazione:Cranio +
## Lunghezza:Cranio, data = Neonati_df_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1158.05 -182.66 -12.72 163.24 1303.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.42565 1192.23828 0.082 0.934879
## Gestazione 165.03521 74.32820 2.220 0.026484 *
## Lunghezza -5.29453 5.66471 -0.935 0.350058
## Cranio -23.20248 10.07304 -2.303 0.021337 *
## Sesso_num 73.76967 10.90636 6.764 1.67e-11 ***
## I(Gestazione^2) -4.80866 1.09043 -4.410 1.08e-05 ***
## I(Lunghezza^2) 0.02655 0.00775 3.426 0.000623 ***
## I(Cranio^2) 0.02961 0.01953 1.516 0.129619
## Gestazione:Cranio 0.71174 0.21778 3.268 0.001097 **
## Lunghezza:Cranio -0.02875 0.02252 -1.276 0.201952
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.3 on 2489 degrees of freedom
## Multiple R-squared: 0.7434, Adjusted R-squared: 0.7424
## F-statistic: 801 on 9 and 2489 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df_clean$Peso - predict(lin_mod_clean))^2)), "\n")## RMSE: 265.7571La sola rimozione di quel punto ha migliorato le metriche del modello ed inoltre sembra aver ridotto l’influenza di alcuni effetti (Cranio^2 e Lunghezza-Cranio).
A questo punto si può provare ad applicare un ultima correzione a tali effetti rimuovendoli.
# Riadattamento del modello
lin_red_mod_clean <- lm(Peso ~ Gestazione + Lunghezza + Cranio +
Sesso_num + I(Gestazione^2) + I(Lunghezza^2) +
Gestazione:Cranio, data = Neonati_df_clean)
summary(lin_red_mod_clean)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso_num +
## I(Gestazione^2) + I(Lunghezza^2) + Gestazione:Cranio, data = Neonati_df_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1163.23 -182.39 -13.11 165.10 1317.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.786e+02 1.086e+03 -0.349 0.7273
## Gestazione 1.586e+02 6.668e+01 2.379 0.0174 *
## Lunghezza -9.287e+00 5.089e+00 -1.825 0.0681 .
## Cranio -1.385e+01 5.949e+00 -2.327 0.0200 *
## Sesso_num 7.356e+01 1.091e+01 6.745 1.89e-11 ***
## I(Gestazione^2) -4.346e+00 1.031e+00 -4.216 2.58e-05 ***
## I(Lunghezza^2) 2.070e-02 5.183e-03 3.993 6.71e-05 ***
## Gestazione:Cranio 6.242e-01 1.538e-01 4.060 5.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7431, Adjusted R-squared: 0.7424
## F-statistic: 1029 on 7 and 2491 DF, p-value: < 2.2e-16cat("RMSE:", sqrt(mean((Neonati_df_clean$Peso - predict(lin_red_mod_clean))^2)), "\n")## RMSE: 265.8958# Plot dei grafici di analisi del modello
par(mfrow = c(2, 2))
plot(lin_red_mod_clean)
A questo punto il modello si può considerare ultimato. Le variabili prese in esame sono: Gestazione, Lunghezza, Cranio e Sesso_num. Inoltre sono significativi anche i seguenti effetti di secondo ordine: Gestazione^2, Lunghezza^2 e Gestazione-Cranio.
#3) Esempi di predizione del modello
Facciamo adesso delle previsioni pratiche fornendo un set di dati nuovi e vediamo cosa predice il modello.
# Creazione di un dataset con più osservazioni
new_data_batch <- data.frame(
Anni.madre = c(13, 25, 28, 32, 46),
N.gravidanze = c(0, 0, 1, 2, 3),
Fumatrici = c(0, 1, 0, 0, 1),
Gestazione = c(25, 38, 39, 40, 43),
Lunghezza = c(310, 480, 500, 510, 550),
Cranio = c(235, 330, 340, 350, 380),
Tipo.parto = c("Nat", "Ces", "Ces", "Nat", "Nat"),
Ospedale = c("osp1", "osp3", "osp2", "osp3", "osp2"),
Sesso = c("F", "M", "F", "M", "M"),
Tipo.parto_num = c(0, 1, 1, 0, 0),
Sesso_num = c(0, 1, 0, 1, 1)
)
# Previsioni
predicted_weights <- predict(lin_bic_ord2_mod, newdata = new_data_batch)
# Visualizzazione dei risultati
print(round(predicted_weights, 2))## 1 2 3 4 5
## 370.00 2993.01 3291.43 3634.15 4601.28Dagli esempi è evidente che più tendiamo a valori estremi, specialemnte limiti inferiori, più il modello tende ad avere predizioni improbabili (vedi punto il 1 che riporta un peso di 370 grammi).
#4) Visualizzazione delle relazioni del modello
Per concludere questo progetto è bene dare delle informazioni rapide e visive su come il modello opera. Per fare ciò eseguiremo delle rappresentazioni che mettono in relazioni le varaibili selezionate con la varaibile target peso.
par(mfrow = c(1, 3))
# Peso vs Gestazione
ggplot(data = Neonati_df_clean) +
geom_point(aes(x=Gestazione, y=Peso, col=Sesso)) +
geom_smooth(aes(x=Gestazione, y=Peso, col=Sesso),
method = "lm", formula = y ~ poly(x, 2)) +
labs(title = "Relazione tra Gestazione e Peso",
x = "Settimane di gestazione", y = "Peso (g)")
# Peso vs Lunghezza Neonato
ggplot(data = Neonati_df_clean) +
geom_point(aes(x=Lunghezza, y=Peso, col=Sesso)) +
geom_smooth(aes(x=Lunghezza, y=Peso, col=Sesso),
method = "lm", formula = y ~ poly(x, 2)) +
labs(title = "Relazione tra Lunghezza e Peso",
x = "Lunghezza misurata (mm)", y = "Peso (g)")
# Peso vs Dimensione Cranio
ggplot(data = Neonati_df_clean) +
geom_point(aes(x=Cranio, y=Peso, col=Sesso)) +
geom_smooth(aes(x=Cranio, y=Peso, col=Sesso),
method = "lm", formula = y ~ poly(x, 2)) +
labs(title = "Relazione tra Cranio e Peso",
x = "Dimensione Cranio misurata (mm)", y = "Peso (g)")
Da questi primi tre grafici sono più evidenti le relazioni non lineari delle varaibili e di come esse siano distinte in base al sesso del neonato. Inoltre si nota quello che avevamo ipotizzato durante gli esempi di predizione, cioè che il modello tende ad essere meno accurato in prossimità di valori estremi inferiori. questo è dovuto anche al fatto che la maggior parte delle nostre osservazioni con cui è stato costruito il modello era concentrato su dei valori medi.
Come ulteriore informazione da rappresentare, oltre a quelle già esposte è l’effetto composto tra Gestazione e Dimensione del Cranio come riporato all’interno del modello. Per fare si può sfruttare un grafico a contorno 2D distinto per entrambi i sessi in cui il valore della Lunghezza è fissata al valore medio.
# Griglia di valori
gest_range <- seq(25, 50, length.out = 50)
cran_range <- seq(200, 500, length.out = 50)
grid <- expand.grid(Gestazione = gest_range,
Cranio = cran_range,
Sesso_num = c(0, 1)) # 0 = femmina, 1 = maschio
# Aggiungi le altre variabili fisse
grid$Lunghezza <- median(Neonati_df_clean$Lunghezza)
grid$`I(Gestazione^2)` <- grid$Gestazione^2
grid$`I(Lunghezza^2)` <- grid$Lunghezza^2
grid$`Gestazione:Cranio` <- grid$Gestazione * grid$Cranio
# Calcola le predizioni
grid$Peso <- predict(lin_red_mod_clean, newdata = grid)
# Aggiungi etichetta per sesso
grid$Sesso <- ifelse(grid$Sesso_num == 1, "Maschio", "Femmina")
# Grafico a contorno con faceting
ggplot(grid, aes(x = Gestazione, y = Cranio, z = Peso)) +
geom_contour_filled() +
facet_wrap(~ Sesso) +
labs(title = "Peso previsto in funzione di Gestazione e Cranio per Sesso",
x = "Gestazione (sett.)", y = "Cranio (mm)", fill = "Peso (g)") +
theme_minimal()