Modello Statistico per la Previsione del Peso Neonatale (v2)
Luca Marletta
06/01/2025
A. TESTO DELLA CONSEGNA
(PROFESSION AI)
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) e il peso del neonato.
Creazione del Modello di Regressione
Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato. Ad esempio, ci aspettiamo che il fumo materno e le gravidanze multiple abbiano effetti negativi sul peso alla nascita, mentre una maggiore durata della gestazione potrebbe aumentare il peso del neonato.
Selezione del Modello Ottimale
Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.
Analisi della Qualità del Modello
Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Mean Squared Error (MSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni.
Previsioni e Risultati
Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice, che partorirà alla 39esima settimana.
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 fumo sul peso previsto.
—————————————
B. SVOLGIMENTO DEL PROGETTO
(LUCA MARLETTA)
1. ANALISI PRELIMINARE
1.1 Import del dataset e delle librerie necessarie
#library(MASS)
library(GGally)
library(patchwork)
library(ggpubr)
library(kableExtra)
library(knitr)
library(RColorBrewer)
library(colorspace)
library(summarytools)
library(dplyr)
library(ggplot2)
library(moments)
library(gridExtra)
library(cowplot)
library(tidyr)
library(TeachingDemos)
library(lmtest)
library(car)
library(purrr)
#import dataset
dati <- read.csv("neonati.csv")
attach(dati)Compiliamo una funziona rapida per formattare rapidamente le tabelle dati in formato “html”, tramite la libreria kable():
kable_fun <- function(tabella, description){
kable(tabella,
"html",
caption = description) %>%
kable_styling(full_width = F,
position ="center",
bootstrap_options = c(
"striped",
"hover",
"condensed",
"responsive"))}1.2 - Analisi descrittiva univariata del dataset
Con il blocco di codice sottostante, verifico l’eventuale presenza di valori nulli all’interno del dataset:
Na_table <- as.data.frame(colSums(is.na(dati))) colnames(Na_table) <- c("Conteggio valori nulli") kable_fun(Na_table, "Tab.0 - Verifica di valori nulli all'interno della dataset")Tab.0 - Verifica di valori nulli all’interno della dataset
Conteggio valori nulli
Anni.madre
0
N.gravidanze
0
Fumatrici
0
Gestazione
0
Peso
0
Lunghezza
0
Cranio
0
Tipo.parto
0
Ospedale
0
Sesso
0
dall’analisi iniziale si evince che non sono presenti valori nulli all’interno del dataset, quindi procediamo con l’analisi della struttura:
n<-nrow(dati) # Cattura l'output della funzione str() str_output <- capture.output(str(dati)) # Converte l'output in un data frame str_df <- as.data.frame(matrix(str_output, nrow = length(str_output), byrow = TRUE)) # Imposta i nomi delle colonne (opzionale) colnames(str_df) <- c(" ") kable_fun(str_df,"Tab1. - Panoramica generale della struttura dataset")Tab1. - Panoramica generale della struttura dataset
‘data.frame’: 2500 obs. of 10 variables:
$ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 …
$ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 …
$ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 …
$ Gestazione : int 42 39 38 41 38 40 39 40 40 41 …
$ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 …
$ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 …
$ Cranio : int 325 345 375 365 335 340 345 349 335 362 …
$ Tipo.parto : chr “Nat” “Nat” “Nat” “Nat” …
$ Ospedale : chr “osp3” “osp1” “osp2” “osp2” …
$ Sesso : chr “M” “F” “M” “M” …
1.2.1 - Descrizione delle variabili
Il dataset è composto da 10 variabili, con 2500 per ogni variabile:
Anni Madre: variabile quantitativa continua su scala di rapporti, ed intuitivamente è la variabile che indica il numero di anni della madre espressa in numeri interi;
Fumatrici: variabile categoriali, espressa con numeri interi e codificata in codice binario, Fumatrice=0 (baseline) indica madre NON fumatrice, ,mentre la codifica Fumatrice=1 indica una madre fumatrice;
N.di gravidanze: variabile quantitativa su scala di rapporti, espressa con numeri interi. Indica il numero di gravidanze pregresse sostenute dalla madre al momento del parto;
Gestazione: variabile quantitativa su scala di rapporti, espressa in numeri interi, ed indica il numero delle settimane di gestazione totali della madre al momento del parto.
Peso del neonato: variabile quantitativa continua su scala di rapporti, espressa in grammi ed indica il peso del neonato al momento del parto
Lunghezza : variabile quantitativa continua su scala di rapporti, espressa in millimetri, che indica la lunghezza del neonato rinvenuta tramite misure ecografiche;
Cranio: variabile quantitativa continua su scala di rapporti, espressa in millimetri, che indica il diametro craniale del neonato rinvenuta tramite misure ecografiche;
Tipo di parto: variabile qualitativa categorica, che può assumere due valori codificati in lettere: indicano il tipo di parto sostenuto dalla madre il giorno del parto, ovvero “Nat” (parto naturale) o “Ces” (parto cesareo).
Ospedale di nascita: variabile qualitativa categorica, che può assumere tre valori, codificata in lettere, “Osp1”,“Osp2”,“Osp3”, che identifica sequenzialmente in quale struttura ospedaliera le varie madre registrate nel dataset hanno partorito;
Sesso del neonato: variabile qualitativa categorica, che può assumere due valori codificati in lettere , “M” (neonato maschio), “F” (neonato femmina).
Per una prima analisi esplorativa del dataset, tramite la libreria dplyr, calcoliamo per ogni variabile numerica quantitativa(Anni.madre, N. gravidanze, Gestazione, Peso, Lunghezza, Cranio) gli indici di posizione e di variabilità, aggregando i risultati nella tabella1. A seguire verranno calcolate le distribuzioni di frequenza e i relativi grafici. Per le variabili categoriche restanti (Fumatrici, Tipo.parto, Ospedale, Sesso) veranno costruite le distribuzioni di frequenza e calcolati gli indici di eterogeneità di Gini dei rispettivi vettori numerici.
1.2.2 - La prima analisi esplorativa sul Dataset
quant.variables <- c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio", "Peso")
tab1 <- dati %>%
select(all_of(quant.variables)) %>%
reframe(
min = round(sapply(., min), 2),
Q1= round(sapply(., function(x) quantile(x, 0.25)), 2),
median= round(sapply(., median), 2),
mean= round(sapply(., mean), 2),
Q3= round(sapply(., function(x) quantile(x, 0.75)), 2),
max = round(sapply(., max), 2),
STD_dev= round(sapply(., sd), 2),
skewness = round(sapply(., skewness), 2),
kurtosis= round(sapply(., function(x) kurtosis(x) - 3), 2))%>%
mutate(CV=round((STD_dev/mean),2),
Range=max-min,
IQR=Q3-Q1)%>%
relocate(CV, .before=skewness)
tab1<-t(tab1)
colnames(tab1) <- c("Anni.madre",
"N.Gravidanze",
"Settimane Gestazione",
"Lunghezza(mm)",
"Diametro craniale(mm)",
"Peso (g)")
kable_fun(tab1, "Tab.2 Indici di posizione/variabilità calcolati sulle variabili quantitative del dataset")| Anni.madre | N.Gravidanze | Settimane Gestazione | Lunghezza(mm) | Diametro craniale(mm) | Peso (g) | |
|---|---|---|---|---|---|---|
| min | 0.00 | 0.00 | 25.00 | 310.00 | 235.00 | 830.00 |
| Q1 | 25.00 | 0.00 | 38.00 | 480.00 | 330.00 | 2990.00 |
| median | 28.00 | 1.00 | 39.00 | 500.00 | 340.00 | 3300.00 |
| mean | 28.16 | 0.98 | 38.98 | 494.69 | 340.03 | 3284.08 |
| Q3 | 32.00 | 1.00 | 40.00 | 510.00 | 350.00 | 3620.00 |
| max | 46.00 | 12.00 | 43.00 | 565.00 | 390.00 | 4930.00 |
| STD_dev | 5.27 | 1.28 | 1.87 | 26.32 | 16.43 | 525.04 |
| CV | 0.19 | 1.31 | 0.05 | 0.05 | 0.05 | 0.16 |
| skewness | 0.04 | 2.51 | -2.07 | -1.51 | -0.79 | -0.65 |
| kurtosis | 0.38 | 10.99 | 8.26 | 6.49 | 2.95 | 2.03 |
| Range | 46.00 | 12.00 | 18.00 | 255.00 | 155.00 | 4100.00 |
| IQR | 7.00 | 1.00 | 2.00 | 30.00 | 20.00 | 630.00 |
Dalla analisi preliminare di tabella1 risalta subito all’occhio che Anni madre ha un range di 46 anni con valori minimi di 0 anni (dato impossibile) a 46 anni (età massima registrata). Decido di indagare con boxplot la distribuzione del vettore:
boxplot(Anni.madre,
main="Boxplot di Anni.madre",
xlab="Anni madre",
ylab="Distribuzione",
col="skyblue",
border="black")
In effetti si osservano 2 punti inferiori della distribuzione (prossimi allo 0 della scala in anni), che mi fa presumere con molta probabilità che all’interno del dataset vi sia stato un’errore di compilazione delle osservazioni, data che è impossibile che madri di 0 anni possano partorire un neonato. Nel passaggio successivo compilo una funzione che identifica le osservazioni outliers del vettore Anni. Madre
1.2.3 - Analisi outlier sul vettore Anni.Madre
- Compilazione di una funzione identity_outliers() che permette di identificare gli outliers (in modo univariato) inferiori (lowerboud) e superiori (upperbound) per ogni distribuzione numerica, al fine di quantificare quanto visto graficamente tramite i boxplots appena illustrato:
outliers_analysis <- as.data.frame(cbind(Anni.madre, Gestazione,Peso, Lunghezza,Cranio, N.gravidanze))
#definizione funzione outliers lowerbound
identify_outliers_lower <- function(df, col) {
q1 <- quantile(df[[col]], 0.25)
q3 <- quantile(df[[col]], 0.75)
iqr <- IQR(df[[col]])
outs_lower <- df[df[, col] < q1 - 1.5 * iqr,]
return(outs_lower)}
#definizione funzione outliers upperbound
identify_outliers_upper <- function(df, col) {
q1 <- quantile(df[[col]], 0.25)
q3 <- quantile(df[[col]], 0.75)
iqr <- IQR(df[[col]])
outs_upper <- df[df[, col] > q3 + 1.5 * iqr,]
return(outs_upper)}
kable_fun(identify_outliers_lower(outliers_analysis,"Anni.madre"), "Tab.3 - Outliers inferiori vettore Anni.Madre")| Anni.madre | Gestazione | Peso | Lunghezza | Cranio | N.gravidanze | |
|---|---|---|---|---|---|---|
| 138 | 13 | 38 | 2760 | 470 | 325 | 0 |
| 1075 | 14 | 39 | 3510 | 490 | 365 | 1 |
| 1152 | 1 | 41 | 3250 | 490 | 350 | 1 |
| 1380 | 0 | 39 | 3060 | 490 | 330 | 0 |
| 1532 | 14 | 39 | 3550 | 500 | 355 | 0 |
Attenzionando l’analisi degli outliers inferiori al lowerbound (valori inferiori a Q1*1.5IQR) del relativo vettore dati, appare evidente che ci siano almeno 2 osservazioni che riportano delle età delle madri pari a 1 e 0 (rispettivamente la 1152 e 1380), un età che non è possibile per ovvi motivi.
Decido quindi di costruire un nuovo dataset a partire da quello originale, eliminando queste due righe dal dataset originale, e porre questo nuovo set di dati come base per tutte le osservazioni che seguiranno.
1.2.4 VARIABILI QUANTITATIVE - Indici di posizione, variabilità e distribuzione di frequenze per le variabili numeriche quantitative (Anni. madre, N.gravidanze, Gestazione, Lunghezza, Cranio e Peso)
dati <- subset(dati, !(Anni.madre %in% c(0, 1)))
attach(dati)
tab2 <- dati %>%
select(all_of(quant.variables)) %>%
reframe(
min = round(sapply(., min), 2),
Q1= round(sapply(., function(x) quantile(x, 0.25)), 2),
median= round(sapply(., median), 2),
mean= round(sapply(., mean), 2),
Q3= round(sapply(., function(x) quantile(x, 0.75)), 2),
max = round(sapply(., max), 2),
STD_dev= round(sapply(., sd), 2),
skewness = round(sapply(., skewness), 2),
kurtosis= round(sapply(., function(x) kurtosis(x) - 3), 2))%>%
mutate(CV=round((STD_dev/mean),2),
Range=max-min,
IQR=Q3-Q1)%>%
relocate(CV, .before=skewness)
tab2<-t(tab2)
colnames(tab2) <- c("Anni.madre",
"N.Gravidanze",
"Settimane Gestazione",
"Lunghezza(mm)",
"Diametro craniale(mm)",
"Peso (g)")
kable_fun(tab2, "Tab.4 - Indici di posizione e di variabilità aggiornati sulle variabili quantitative del dataset")| Anni.madre | N.Gravidanze | Settimane Gestazione | Lunghezza(mm) | Diametro craniale(mm) | Peso (g) | |
|---|---|---|---|---|---|---|
| min | 13.00 | 0.00 | 25.00 | 310.00 | 235.00 | 830.00 |
| Q1 | 25.00 | 0.00 | 38.00 | 480.00 | 330.00 | 2990.00 |
| median | 28.00 | 1.00 | 39.00 | 500.00 | 340.00 | 3300.00 |
| mean | 28.19 | 0.98 | 38.98 | 494.70 | 340.03 | 3284.18 |
| Q3 | 32.00 | 1.00 | 40.00 | 510.00 | 350.00 | 3620.00 |
| max | 46.00 | 12.00 | 43.00 | 565.00 | 390.00 | 4930.00 |
| STD_dev | 5.22 | 1.28 | 1.87 | 26.33 | 16.43 | 525.23 |
| CV | 0.19 | 1.31 | 0.05 | 0.05 | 0.05 | 0.16 |
| skewness | 0.15 | 2.51 | -2.07 | -1.51 | -0.79 | -0.65 |
| kurtosis | -0.11 | 10.98 | 8.26 | 6.48 | 2.94 | 2.03 |
| Range | 33.00 | 12.00 | 18.00 | 255.00 | 155.00 | 4100.00 |
| IQR | 7.00 | 1.00 | 2.00 | 30.00 | 20.00 | 630.00 |
LA VARIABILE ANNI.MADRE
il vettore numerico Anni.Madre ha una media di 28,19 anni e dev.standard di +-5.22 anni. La distribuzione presenta una forma leggermente asimmetrica positiva (skewness 0.15) e platicurtica (-0.11), con un coefficiente di variazione del 19%. Il valore più frequente (moda) è 30 anni (con 200 osservazioni registrate), cioè l’8% delle osservazioni del dataset neonati;
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
#nbins <- (1+log2(n))
bins_amp <- 1
#vector_classes <-cut(Anni.madre,(seq(min(Anni.madre),max(Anni.madre),bins_amp)))
ni <- table(Anni.madre)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_Anni.madre<- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_Anni.madre) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
colnames(distr_freq_Anni.madre) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)","Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_Anni.madre, "Tab.5 - Distristribuzione di frequenza - Anni della Madre")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| 13 | 1 | 0.00 | 1 | 0.00 |
| 14 | 2 | 0.00 | 3 | 0.00 |
| 15 | 6 | 0.00 | 9 | 0.00 |
| 16 | 13 | 0.01 | 22 | 0.01 |
| 17 | 18 | 0.01 | 40 | 0.02 |
| 18 | 24 | 0.01 | 64 | 0.03 |
| 19 | 45 | 0.02 | 109 | 0.04 |
| 20 | 66 | 0.03 | 175 | 0.07 |
| 21 | 74 | 0.03 | 249 | 0.10 |
| 22 | 100 | 0.04 | 349 | 0.14 |
| 23 | 115 | 0.05 | 464 | 0.19 |
| 24 | 131 | 0.05 | 595 | 0.24 |
| 25 | 180 | 0.07 | 775 | 0.31 |
| 26 | 184 | 0.07 | 959 | 0.38 |
| 27 | 197 | 0.08 | 1156 | 0.46 |
| 28 | 172 | 0.07 | 1328 | 0.53 |
| 29 | 174 | 0.07 | 1502 | 0.60 |
| 30 | 200 | 0.08 | 1702 | 0.68 |
| 31 | 147 | 0.06 | 1849 | 0.74 |
| 32 | 159 | 0.06 | 2008 | 0.80 |
| 33 | 110 | 0.04 | 2118 | 0.85 |
| 34 | 96 | 0.04 | 2214 | 0.89 |
| 35 | 66 | 0.03 | 2280 | 0.91 |
| 36 | 64 | 0.03 | 2344 | 0.94 |
| 37 | 41 | 0.02 | 2385 | 0.95 |
| 38 | 38 | 0.02 | 2423 | 0.97 |
| 39 | 27 | 0.01 | 2450 | 0.98 |
| 40 | 19 | 0.01 | 2469 | 0.99 |
| 41 | 13 | 0.01 | 2482 | 0.99 |
| 42 | 8 | 0.00 | 2490 | 1.00 |
| 43 | 2 | 0.00 | 2492 | 1.00 |
| 44 | 4 | 0.00 | 2496 | 1.00 |
| 45 | 1 | 0.00 | 2497 | 1.00 |
| 46 | 1 | 0.00 | 2498 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = Anni.madre)) +
geom_histogram(binwidth = bins_amp, fill = "skyblue", color = "white") +
scale_x_continuous(name = "Anni della Madre", breaks = seq(0,60, by=5)) +
scale_y_continuous(name = "Frequenza assoluta", breaks = seq(0, 1200, by = 100)) +
ggtitle("Istogramma di frequenza di Anni.madre") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
boxplot(Anni.madre,
main="Boxplot di Anni.madre",
xlab="Anni madre",
ylab="Distribuzione",
col="skyblue",
border="black")
LA VARIABILE N. DI GRAVIDANZE
Il vettore numerico in questione presenta una media di 0.98+-1.28 gravidanze (dev. standard) e un coefficiente di variazione CV del 130%. Il numero di gravidanze delle madri può essere rappresentato efficacemente come variabile numerica discreta, che varia da 0 (per le madri alla prima gravidanza) fino a un massimo di 12 gravidanze. In accordo con le distribuzioni di frequenza rappresentate in tabella, la maggior parte delle madri considerate è in attesa del primo o del secondo figlio ed, in misura minore, il attesa del terzo (circa il 96% delle osservazioni registrate delle madri è racchiusa in questo range di gravidanze 0-2). La moda è 0 (madri alla prima gravidanza). Graficamente il numero di gravidanze è stato realizzato tramite un grafico a barre, data la natura discreta con pochi valori distinti. L’altezza della barra rappresenta la frequenza assoluta del rispettivo numero di gravidanza.
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
nbins <- (1+log2(n))
bins_amp <- (max(N.gravidanze)-min(N.gravidanze))/nbins
#vector_classes <-cut(N.gravidanze,(seq(min(N.gravidanze),max(N.gravidanze),bins_amp)))
ni <- table(N.gravidanze)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_N.gravidanze <- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_N.gravidanze) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_N.gravidanze, "Tab.6 - Distr. frequenza - Numero di gravidanze delle madri")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| 0 | 1095 | 0.44 | 1095 | 0.44 |
| 1 | 817 | 0.33 | 1912 | 0.77 |
| 2 | 340 | 0.14 | 2252 | 0.90 |
| 3 | 150 | 0.06 | 2402 | 0.96 |
| 4 | 48 | 0.02 | 2450 | 0.98 |
| 5 | 21 | 0.01 | 2471 | 0.99 |
| 6 | 11 | 0.00 | 2482 | 0.99 |
| 7 | 1 | 0.00 | 2483 | 0.99 |
| 8 | 8 | 0.00 | 2491 | 1.00 |
| 9 | 2 | 0.00 | 2493 | 1.00 |
| 10 | 3 | 0.00 | 2496 | 1.00 |
| 11 | 1 | 0.00 | 2497 | 1.00 |
| 12 | 1 | 0.00 | 2498 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = as.factor(N.gravidanze))) +
geom_bar(fill = "lightblue",
alpha = 1,
color = "white") +
scale_y_continuous(name = "Frequenza assoluta", breaks = seq(0, 1200, by = 100))+
labs(title = "Numero di gravidanze \n (grafico a barre)",
x = "numero di gravidanze madre",
y = "frequenza assoluta")+
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
LA VARIABILE GESTAZIONE (SETTIMANE DI GESTAZIONE)
Le settimane di gestazione hanno una media di circa 38.98 +- 1.87 settimane (dev standard), e un coefficiente di variazione CV del 5%.La forma della distribuzione è caratterizzata da asimmetria negativa e leptocurtosi. Questa asimmetria negativa può essere indice di presenza di casi sporadici di durate di gravidanze significativemente più corte rispetto alla media delle osservazioni. Il valore di moda è 40 settimane con il 30% di osservazione totali nel dataset, indicando una maggiore frequenza assolute di parti registrati in un tempo di gestazione pari a 40 settimane.
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
nbins <- (1+log2(n))
bins_amp <- 1
#vector_classes <-cut(Gestazione,(seq(min(Gestazione),max(Gestazione),bins_amp)))
ni <- table(Gestazione)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_Gestazione <- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_Gestazione) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_Gestazione, "Tab.7 - Distristribuzione di frequenza - Settimane di Gestazione")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| 25 | 1 | 0.00 | 1 | 0.00 |
| 26 | 1 | 0.00 | 2 | 0.00 |
| 27 | 2 | 0.00 | 4 | 0.00 |
| 28 | 4 | 0.00 | 8 | 0.00 |
| 29 | 3 | 0.00 | 11 | 0.00 |
| 30 | 5 | 0.00 | 16 | 0.01 |
| 31 | 8 | 0.00 | 24 | 0.01 |
| 32 | 9 | 0.00 | 33 | 0.01 |
| 33 | 18 | 0.01 | 51 | 0.02 |
| 34 | 16 | 0.01 | 67 | 0.03 |
| 35 | 33 | 0.01 | 100 | 0.04 |
| 36 | 62 | 0.02 | 162 | 0.06 |
| 37 | 192 | 0.08 | 354 | 0.14 |
| 38 | 437 | 0.17 | 791 | 0.32 |
| 39 | 580 | 0.23 | 1371 | 0.55 |
| 40 | 741 | 0.30 | 2112 | 0.85 |
| 41 | 328 | 0.13 | 2440 | 0.98 |
| 42 | 56 | 0.02 | 2496 | 1.00 |
| 43 | 2 | 0.00 | 2498 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = Gestazione)) +
geom_histogram(binwidth = 1, fill = "skyblue", color = "white") +
scale_x_continuous(name = "Gestazione (settimane)", breaks = seq(0,60, by=2)) +
scale_y_continuous(name = "Frequenza assoluta", breaks = seq(0, 1400, by = 200)) +
ggtitle("Istogramma di frequenza \n delle settimane di Gestazione") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
boxplot(Gestazione,
main="Boxplot di Gestazione",
xlab="Gestazione",
ylab="Distribuzione",
col="skyblue",
border="black")
LA VARIABILE LUNGHEZZA (LUNGHEZZA NEONATO)
La lunghezza media dei neonati è pari a 494,70 mm (deviazione standard ±26,33 mm) , e un coefficiente di variazione CV= 5%. La distribuzione è caratterizzata fa forma asimmetrica negativa e leptocurtosi. La Classe modale con maggiore frequenza assoluta è quella della lunghezza (497,518] millimetri (518mm escluso). Anche qui sono presenti diversi neonati la cui lunghezza registrata tramite controllo ecografico è abbastanza inferiore alla media, e tali osservazioni spostano le code delle distribuzioni verso valori inferiori.
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
nbins <- (1+log2(n))
bins_amp <- (max(Lunghezza)-min(Lunghezza))/nbins
vector_classes <-cut(Lunghezza,(seq(min(Lunghezza),max(Lunghezza),bins_amp)))
ni <- table(vector_classes)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_Lunghezza <- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_Lunghezza) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_Lunghezza, "Tab.8 - Distribuzione di Frequenza - Lunghezza Neonato")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| (310,331] | 3 | 0.00 | 3 | 0.00 |
| (331,352] | 2 | 0.00 | 5 | 0.00 |
| (352,372] | 7 | 0.00 | 12 | 0.00 |
| (372,393] | 9 | 0.00 | 21 | 0.01 |
| (393,414] | 15 | 0.01 | 36 | 0.01 |
| (414,435] | 19 | 0.01 | 55 | 0.02 |
| (435,455] | 89 | 0.04 | 144 | 0.06 |
| (455,476] | 313 | 0.13 | 457 | 0.18 |
| (476,497] | 737 | 0.30 | 1194 | 0.48 |
| (497,518] | 886 | 0.35 | 2080 | 0.83 |
| (518,538] | 344 | 0.14 | 2424 | 0.97 |
| (538,559] | 70 | 0.03 | 2494 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = Lunghezza)) +
geom_histogram(binwidth = bins_amp, fill = "skyblue", color = "white") +
scale_x_continuous(name = "Lunghezza neonato (mm)", breaks = seq(0,800, by=50)) +
scale_y_continuous(name = "Frequenza", breaks = seq(0, 1200, by = 100)) +
ggtitle("Istogramma di frequenza di Lunghezza neonato") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#boxplot
boxplot(Lunghezza,
main="Lunghezza",
xlab="Gestazione",
ylab="Distribuzione",
col="skyblue",
border="black")
LA VARIABILE CRANIO
Il diametro craniale medio è di 340,03mm (deviazione standard ±16,43 mm), e un coefficiente di variazione CV del 5%.La distribuzione è caratterizzata fa forma asimmetrica negativa e leptocurtosi. La Classe modale con maggiore frequenza assoluta del diametro craniale è (336-349] millimetri (349mm escluso). Anche qui sono presenti diversi neonati il cui diametro craniale registrato tramite controllo ecografico è abbastanza inferiore alla media, e tali osservazioni spostano le code delle distribuzioni verso valori inferiori.
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
nbins <- (1+log2(n))
bins_amp <- (max(Cranio)-min(Cranio))/nbins
vector_classes <-cut(Cranio,(seq(min(Cranio),max(Cranio),bins_amp)))
ni <- table(vector_classes)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_Cranio <- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_Cranio) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_Cranio, "Tab.7 - Distr.Frequenza - Cranio Neonato")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| (235,248] | 1 | 0.00 | 1 | 0.00 |
| (248,260] | 2 | 0.00 | 3 | 0.00 |
| (260,273] | 5 | 0.00 | 8 | 0.00 |
| (273,285] | 12 | 0.00 | 20 | 0.01 |
| (285,298] | 15 | 0.01 | 35 | 0.01 |
| (298,311] | 73 | 0.03 | 108 | 0.04 |
| (311,323] | 198 | 0.08 | 306 | 0.12 |
| (323,336] | 629 | 0.25 | 935 | 0.37 |
| (336,349] | 811 | 0.32 | 1746 | 0.70 |
| (349,361] | 575 | 0.23 | 2321 | 0.93 |
| (361,374] | 144 | 0.06 | 2465 | 0.99 |
| (374,386] | 28 | 0.01 | 2493 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = Cranio)) +
geom_histogram(binwidth = bins_amp, fill = "skyblue", color = "white") +
scale_x_continuous(name = "Cranio neonato (mm)", breaks = seq(0,800, by=25)) +
scale_y_continuous(name = "Frequenza", breaks = seq(0, 1200, by = 100)) +
ggtitle("Istogramma di frequenza del diametro craniale") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
boxplot(Cranio,
main="Boxplot di Cranio",
xlab="Gestazione",
ylab="Distribuzione",
col="skyblue",
border="black")
LA VARIABILE PESO (TARGET)
graficamente la varaibile target Peso somiglia vagamente ad una distribuzione normale, ma il test di Shapiro-wilk indica che per un p-value inferiore a 2.2e-16 viene rigettata l’ipotesi H0 nulla di normalità distributiva della variabile a favore di un’ipotesi alternativa H1 di distribuzione NON normale. I parametri di forma skewness e curtosi (rispettivamente -0.64 e 2.0) che ci indicano che la distribuzione di tipo leptocurtica ed asimmetrica negativa , probabilmente a causa della presenza di molti outliers verso in valori inferiori della distribuzione.
La variabile Peso ha una media di Peso: media 3284.18 +- 525.23g (dev standard) ed un Coeff.di variazione CV del 16%. Il range della variabile oscilla fra valori 830g e 4930g, quindi un range di 4100g molto elevato. La classe modale della variabile peso è (3170-3500] grammi (3500g escluso).
#CALCOLO DISTRIBUZIONI DI FREQUENZA
n <- nrow(dati)
nbins <- (1+log2(n))
bins_amp <- (max(Peso)-min(Peso))/nbins
vector_classes <-cut(Peso,(seq(min(Peso),max(Peso),bins_amp)))
ni <- table(vector_classes)
fi <- ni / n
Ni <- cumsum(ni)
Fi <- cumsum(ni) / n
distr_freq_Peso <- round(as.data.frame(cbind(ni, fi, Ni, Fi)), 2)
colnames(distr_freq_Peso) <- c("Freq.assoluta (ni)", "Freq.relativa (fi)",
"Freq.cumulata (Ni)", "Freq.cum.relativa (Fi)")
#tabella
kable_fun(distr_freq_Peso, "Tab.9 - Distribuzione di Frequenza - Variabile Target Peso del neonato")| Freq.assoluta (ni) | Freq.relativa (fi) | Freq.cumulata (Ni) | Freq.cum.relativa (Fi) | |
|---|---|---|---|---|
| (830,1.16e+03] | 6 | 0.00 | 6 | 0.00 |
| (1.16e+03,1.5e+03] | 15 | 0.01 | 21 | 0.01 |
| (1.5e+03,1.83e+03] | 18 | 0.01 | 39 | 0.02 |
| (1.83e+03,2.16e+03] | 31 | 0.01 | 70 | 0.03 |
| (2.16e+03,2.5e+03] | 66 | 0.03 | 136 | 0.05 |
| (2.5e+03,2.83e+03] | 248 | 0.10 | 384 | 0.15 |
| (2.83e+03,3.17e+03] | 566 | 0.23 | 950 | 0.38 |
| (3.17e+03,3.5e+03] | 675 | 0.27 | 1625 | 0.65 |
| (3.5e+03,3.83e+03] | 545 | 0.22 | 2170 | 0.87 |
| (3.83e+03,4.17e+03] | 242 | 0.10 | 2412 | 0.97 |
| (4.17e+03,4.5e+03] | 67 | 0.03 | 2479 | 0.99 |
| (4.5e+03,4.83e+03] | 16 | 0.01 | 2495 | 1.00 |
#istogramma di frequenza
ggplot(dati, aes(x = Peso)) +
geom_histogram(binwidth = bins_amp, fill = "skyblue", color = "white") +
scale_x_continuous(name = "Peso neonato (grammi)", breaks = seq(0,5000, by=500)) +
scale_y_continuous(name = "Frequenza", breaks = seq(0, 1200, by = 100)) +
ggtitle("Istogramma di frequenza \n variabile Peso (target)") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#plot variabile Peso vs Teorica Normale
dens_Peso <- density(Peso)
plot(dens_Peso,
main = "Kernel Density Peso \ vs Normale Teorica",
xlab = "Distribuzione dei valori Peso",
ylab = "Densità", col = "blue", lwd = 2)
#dati della curva teorica
mean_Peso <- mean(Peso)
sd_Peso <- sd(Peso)
curve(dnorm(x,
mean = mean_Peso,
sd = sd_Peso),
col = "black",
lwd = 5,
add = TRUE)
legend("topleft",
legend = c("Distribuzione Variabile Peso",
"Normale teorica"),
col = c("blue", "black"), lwd = 2)
#boxplot Peso
boxplot(Peso,
xlab="Peso neonati (g)",
ylab="Distribuzione",
main="Boxplot (Peso)",
col="skyblue",
border="black")
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97068, p-value < 2.2e-16- OUTLIERS variabile target Peso
outliers_analysis <- as.data.frame(cbind(Anni.madre, Gestazione,Peso, Lunghezza,Cranio, N.gravidanze))
#definizione funzione outliers lowerbound
identify_outliers_lower <- function(df, col) {
q1 <- quantile(df[[col]], 0.25)
q3 <- quantile(df[[col]], 0.75)
iqr <- IQR(df[[col]])
outs_lower <- df[df[, col] < q1 - 1.5 * iqr,]
return(outs_lower)}
#definizione funzione outliers upperbound
identify_outliers_upper <- function(df, col) {
q1 <- quantile(df[[col]], 0.25)
q3 <- quantile(df[[col]], 0.75)
iqr <- IQR(df[[col]])
outs_upper <- df[df[, col] > q3 + 1.5 * iqr,]
return(outs_upper)}
#INDAGINE
kable_fun(identify_outliers_lower(outliers_analysis,"Peso")," Tab.10 - outliers inferiori della variabile Peso ")| Anni.madre | Gestazione | Peso | Lunghezza | Cranio | N.gravidanze | |
|---|---|---|---|---|---|---|
| 101 | 31 | 34 | 1370 | 390 | 287 | 0 |
| 106 | 29 | 30 | 1340 | 400 | 273 | 4 |
| 206 | 39 | 31 | 1500 | 405 | 295 | 1 |
| 295 | 18 | 40 | 1850 | 460 | 305 | 0 |
| 310 | 40 | 28 | 1560 | 420 | 379 | 3 |
| 312 | 26 | 32 | 1280 | 360 | 276 | 1 |
| 322 | 25 | 37 | 1750 | 430 | 305 | 1 |
| 378 | 27 | 28 | 1285 | 400 | 274 | 0 |
| 445 | 27 | 32 | 1550 | 410 | 289 | 0 |
| 492 | 34 | 33 | 1410 | 380 | 295 | 2 |
| 587 | 16 | 31 | 1900 | 440 | 300 | 1 |
| 638 | 25 | 33 | 1720 | 420 | 300 | 0 |
| 724 | 39 | 35 | 1980 | 450 | 308 | 2 |
| 748 | 35 | 33 | 1390 | 390 | 277 | 0 |
| 750 | 24 | 35 | 1450 | 405 | 280 | 0 |
| 765 | 26 | 33 | 1970 | 445 | 300 | 1 |
| 805 | 30 | 29 | 1190 | 360 | 272 | 2 |
| 838 | 32 | 35 | 2000 | 430 | 312 | 0 |
| 928 | 25 | 28 | 830 | 310 | 254 | 0 |
| 947 | 34 | 32 | 1615 | 390 | 297 | 3 |
| 1067 | 26 | 31 | 1960 | 420 | 300 | 3 |
| 1091 | 30 | 33 | 1770 | 410 | 275 | 1 |
| 1096 | 37 | 34 | 1750 | 420 | 306 | 0 |
| 1247 | 26 | 30 | 1170 | 370 | 266 | 1 |
| 1272 | 32 | 33 | 2040 | 480 | 307 | 1 |
| 1274 | 24 | 38 | 1980 | 430 | 305 | 0 |
| 1310 | 40 | 34 | 1840 | 430 | 305 | 6 |
| 1383 | 33 | 29 | 1620 | 410 | 292 | 0 |
| 1426 | 30 | 36 | 1280 | 385 | 292 | 1 |
| 1427 | 24 | 29 | 1280 | 390 | 355 | 4 |
| 1508 | 31 | 37 | 2040 | 465 | 305 | 0 |
| 1591 | 41 | 35 | 1500 | 420 | 304 | 3 |
| 1608 | 37 | 33 | 2000 | 470 | 293 | 3 |
| 1617 | 31 | 31 | 990 | 340 | 278 | 0 |
| 1684 | 27 | 31 | 1800 | 430 | 308 | 0 |
| 1699 | 22 | 32 | 1430 | 380 | 301 | 0 |
| 1741 | 19 | 38 | 1950 | 445 | 305 | 1 |
| 1753 | 34 | 36 | 1970 | 450 | 303 | 0 |
| 1778 | 25 | 25 | 900 | 325 | 253 | 2 |
| 1807 | 35 | 32 | 1780 | 420 | 277 | 0 |
| 2087 | 32 | 33 | 1780 | 400 | 305 | 1 |
| 2112 | 36 | 31 | 1180 | 355 | 270 | 0 |
| 2113 | 35 | 32 | 1890 | 500 | 309 | 1 |
| 2118 | 32 | 27 | 1140 | 370 | 267 | 0 |
| 2138 | 30 | 33 | 1600 | 410 | 290 | 2 |
| 2147 | 39 | 30 | 1300 | 380 | 276 | 3 |
| 2173 | 37 | 28 | 930 | 355 | 235 | 8 |
| 2198 | 33 | 30 | 1750 | 410 | 294 | 0 |
| 2222 | 41 | 33 | 2000 | 425 | 312 | 1 |
| 2255 | 40 | 34 | 1580 | 400 | 300 | 0 |
| 2305 | 26 | 30 | 1170 | 370 | 273 | 1 |
| 2406 | 37 | 31 | 1690 | 405 | 290 | 2 |
| 2435 | 28 | 27 | 980 | 320 | 265 | 1 |
| 2450 | 28 | 26 | 930 | 345 | 245 | 0 |
| 2456 | 31 | 31 | 1730 | 430 | 300 | 0 |
kable_fun(identify_outliers_upper(outliers_analysis,"Peso"),"Tab.11 - outliers superiori della variabile Peso ")| Anni.madre | Gestazione | Peso | Lunghezza | Cranio | N.gravidanze | |
|---|---|---|---|---|---|---|
| 126 | 27 | 41 | 4680 | 545 | 370 | 0 |
| 368 | 30 | 41 | 4600 | 540 | 370 | 0 |
| 1292 | 30 | 38 | 4600 | 485 | 380 | 3 |
| 1305 | 23 | 41 | 4900 | 510 | 352 | 0 |
| 1431 | 31 | 41 | 4810 | 530 | 364 | 4 |
| 1511 | 25 | 40 | 4620 | 560 | 367 | 0 |
| 1637 | 39 | 40 | 4760 | 550 | 365 | 3 |
| 1836 | 34 | 40 | 4580 | 515 | 360 | 3 |
| 1918 | 26 | 39 | 4930 | 550 | 350 | 0 |
| 1961 | 27 | 42 | 4700 | 540 | 362 | 0 |
| 2021 | 27 | 39 | 4650 | 510 | 354 | 1 |
| 2074 | 23 | 40 | 4720 | 540 | 360 | 0 |
| 2233 | 30 | 41 | 4690 | 540 | 373 | 2 |
| 2390 | 28 | 40 | 4720 | 540 | 355 | 1 |
E’ da notare che una grande quantità di outliers inferiori al lowerbound (Q1-1.5IQR) della distribuzione Peso riportati in Tab.10, siano relativi al peso di neonati che difficilmente hanno superato le 35 settimane di gestazione in fase di gravidanza della madre.
1.2.4 VARIABILI QUALITATIVE - Frequenze ed Indice di Gini delle (Fumatrici, Tip.parto, Ospedale, Sesso)
A seguire, riporterò le analisi sulle variabili qualitative del dataset e commenterò i risultati:
n <- nrow(dati) #nuove codficihe Fumatrici_categoriale <- factor(Fumatrici, levels = c(0, 1), labels = c("NON fumatrice", "Fumatrice")) parto_categoriale <- factor(Tipo.parto, levels = c("Ces", "Nat"), labels = c("Parto Cesario", "Parto Naturale")) Sesso_categoriale <- factor(Sesso, levels = c("M", "F"), labels = c("Maschio", "Femmina")) #Fumatrici ni <- table(Fumatrici_categoriale) fi <- ni/n*100 ni_sum <- sum(ni) fi_sum <- sum(fi) sum_cbind<- as.data.frame(cbind(ni_sum,fi_sum)) rownames(sum_cbind) <- c("TOTALE CONTEGGIO FUMATRICE") colnames(sum_cbind) <- c("ni", "fi") distr_freq1 <- round((as.data.frame(cbind(ni,fi))),2) distr_freq1 <- as.data.frame(rbind(distr_freq1, sum_cbind)) #Tipo.parto ni <- table(parto_categoriale) fi <- ni/n*100 distr_freq2 <- round((as.data.frame(cbind(ni,fi))),2) ni_sum <- sum(ni) fi_sum <- sum(fi) sum_cbind<- as.data.frame(cbind(ni_sum,fi_sum)) rownames(sum_cbind) <- c("TOTALE CONTEGGIO PARTO") colnames(sum_cbind) <- c("ni", "fi") distr_freq2 <- round((as.data.frame(cbind(ni,fi))),2) distr_freq2 <- as.data.frame(rbind(distr_freq2, sum_cbind)) #Ospedale ni <- table(Ospedale) fi <- ni/n*100 ni_sum <- sum(ni) fi_sum <- sum(fi) sum_cbind<- as.data.frame(cbind(ni_sum,fi_sum)) rownames(sum_cbind) <- c("TOTALE CONTEGGIO OSPEDALE") colnames(sum_cbind) <- c("ni", "fi") distr_freq3 <- round((as.data.frame(cbind(ni,fi))),2) distr_freq3 <- as.data.frame(rbind(distr_freq3, sum_cbind)) #Sesso ni <- table(Sesso_categoriale) fi <- ni/n*100 ni_sum <- sum(ni) fi_sum <- sum(fi) sum_cbind<- as.data.frame(cbind(ni_sum,fi_sum)) rownames(sum_cbind) <- c("TOTALE CONTEGGIO SESSO") colnames(sum_cbind) <- c("ni", "fi") distr_freq4 <- round((as.data.frame(cbind(ni,fi))),2) distr_freq4 <- as.data.frame(rbind(distr_freq4, sum_cbind)) tabella_unica <- rbind(distr_freq1,distr_freq2,distr_freq3,distr_freq4) kable_fun(tabella_unica, "Tab.12 - Riassuntiva delle frequenze assolute (ni) e relative (fi) delle variabili categoriali")Tab.12 - Riassuntiva delle frequenze assolute (ni) e relative (fi) delle variabili categoriali
ni
fi
NON fumatrice
2394
95.84
Fumatrice
104
4.16
TOTALE CONTEGGIO FUMATRICE
2498
100.00
Parto Cesario
728
29.14
Parto Naturale
1770
70.86
TOTALE CONTEGGIO PARTO
2498
100.00
osp1
816
32.67
osp2
848
33.95
osp3
834
33.39
TOTALE CONTEGGIO OSPEDALE
2498
100.00
Maschio
1243
49.76
Femmina
1255
50.24
TOTALE CONTEGGIO SESSO
2498
100.00
#INDICE DI ETEROGENEITA' di GINI: #per variabili qualitative gini.index <- function(x){ ni=table(x) fi=ni/length(x) fi2=fi^2 J=length(table(x)) gini= 1-sum(fi2) gini.normalizzato=round(gini/((J-1)/J),2) return(gini.normalizzato) } Gini_df <- as.data.frame(cbind(Fumatrici,Tipo.parto, Sesso,Ospedale)) Gini_df<- (sapply(Gini_df,gini.index)) Gini_df <- as.data.frame(Gini_df) colnames(Gini_df) <- c("Indice di Gini") kable_fun(Gini_df, "Tab.13 - Indici di Gini per le variabili qualitative categoriche")Tab.13 - Indici di Gini per le variabili qualitative categoriche
Indice di Gini
Fumatrici
0.16
Tipo.parto
0.83
Sesso
1.00
Ospedale
1.00
#grafici a barre freq.N.fumatrici <- as.data.frame(table(Fumatrici)) grafico1 <- ggplot(data=freq.N.fumatrici, aes(x=Fumatrici, y=Freq))+ geom_bar(stat = "identity", fill="green", color="black")+ labs(title="frequenza madri fumatrici", x="Fumatrice", y="Frequenze assolute")+ theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) freq.tipo.parto <-as.data.frame(table(Tipo.parto)) grafico2 <- ggplot(data=freq.tipo.parto, aes(x=Tipo.parto,y=Freq))+ geom_bar(stat="identity", fill="blue", color="black")+ labs(title="Frequenza parti suddivisi per \ntipologia", x="tipo parto", y="Frequenze assolute")+ theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) tipo.ospedale <-as.data.frame(table(Ospedale)) grafico3 <- ggplot(data=tipo.ospedale, aes(x=Ospedale,y=Freq))+ geom_bar(stat="identity", fill="yellow", color="black")+ labs(title="Tipo Ospedale", x="Unità Ospedaliera", y="Frequenze assolute")+ theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) freq.Sesso <-as.data.frame(table(Sesso)) grafico4 <- ggplot(data=freq.Sesso, aes(x=Sesso,y=Freq))+ geom_bar(stat="identity", fill="orange", color="black")+ labs(title="Frequenza Sesso", x="Sesso Nascituro", y="Frequenze assolute")+ theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) grid.arrange(grafico1, grafico2, grafico3, grafico4)
Moda Fumatrici= 0 (madre non fumatrice), con 2396 osservazioni che corrispondono al 96% dell’intero dataset. Il 4% restante è costituito da madri non fumatrici codificate con Fumatrice=1;
Moda Tipologia Parto = Parto naturale con 1773 osservazioni, circa il 71% delle osservazioni nel dataset), mentre il restante 29% del dataset è costituito da parti di tipo Cesareo;
Moda Tipo.Ospedale = Ospedale 2 con 849 osservazione (34% delle osservazione), mentre il resto delle osservazioni sono equi-distribuite con ospedale fra ospedale 1 e ospedale 3
Sesso neonato = variabile categorica equi-distribuita al 50% fra i 2 sessi del neonati.
Calcolo dell’INDICE DI ETEROGENEITA’ DI GINI per le variabili categoriche:
L’indice di Gini può assumere i seguenti valori:
G = 0, massima concentrazione, corrispondente ad una distribuzione statistica eterodistribuita, con MINIMA eterogeneità fra i valori assunti da una variabile
G =1, minima concentrazione, corrispondente ad una distribuzione statististica equidistribuita, con la MASSIMA massima eterogeneità dei valori
In accordo con quanto sopra, la variabile Fumatrici presentando valori di Gini prossimi allo 0 (0.159), è costituta da un’alta concentrazione di una delle 2 modalità possibili che, come visto dai grafici precedenti è la modalità 0 (madre NON fumatrici).
La variabile con indice di Gini medio-alto è il vettore Tipo.parto (Gini=0.82) che risulta essere eterodistribuita per il 70% circa nella modalità parto Naturale.
Infine i sesso del neonato con il valore di Gini più elevato, e prossimo a 1, indicando uniformità nella distribuzione delle osservazioni.
1.3 - Valutazione delle pratiche ospedaliere
Attraverso un’analisi comparativa tra i tre ospedali coinvolti, valuteremo se vi è una maggiore incidenza di parti cesarei in una determinata struttura. Questo consente di monitorare le qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.
Per testare l’ipotesi nulla che la distribuzione dei cesarei sia uguale nei tre ospedali presi in considerazione si utilizza un test del Chi-quadro.
Dalla tabella 12, sopra riportata, notiamo che le percentuali di frequenza relativa ai parti cesarei e natuali condizionatamente all’ospedale sono molto simili nelle 3 unità ospedaliera . Per avvalorare questa nostra prima analisi eseguiremo un test statistica del chi-quadrato.
Test del chi-quadrato
Un test che potremmo usare che saggiare una generale tendenza all’associazione fra le variabili Ospedale e tipo parto è quello del CHI-quadrato, eseguito sulla base di tabelle di contingenza distribuzioni statistiche di tipo chi-quadro.
Per eseguire il test si raggruppano per prima le frequenze assolute di osservazioni dei pesi in funzione dell’età della madre, e poi si saggia se esiste associazione o indipendenza fra il peso è gli anni della madre. Il test assume H0 ipotesi nulla come indipendenza fra le variabili (assenza di associazione), di converso l’ipotesi H1 alternativa riguarda l’associazione fra le variabili.
tabella <- table(Ospedale,Tipo.parto)
tabella_perc <- round(prop.table(tabella,margin =1) * 100,2)
tabella_perc <- as.data.frame(tabella_perc)
tabella_wide <- pivot_wider(
data = tabella_perc,
names_from = Tipo.parto,
values_from = Freq
)
colnames(tabella_wide) <- c("Unità ospedaliera", "Cesario", " Naturale")
kable_fun(tabella_wide, "Tab.13 - Distribuzione(%) \n tipologie parto nelle 3 unità ospedaliere")| Unità ospedaliera | Cesario | Naturale |
|---|---|---|
| osp1 | 29.66 | 70.34 |
| osp2 | 29.95 | 70.05 |
| osp3 | 27.82 | 72.18 |
#test
chisq.test(tabella)##
## Pearson's Chi-squared test
##
## data: tabella
## X-squared = 1.083, df = 2, p-value = 0.5819Per un valore p di circa 0,58, superiore al livello di significatività (in genere posto allo alpha=0.01), il test chi-quadro appena eseguito non rigetta l’ipotesi nulla (H₀) di indipendenza. Questo implica l’assenza di un’associazione significativa tra le tipologie di parto registrate nelle diverse aziende ospedaliere. In altre parole, il test conferma quanto osservato nella tabella 13 di riepilogo : non vi sono evidenze significative che suggeriscano che determinate tipologie di parto siano associate in modo specifico alle diverse unità ospedaliere.
1.4 - Variabili materne (età,fumo,n.gravidanze) vs Peso nel neonato
A questo punto il focus sarà rivolto alla relazione tra le variabili materne (età, fumo, gravidanze precedenti) e il peso del neonato. Produrremo scatterplot bidimensionali al fine di cercare relazioni grafiche all’interno delle distribuzioni del dataset, ed eseguiremo i test di correlazioni con le variabili target.
1.4.1. Peso del neonato ~ Anni della madre
In questa rassegna di test statistici andremo ad indagare se vi è una relazione di dipendenza fra il peso del neonato e gli anni della madre. Per fare questo tipo possiamo iniziare ad indagare gli scatterplot bidimensionali delle due variabili e visualizzarne la relazione, per poi passare all’analisi numerica tramite il coefficiente di correlazione lineare di Bravais-Pearson, ed infine i relativi:
plot(Anni.madre, Peso,
main="Scatterplot \nAnni.madre vs Peso neoanati",
xlab="Anni madre", ylab="Pesi neonati")
Dallo scatterplot bidimensionale NON si visualizzano particolari correlazioni fra la variabile Anni.madre (indipendente) e la potenziale variabile target Peso del neonato (dipendente ), indicando una netta potenziale decorrelazione fra le variabili, almeno visiva.
Continuando con la disamina numerica, cerchiamo di utilizzare uno strumento numerico statistico come il coefficiente di correlazione lineare di Bravais-Pearson, fra le variabili Anni.madre (indipendente) e Peso del neonato (target).
Volendo sintetizzare, a parole, cos’è il coefficiente di Pearson, quest’ultimo è un parametro numerico standardizza la covarianza fra 2 variabili numeriche, nel campo numerico fra -1 e 1 (0 compreso), ed indica la concordanza o discordanza fra due variabili quantitative, e cioè se le 2 variabili variano numericamente contemporanemente o meno . Nel particolare se il coefficiente di correlazione è pari a 0 le due variabili sono decorrelate e i loro valori possono variare in modo indipendente l’una dall’altra; man mano che il valore di coefficiente di correlazione delle variabili si avvicina in valore assoluto a 1, le due variabili considerate variano insieme in modo concorde (aumenta l’una e quindi aumenta l’altra, con r positivo) o discorde (aumenta l’una e quindi diminuisce anche l’altra, con r negativo), in modo sempre più marcato.
cat("coefficiente di Pearson: ", cor(Peso,Anni.madre))## coefficiente di Pearson: -0.02378138Di seguito viene effettuato anche il test di Pearson che valuta se la correlazione lineare tra due variabili quantitative è significativamente diversa da 0. Si basa sul calcolo del coefficiente di correlazione di Pearson, e ha formule le seguenti ipotesi:
Ipotesi nulla (H0): non c’è correlazione tra le variabili (r=0)
Ipotesi alternativa (H1): esiste una correlazione (r diverso da 0).
confronta i risultati della statistica t, con la relativa distribuzione e ne calcola il p-valore (per un predeterminato intervallo di confidenza alpha):
- Se p<alpha (tipicamente 0.05), si rifiuta H0 e si conclude che la correlazione è significativa.
test_pearson <- cor.test(Anni.madre, Peso, method = "pearson")
test_pearson##
## Pearson's product-moment correlation
##
## data: Anni.madre and Peso
## t = -1.1885, df = 2496, p-value = 0.2348
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06294109 0.01545144
## sample estimates:
## cor
## -0.02378138In questo caso il pvalue è del 23.48%, quindi non possiamo rifiutare H0, e ne consegue che Anni.madre e Peso non sono significativamente correlate in modo lineare.
1.4.2 - Peso del neonato ~ Numero di Gravidanze
Indaghiamo le relazioni fra le variabile indipendente Numero di gravidanze e la variabile di riposta peso neonato:
plot(N.gravidanze, Peso, main="Scatterplot \nN.gravidanze vs Peso neoanati",
xlab="Gravidanze Madre",
ylab="Peso Neonato (g)")
#correlazioni lineare
cor(Peso, N.gravidanze)## [1] 0.002277118test_pearson <- cor.test(N.gravidanze, Peso, method = "pearson")
test_pearson##
## Pearson's product-moment correlation
##
## data: N.gravidanze and Peso
## t = 0.11377, df = 2496, p-value = 0.9094
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.03694459 0.04149182
## sample estimates:
## cor
## 0.002277118Anche in questo caso, studiando la relazione Numero di gravidanze pregresse della madre e peso si conclude che non è presente alcuna correlazione lineare , dato un alto valore di p-value di 0.91 relativo al test di Pearson, e un coefficiente di correlazione lineare di 0.0022.
1.4.3 Peso del neonato ~ Fumatrici
Indagheremo con il prossimo test, come le abitudini delle madri legate al fumo, possano influenzare il peso previsto dei neonati, o se più in generale riusciamo a trovare associazione sulla base dei dati ottenuti. Iniziamo con l’eseguire un test t fra le medie dei pesi di neonati, calcolando le medie raggruppate per madre fumatrice (0) e madre NON fumatrice(1):
#analisi grafica
boxplot(Peso~Fumatrici,
main="Boxplot - Distribuzione Pesi Neonati \ncondizionati alle abitudini della madre (Fumo)",
xlab="(0) NON FUMATRICE, (1) FUMATRICE ",
ylab="Peso Neonato (g)")
#analisi numerica
df_fumatrici<-as.data.frame(cbind(NO_FUMATRICE=mean(Peso[Fumatrici==0]),
SI_FUMATRICE=mean(Peso[Fumatrici==1])))
rownames(df_fumatrici) <- c("Peso(grammi)")
#visualizzazione
kable_fun(df_fumatrici,"Tab.13 - \nMedie dei pesi dei neonati - confronto da Madre Fumatrice e NON Fumatrice")| NO_FUMATRICE | SI_FUMATRICE | |
|---|---|---|
| Peso(grammi) | 3286.262 | 3236.346 |
Risultati:
Medie peso neonati da madri non fumatrici: 3286.153 grammi.
Medie peso neonati da madri fumatrici: 3236.346 grammi.
Differenza media stimata: circa -50 grammi (3236g - 3286g).
#inferenza
t.test(Peso~Fumatrici)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.0362, df = 114.12, p-value = 0.3023
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.5076 145.3399
## sample estimates:
## mean in group 0 mean in group 1
## 3286.262 3236.346alpha=0.05
Intervallo di confidenza (95%): Da -45.61354 a +145.22674.
p-value = 0.3033
Con un intervallo di confidenza del (95%): Da -45.61354 a 145.22674, e un p-valore del 0.3033, possiamo affermare che dai dati ottenuti non possiamo rifiutare l’ipotesi nulla per la quale le differenze fra le due medie sia pari a 0, ed in base all’intervallo di confidenza delle (per alpha 0.05), il p-valore indica una probabilità di circa il 30% di osservare una differenza tra i gruppi almeno grande quanto quella trovata, assumendo che la differenza reale sia zero. In altre parole la differenze fra le due medie non è significativamente diversa da 0.
1.5 - Variabili di controllo neonato (Sesso, Lunghezza, Cranio) vs Peso Neonato
1.5.1 - Peso neonato vs Lunghezza e Cranio
Verranno indagate di seguito le relazioni di numeriche, grafiche e di inferenza statistica fra le variabili indipendenti registrate tramite ecografia pre-natale (LUNGHEZZA o CRANIO) e la variabile target PESO:
plot(Lunghezza, Peso,
main="Scatterplot\n Lunghezza Ecografica vs Peso Neonato",
xlab="Lunghezza neonato in (millimetri)") 
plot(Cranio,Peso,
main="Scatterplot\n Cranio vs Peso Neonato",
xlab="Diametro craniale (millimetri)") 
test_pearson_LvsPeso<-cor.test(Peso, Lunghezza, method = "pearson")
test_pearson_LvsPeso ##
## Pearson's product-moment correlation
##
## data: Peso and Lunghezza
## t = 65.71, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7812121 0.8099730
## sample estimates:
## cor
## 0.7960415test_pearson_CvsPeso<-cor.test(Peso,Cranio, method = "pearson")
test_pearson_CvsPeso##
## Pearson's product-moment correlation
##
## data: Peso and Cranio
## t = 49.642, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6845483 0.7240475
## sample estimates:
## cor
## 0.7048438Appare chiara e significativa (p-value<2.2e-16) la relazione lineare fra Peso e le variabili di controllo ecografiche del neonato prima del parto:
Correlazione lineare positiva Lunghezza vs Peso: 0.796
Correlazione lineare positiva Cranio vs Peso: 0.705
indicando che all’aumentare di una delle variabili indipendenti fra la lunghezza e il cranio del neonato, aumenta anche il peso del neonato (variabile target).
1.5.1 Peso del neonato ~ Sesso Neonato
La variabile Sesso è sicuramente un fattore di controllo importante, e probabilmente fornirà informazioni interessanti sull’andamento dei peso dei neonati. Iniziamo con eseguire un boxplot delle distribuzioni dei pesi neonati condizionati ai 2 sessi dei neonati, e nel frattempo calcoliamo la differenza numerica fra le medie dei 2 sessi:
#analisi grafica
boxplot(Peso~Sesso) 
#analisi numerica
df_sex<-as.data.frame(cbind(Maschio=mean(Peso[Sesso=="M"]),
Femmina=mean(Peso[Sesso=="F"])))
rownames(df_sex) <- c("Peso(grammi)")
#visualizzazione
kable_fun(df_sex,"Tab.14 - Medie pesi dei neonati - confronto sulla base del sesso del nascituro")| Maschio | Femmina | |
|---|---|---|
| Peso(grammi) | 3408.496 | 3161.061 |
E’ possibile notare che vi è una differenza numerica fra le medie dei pesi dei neonati, quando questi vengono raggruppati per sesso.
media pesi neonati maschi: 3408.496
media pesi neonate femmine: 3161.061 g
In valore assoluto tale differenza risulta essere pari 247.405grammi a favore dei neonati di sesso maschile. Studiamo se dal punto di vista statistico è possibile affermare che le due medie ottenute nei due sessi differiscono significativamente fra loro:
t.test(Peso~Sesso)##
## Welch Two Sample t-test
##
## data: Peso by 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.496Il test t applicato ai due sessi dei neonati ha fornito che le due medie differiscono significativamente, con una differenza in valore assoluto di 247.405g a favore dei maschi, e un p-valore prossimo allo 0 (< 2.2e-16) che decreta la significatività della differenza fra le medie, quindi rigettando l’ipotesi nulla H0 per la quale le due medie nei sessi non differiscano fra loro. Statisticamente, quindi, il peso dei neonati maschi è maggiore di quello femminile.
1.6 - Variabili di controllo parto (Durata Gestazione e Tipo di Parto) vs Peso Neonato
1.6.1 - Peso del neonato ~ Durata della Gestazione
Varrà sicuramente la pena indagare circa le informazioni che potrebbero restiturici informazioni sulla relazione numerica le settimane di gestazione del feto e il peso del dello stesso neonato alla nascita:
plot(Gestazione, Peso, main = "Gestazione vs Peso", xlab="Settimane di Gestazione", ylab="Peso Neonato (grammi)")
Nei relativi grafici in cui vengono messe in relazioni 2 queste variabili, e si nota che con il passare delle settimane di gestazione, i pesi delle distribuzioni di neonati tendono chiaramente a diventare sempre maggiori con il passare del tempo di gestazione. Vale la pena indagare numericamente le correlazioni lineari fra settimana di gestazione e relativo peso del neonato:
test_pearson <- cor.test(Gestazione, Peso, method = "pearson")
test_pearson##
## Pearson's product-moment correlation
##
## data: Gestazione and Peso
## t = 36.694, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5658780 0.6168568
## sample estimates:
## cor
## 0.5919592Il test R^2 di Pearson conferma una relazione significativa (p-value < 2.2e-16) fra la settimana di Gestazione e il peso del neonato, con un indice di correlazione lineare di Pearson di 0.59, rigettando l’ipotesi nulla H0 di assenza di correlazione. Il test t, conferma che l’associazione fra le due variabili è di tipo lineare fra Gestazione vs Peso e nello specifico positiva, quindi ed è lecito aspettarsi che all’aumentare di una variabile indipendente gestazione possa anche aumentare (in media) la variabile target dipendente come il peso.
1.6.2 Peso del neonato ~ Tipo Parto
Indaghiamo la relazione il peso dei neonati registrato alla nascita e la tipologia di parto che la madre del neonato ha dovuto sostenere:
boxplot(Peso~Tipo.parto,
main="Boxplot Pesi dei neonati \ncondizionati alla tipologia di Parto",
xlab="(Ces) Cesario - (Nat) Naturale",
ylab="Peso neonato (g)")
#analisi numerica
df_parto<-as.data.frame(cbind(Parto_Cesario=mean(Peso[Tipo.parto=="Ces"]),
Parto_Naturale=mean(Peso[Tipo.parto=="Nat"])))
rownames(df_parto) <- c("Media(grammi)")
#visualizzazione
kable_fun(df_parto,
"Tab.15 - Medie dei pesi dei neonati a confronto per tipologia di parto") | Parto_Cesario | Parto_Naturale | |
|---|---|---|
| Media(grammi) | 3282.047 | 3285.063 |
Dal punto di vista numerico (Tab.15) le due medie sono praticamente uguali, ma possiamo ugualmente saggiare le differenze delle medie fra le due categorie di parto:
t.test(Peso~Tipo.parto)##
## Welch Two Sample t-test
##
## data: Peso by Tipo.parto
## t = -0.13626, df = 1494.4, p-value = 0.8916
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
## -46.44246 40.40931
## sample estimates:
## mean in group Ces mean in group Nat
## 3282.047 3285.063- NON si osservano differenze significative fra le medie dei bambini nati da parti Naturali o da parto cesario in base ai dati del dataset fornito in partenza (p-value 89%). Non rigetta l’ipotesi nulla di uguaglianza fra le medie dei pesi dei neonati. In altre parole, la tipologia di parto non si ripercuote direttamente sul Peso del neonato registrato alla nascita.
2. CREAZIONE DEL MODELLO DI REGRESSIONE
2.1 - Matrice di correlazione lineare del dataset
Nei punti precedenti abbiamo risposto ad eventuali quesiti posti dalla consegna, e fatto chiarezza sull’associazione numerica fra variabili indipendenti e di risposta del futuro modello di regressione che andremo a costruire. In questa sezione facciamo un recap generale delle correlazioni di tutte le variabili all’interno del dataset (Tabella 16), forniremo una visione grafica di insieme e poi inizieremo con la costruzione del modelli lineare di regressione multipla a scopi previsionali
dati2 <- dati
dati2$Sesso2 <- as.numeric(factor(dati$Sesso, levels = c("F", "M")))
dati2$Ospedale2 <- as.numeric(factor(dati$Ospedale, levels = c("osp1", "osp2","osp3")))
dati2$Tipo.parto2 <- as.numeric(factor(dati$Tipo.parto, levels=c("Nat","Ces")))
dati2 <- as.data.frame(round(cor(dati2[c("Peso", "Lunghezza", "Cranio", "Anni.madre", "Gestazione","Fumatrici", "N.gravidanze", "Sesso2", "Ospedale2", "Tipo.parto2")]),2))
kable_fun(dati2, "Tab.16 - Matrice delle correlazioni fra le variabili del dataset neonati.csv")| Peso | Lunghezza | Cranio | Anni.madre | Gestazione | Fumatrici | N.gravidanze | Sesso2 | Ospedale2 | Tipo.parto2 | |
|---|---|---|---|---|---|---|---|---|---|---|
| Peso | 1.00 | 0.80 | 0.70 | -0.02 | 0.59 | -0.02 | 0.00 | 0.24 | 0.03 | 0.00 |
| Lunghezza | 0.80 | 1.00 | 0.60 | -0.06 | 0.62 | -0.02 | -0.06 | 0.19 | 0.01 | 0.04 |
| Cranio | 0.70 | 0.60 | 1.00 | 0.02 | 0.46 | -0.01 | 0.04 | 0.15 | 0.01 | 0.00 |
| Anni.madre | -0.02 | -0.06 | 0.02 | 1.00 | -0.13 | 0.01 | 0.38 | 0.01 | 0.03 | 0.00 |
| Gestazione | 0.59 | 0.62 | 0.46 | -0.13 | 1.00 | 0.03 | -0.10 | 0.13 | 0.02 | 0.02 |
| Fumatrici | -0.02 | -0.02 | -0.01 | 0.01 | 0.03 | 1.00 | 0.05 | 0.01 | -0.02 | -0.02 |
| N.gravidanze | 0.00 | -0.06 | 0.04 | 0.38 | -0.10 | 0.05 | 1.00 | 0.02 | 0.01 | 0.02 |
| Sesso2 | 0.24 | 0.19 | 0.15 | 0.01 | 0.13 | 0.01 | 0.02 | 1.00 | 0.01 | 0.00 |
| Ospedale2 | 0.03 | 0.01 | 0.01 | 0.03 | 0.02 | -0.02 | 0.01 | 0.01 | 1.00 | -0.02 |
| Tipo.parto2 | 0.00 | 0.04 | 0.00 | 0.00 | 0.02 | -0.02 | 0.02 | 0.00 | -0.02 | 1.00 |
# Creazione del dataframe
correlazioni <- data.frame(
Variabile = c("Lunghezza", "Cranio", "Gestazione", "Sesso", "Ospedale",
"N.gravidanze", "Tipo.parto", "Anni.madre", "Fumatrici"),
Correlazione = c(0.8, 0.7, 0.59, 0.24, 0.03, 0, -0.02, -0.02, -0.02)
)
colnames(correlazioni) <- c("variabile", "Correlazione con Peso")
# Visualizzazione del dataframe
kable_fun(correlazioni, "Tab.17-Correlazioni Lineari con il target")| variabile | Correlazione con Peso |
|---|---|
| Lunghezza | 0.80 |
| Cranio | 0.70 |
| Gestazione | 0.59 |
| Sesso | 0.24 |
| Ospedale | 0.03 |
| N.gravidanze | 0.00 |
| Tipo.parto | -0.02 |
| Anni.madre | -0.02 |
| Fumatrici | -0.02 |
La tabella 16 rappresenta una matrice di correlazione generale a coppie fra le variabili, in cui ogni ad ogni incrocio di elemento nella tabella, descrive il valore di correlazione lineare tra due variabili. Mentre la tabella 17 mostra una versione sintetica ed ordinata delle correlazioni delle variabili indipendenti con la variabile target Peso.
Osservazioni
Le variabili fisiologiche del neonato (Peso, Lunghezza, Cranio) e le settimane di gestazione risultano correlate positivamente tra loro:
Correlazioni della variabile con il target (Peso Neonato): Correlato positivamente con Lunchezza (0.80), Cranio (0.70) e con Gestazione (0.59). In misura minore Peso correlato con Sesso del neonato (0.24).
Correlazione positiva fra la variabile Gestazione con la Lunghezza (0.62) e con il cranio
Correlazione positiva fra la variabile Gestazione con la Lunghezza (0.62) e con Cranio (0.46)
Anche Lunghezza e Cranio risultano essere correlate fra loro (0.60)
correlazione nulla fra la variabile Peso con ed altri fattori come Anni madre (-0.02), Fumatrici (-0.02) e N.gravidanze (0.00), indicando che queste variabili indipendenti abbiano un’influenza molto marginale sulle variabile Peso.
Anni madre e N.gravidanze mostrano una certa relazione, suggerendo un legame demografico o biologico, probabilmente legato al fatto che madri più anziani hanno avuto più opportunità di avere molteplici gravidanze nel trascorso della propria vita.
2.2 - Approfondimento grafico dei plot bi-dimensionali
Di seguito una sintetica disamina grafica d’insieme
panel.smooth.custom <- function(x, y, col = "red",lwd = 2,cex = 2,...)
{points(x, y, ...)
lines(lowess(x, y),
col = col,
lwd = lwd, ...)}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{par(usr = c(0, 1, 0, 1))
r <- (cor(x, y, use = "complete.obs"))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if (missing(cex.cor)) cex.cor <- 0.8 / strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)}
# Creazione della matrice di scatterplot
pairs(dati[c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")],
upper.panel = panel.smooth.custom,
lower.panel = panel.cor,
main="Correlazioni lineari (parte inferiore) - Scatterplot (parte superiore)")
In aggiunte alle considerazioni precedenti notiamo che la correlazione della variabile Lunghezza e Peso abbia forma quadratica, come se la lunghezza aumentasse in forma quadratica in rapporto al peso
Il rapporto fra Cranio e Lunghezza aumenta in modo lineare durante le prime fasi della gestazione per poi appiattirsi durante le ultime settimane della gestazione che portano al parto della neonato. Questa correlazione fra variabili indipendenti Cranio e Lunghezza dovrà essere sicuramente attenzionata
2.3 La costruzione del modello multivariato
Passiamo quindi alla costruzione del migliore modello lineare multivariato, a partire dai dati a disposizione. Il modello di regressione lineare multivariato è una modellazione matematica della realtà a più fattori, che cerca di approssimare (generalizzare) un fenomeno reale al fine di descrivere matematicamente e/o prevedere fenomeni supposti dipendenti a partire dalle informazioni a disposizione.
Dal punto di vista matematico è utilizzato per descrivere la relazione tra una variabile dipendente (o di risposta) Y e più variabili indipendenti X1,X2,…,Xn. L’obiettivo è prevedere Y o comprenderne le variazioni in base alle X, assumendo una relazione lineare del tipo: Y=β0+β1X1+β2X2+⋯+βnXn+ ε, dove i vari termini rappresentano:
- β0 è un termine costante, detto intercetta, (valore di Y quando tutte le Xi sono zero), e rappresenta il valore medio (o atteso) della variabile dipendente Y, quando tutte le altre variabili indipendenti Xi del modello sono pari a 0;
- β1,β2,…,βn sono i coefficienti che quantificano l’influenza delle rispettive X, vengono detti “regressori del modello”;
- ε, rappresenta la parte erratica del modello, ovvero l’errore residuo sulla variabile dipendente Y che il modello non riesce a generalizzare in modo esatto (con esatta accuratezza), e quindi quella parte di informazione che esso non riesce a spiegare a causa della variabilità statistica, o per la scelta di avere modelli più o meno robusti (Bias-Variance Trade-off). La diagnostica dei residui erratici del modello scelto, è un operazione sistematica.
Calato nel contesto di questo studio, potremmo utilizzare come regressori del modello lineare variabili supposte indipendenti e significative nel modello, Lunghezza, Cranio e Gestazione, Fumatrici, Sesso al fine di prevedere e descrivere il peso teorico dei neonati, e quindi fornire informazioni predittive e pro-attive alla struttura ospedaliera al fine di contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere.
2.4 MODELLO 1 (mod1)
Per puro scopo dimostrativo, iniziamo con il costruire il primo modello lineare che comprende tutte le variabili del dataset, comprendendo all’interno della sua descrizione anche variabili non significative o correlate linearmente con la variabile target ai fini della sua predizione.
Il primo passaggio quindi sarà quello di creare un modello generale con tutte le variabili, spiegare i risultati per poi passare ai modelli successivi eliminando le variabili meno influenti:
mod1<- lm(Peso~.,data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-162.4.1 - DESCRIZIONE SIGNIFICATO GENERICO DELLE VOCI NEL MODELLO
Dal sommario appena illustrato, notiamo le performance del primo modello di regressione lineare multipla, costruito con tutte le variabili del dataset. Descriviamolo per intero la prima volta in modo da enunciare le voci in esso contenuto:
“Call:” lm(formula = Peso ~ ., data = dati): questa è la chiamata al modello. La formula Peso ~ . indica che stai si sta predicendo il target Peso utilizzando tutte le altre variabili presenti nel dataset dati come predittori (indicate con .).
I “residuals” Questi rappresentano le differenze tra i valori osservati e quelli predetti dal modello per la variabile dipendente Peso. Vengono individuati il range e gli indici di posizione della distribuzione dei residui, ovvero la parte erratica del modello che non riesce a generalizzare nel modello.
Il “Coefficients Estimate” è il coefficiente stimato per ciascuna variabile. Indica quanto cambia, in media, la variabile dipendente (Peso) quando la variabile predittiva aumenta di una unità, tenendo le altre variabili costanti.
Lo “Std. Error (Errore standard)“misura l’incertezza nella stima del coefficiente. Coefficienti con un errore standard basso sono stimati con maggiore precisione.
Il “t value” è il rapporto tra il coefficiente stimato e il suo errore standard. Un valore t alto indica che il coefficiente è significativamente diverso da zero. Il”Pr(>|t|)“ è il valore p associato al test t. Indica la probabilità che il coefficiente sia uguale a zero sotto l’ipotesi nulla. Per una lettura rapida del sommario, il modello lineare riporta i codici di significatività sotto forma di asterischi:
- ***molto significativo,**significativo, *moderatamente significativo, . marginalmente significativo, “valore vuoto” non significativo
“Residual standard error” è la dev.standard dei residui, che misura quanto i valori osservati si discostano da quelli generalizzati dal modello. Un errore standard più basso indica un modello più accurtato. In questo caso, è 274 con 2487 gradi di libertà. Il “Multiple R-squared” indica la proporzione di varianza nella variabile dipendente spiegata dal modello. Valori vicini a 1 indicano un buon adattamento al modello. L’”Adjusted R-squared” è una versione corretta dell’ R2 che tiene conto del numero dei predittori nel modello. Serve per evitare sovrastime di R^2, legate all’aggiunta di varianza che potrebbero apportare l’ingresso nel modello di variabili non significative. Per concludere l’ F-statistic è un test globale che valuta se almeno una delle variabili predittive ha un effetto significativo sulla variabile dipendente. Il suo valore accoppiato, con un valore di p-value molto basso, indica che il modello ha capacità predittive sulla variabile dipendente in modo significativo.
INTEPRETAZIONE DEI COEFFICIENTI DEL MODELLO 1 (mod1):
Intercetta: pari a -6735.7960g, rappresenta il valore numerico del peso previsto del neonato quando tutte le variabili indipendenti assumono valore nullo. Questo termine ha un valore poco pratico, dato che alquanto improbabile registrare un peso di un neonato negativo, o che tutte le variabili indipendenti assumano valore nullo;
N.gravidanze: Il coefficiente è 11,38 (p-value di 0.015), suggerendo che per ogni gravidanza aggiuntiva, il peso del neonato aumenta di circa 11.4 grammi. All’interno del contesto di questo modello generale è leggermente significativa, ma come vedremo si andrà eliminare procedendo verso modelli più efficienti;
Gestazione: Il coefficiente di 32.58 (p-value < 2e-16) suggerisce che per ogni settimana in più di gestazione, il peso del neonato aumenta in media di 32.6 grammi. Questo variabile è molto significativa in quanto il p-value è pressoché nullo;
Lunghezza: per ogni millimetro di aumento della lunghezza del neonato si registra in media un aumento di un 10,3 grammi nel peso del neonato (p-value < 2e-16). Questo variabile è molto significativa in quanto il p-value è pressoché nullo;
Cranio: per ogni millimetro di aumento del diametro craniale del neonato, si registra in media un aumento di un 10,5 grammi nel peso del neonato (p-value < 2e-16); Questo variabile è molto significativa in quanto il p-value è pressoché nullo;
Tipo.parto. Il peso dei neonati nati da parto NATURALE è in media 29,63 grammi maggiore rispetto ai neonati nati da parto Cesario (p-value 0.014). Valore di p-value al di sotto della soglia di significatività alpha del 5%;
Ospedale. Il peso dei neonati nati in Ospedale 3 è in media 28.2495 grammi superiore rispetto a quelli nati in Ospedale 1 (p =0.0366=0.0366) mentre l’Ospedale 2 non mostra differenze significative rispetto all’Ospedale 1.
Sesso. Il coefficiente di 77.5723 indica che, in media, i neonati maschi pesano circa 77.6 grammi in più rispetto alle femmine (p <0.001<0.001).
R^2 aggiustato del modello è abbastanza alta (0.7278), con un pvalue della statistica F prossimo a 0.
Non ha senso tenere la variabile tipo Ospedale 1,2,3 all’interno del dataset in quanto poco significative, e probabilmente gli effetti di significatività legati al target Peso sono da ricercarsi nella maggiore presenza di bambini maschi nell’ospedale 3, rispetto alla baseline (osp1). Quindi dalle prossime previsione iniziamo con escludere queste alcune variabili non significative.
2.5 MODELLO 2 (eliminazione Anni.madre e Ospedale dal modello)
Dal precedente modello 1 verranno eliminate 2 variabili non significative, e cioè la variabile Anni.madre e quella dell’ Ospedale di riferimento. Mentre per la prima la ragion è presto detta, a causa della bassa significatività predittiva all’interno del modello 1 (p.value = 0.485), per la seconda si deve fare un discorso leggermente diverso.
Chiaramente non si può associare una dipendenza causa-effetto fra il peso del neonato e l’azienda ospedaliera in cui il neonato è stato partorito, quindi questo basterebbe (logicamente) per non includere tale variabile all’interno del modello, però osserviamo un comportamento alquanto insolito all’interno del summary della primo modello. Infatti , nel mod1 si può leggere che esiste differenza significativa (positiva) di peso, rispetto alla baseline (osp1), con i pesi registrati dei neonati nati in ospedale 2 e 3, la quale è in media di 28.25g maggiore rispetto alle natalità generale registrata nella baseline del 1° ospedale. Una possibile spiegazione di questo “controsenso”, potrebbe ricercarsi in un fattore puramente numerico:
avendo constatato che la media dei pesi dei neonati, raggruppati per sesso, è significativamente differente, e cioè che la media dei pesi dei maschi >media dei pesi delle femmine, probabilmente il fatto che i neonati nell’ospedale 3 abbiano una significatività maggiore rispetto alla baseline ospedale 1, può è legato ad un “effetto matematico”, facilmente dimostrabile con una tabella di contingenza Ospedale~Sesso (trasformata in frequenze relative).
boxplot(Peso~Ospedale,
main="Boxplot Pesi neonati \ncondizionato per Struttura Ospedaliera",
y="Peso neonato (g)")
abline(h=mean(Peso[Ospedale=="osp1"]),col="red")
abline(h=mean(Peso[Ospedale=="osp2"]),col="blue")
abline(h=mean(Peso[Ospedale=="osp3"]), col="green")
cont_Sex_Hosp<- as.data.frame(dati)%>%
group_by(Ospedale,Sesso)%>%
count()
cont_Sex_Hosp_wide <-pivot_wider(
data=cont_Sex_Hosp,
names_from = Sesso,
values_from = n)
kable_fun(cont_Sex_Hosp_wide, "Tab.18 - Contingenza Sesso per unità ospedaliera")| Ospedale | F | M |
|---|---|---|
| osp1 | 409 | 407 |
| osp2 | 434 | 414 |
| osp3 | 412 | 422 |
Nella tabella 18 di contingenza, sintetizzata per sesso e struttura ospedaliera, è possibile osservare il numero di neonati maschi nati nelle ospedali 2 e 3 sono maggiori rispetto a quelli nati nell’ospedale 1 (base line). Qui la “variabile nascosta” è il sesso del neonato che, come dimostrato in precedenza, fa variare intrinsecamente il peso nei maschi, il quale è significativamente in media maggiore rispetto alla media dei pesi di una neonata. E’ probabile, quindi, che questo effetto combinato sesso nascituro e numero maggiori di neonati maschi nati nelle aziende ospedaliere 2 e 3, porti al prossimo modello di regressione lineare a considerare significativa la differenza di peso dei neonati fra le aziende ospedaliere. Giustificanto tale significatività, per motivi puramente di calcolo del modello, oltre alla logica disconnessione causa-effetto fra le varibili, preferisco eliminare anche il vettore Ospedale dal modello 1, aggiornandolo alla versione 2 del modello, ovvero mod2. Si osserva che l’eliminazione
mod2 <- update(mod1,~.-Anni.madre-Ospedale-Fumatrici)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.14 -181.97 -16.26 160.95 2638.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.0171 136.0715 -49.298 < 2e-16 ***
## N.gravidanze 12.7356 4.3385 2.935 0.00336 **
## Gestazione 32.3253 3.7969 8.514 < 2e-16 ***
## Lunghezza 10.2833 0.3009 34.177 < 2e-16 ***
## Cranio 10.5063 0.4263 24.648 < 2e-16 ***
## Tipo.partoNat 30.1601 12.1027 2.492 0.01277 *
## SessoM 77.9171 11.1994 6.957 4.42e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16Dal summary del modello 2 è possibile notare che l’eliminare delle 2 variabili non influenti sul modello abbia comunque mantenuto le medesime capacità del modello 2 di predire la variabile target dello studio, infatti passando dal modello 1 al modello 2, l’R2 adj. è rimasto pressoché identico (da 0.7278 –> 0.727).
2.6 MODELLO 3 (eliminazione variabile Tipo.parto dal modello)
In accordo con i precedenti studi (punto 1.6.2, test t delle medie fra tipologie di parti) non è stato possibile dimostrare una differenza significativa fra le medie dei pesi dei neonati in funzione del tipo di parto, e questo per certi versi potrebbe avere una logica di disconnesione causa effetto fra le variabili. Quindi anche se il modello interpreta una correlazione lineare, anche se non troppo marcata, decido di eliminarla dal ultimo modello a disposizione, per crearne uno nuovo mantenendo tutte le variabili precedenti. Se l’R2 aggiustato del nuovo modello mod3 si mantiene pressoché costante manterrò le modifiche:
mod3 <- update(mod2,~.-Tipo.parto)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16L’R2 adj. è rimasto pressoché lo stesso (da 0.727 -> 0.7265), quindi abbiamo ottenuto un modello più leggero che ha la stessa capacità di spiegare la variabilità dei dati tramite un R2 adj. simile, dove al suo interno vi sono solo variabili indipendenti significative. E’ lecito pensare adesso di aver ottenuto un potenziale modello multivariato. Proviamo anche a testare la sua stabilità tramite analisi VIF delle variabili indipendenti che lo compongono. ll VIF (Variance Inflation Factor) è un indice statistico utilizzato per rilevare la multicollinearità tra le variabili indipendenti in un modello di regressione. La multicollinearità si verifica quando due o più variabili predittive sono altamente correlate tra loro (es. Cranio e Lunghezza come abbiamo visto dai grafici precedenti), il che può influire negativamente sull’interpretazione e sulle stime del modello. Per avere un modello robusto e stabile, il VIF di ogni variabile predittiva non deve superare il valore di 5
kable_fun(round(vif(mod3),2), "Tab.19 - VIF DEL MODELLO 3")| x | |
|---|---|
| N.gravidanze | 1.02 |
| Gestazione | 1.67 |
| Lunghezza | 2.08 |
| Cranio | 1.62 |
| Sesso | 1.04 |
Il terzo modello (mod3) si dimostra stabile ed efficace nel descrivere la variabile target. Nella prossima sezione del progetto, prima di passare all’analisi diagnostica degli errori del modello mod3, effettueremo una verifica parallela utilizzando procedure automatiche di stepwise selection per identificare le variabili significative attraverso tre ulteriori test.
2.6 Procedura di selezione delle variabili
Analisi ANOVA (Analysis of variance) basata sulla statistica del test F
Akaike Information Criterion (AIC)
Bayesian Information Criterior (BIC)
A prescindere dal metoto utilizzato per la selezione automatica delle variabili è bene partire da un modello di regressione lineare completo, che include tutte le variabili disponibili (significative e non) in aggiunta all’intercetta. A ogni passaggio, si decide se rimuovere una variabile (o eventualmente reinserirla dal secondo passaggio in poi) valutando il valore della statistica F, dell’AIC o del BIC. Il processo continua fino a quando non è più possibile migliorare il modello in base al criterio scelto.
2.6.1 Anova - Statistica F
Un metodo per capire se l’aggiunta di una variabile apporta modifiche significative al modello può essere quello dell’ANOVA, che basandosi sulla differenza di varianze di 2 modelli lineari multivariati a confronto, esegue un TEST F (basato sulla distribuzione Fisher-Snedecor) e saggia l’ipotesi di differenza significatività nell’informazione delle varianze a seguito della comparazione di modelli:
- L’ipotesi nulla H0: non c’è significatività per la differenza fra i due modelli, per p-value > 5%
- L’ipotesi alternativa H1: c’è significatività per la differenza fra i due modelli, per p-value < 5%
anova(mod3,mod2,mod1)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
## Model 3: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2491 187574428 1 467626 6.2277 0.01264 *
## 3 2487 186743194 4 831234 2.7675 0.02601 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L’analisi anova indica che nel processo manuale di selezione delle variabili fatte in precedenza, creando i modelli 1,2 e 3, fra le variabili eliminate (Tipo.parto, Ospedale , fumatrici, anni madre) vi è stata una differenza significativa fra i due modelli solo quando non si prende più in considerazione Tipo.parto e Ospedale, ma dato che queste due variabili potrebbero essere escluse logicamente a priori, continuiamo a prendere per buona la selezione manuale fatta in precedenza e verifichiamo cosa suggeriscono in parallelo a questo studio, i criteri di informazione AIC e BIC.
2.6.2 AIC (Akaike Information Criterion)
L’AIC (Akaike Information Criterion) è un criterio utile per selezionare un modello, tenendo conto sia della capacità del modello di adattarsi ai dati (bontà di adattamento) sia della sua complessità. In particolare, l’AIC applica una penalizzazione costante, pari a 2 per ogni parametro aggiuntivo del modello (k=2), per evitare di favorire modelli troppo complessi.
Tra i modelli confrontati, viene scelto quello con il valore di AIC più basso, perché rappresenta il miglior compromesso tra semplicità e accuratezza nell’adattamento ai dati.
stepwise_AIC <- step(mod1, direction = "both", k=2)## Start: AIC=28054.55
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36710 186779904 28053
## - Fumatrici 1 90677 186833870 28054
## <none> 186743194 28055
## - N.gravidanze 1 446244 187189438 28059
## - Tipo.parto 1 451073 187194266 28059
## - Ospedale 2 687555 187430749 28060
## - Sesso 1 3610705 190353899 28100
## - Gestazione 1 5458852 192202046 28125
## - Cranio 1 45318506 232061700 28595
## - Lunghezza 1 87861708 274604902 29016
##
## Step: AIC=28053.05
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 91599 186871503 28052
## <none> 186779904 28053
## + Anni.madre 1 36710 186743194 28055
## - Tipo.parto 1 452049 187231953 28057
## - Ospedale 2 693914 187473818 28058
## - N.gravidanze 1 631082 187410986 28060
## - Sesso 1 3617809 190397713 28099
## - Gestazione 1 5424800 192204704 28123
## - Cranio 1 45569477 232349381 28596
## - Lunghezza 1 87852027 274631931 29014
##
## Step: AIC=28052.27
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 186871503 28052
## + Fumatrici 1 91599 186779904 28053
## + Anni.madre 1 37633 186833870 28054
## - Tipo.parto 1 444404 187315907 28056
## - Ospedale 2 702925 187574428 28058
## - N.gravidanze 1 608136 187479640 28058
## - Sesso 1 3601860 190473363 28098
## - Gestazione 1 5358199 192229702 28121
## - Cranio 1 45613331 232484834 28596
## - Lunghezza 1 88259386 275130889 29017la procedura ha portato all’eliminazione delle variabili Anni.madre e Fumatrici, ritenute dal modello non significative dal punto di vista della loro capacità di migliorare la spiegazione delle variabilità nella variabile target (Peso). Più in particolare il criterio informativo RSS (residual sum of squares) è stato minimizzato da 28054.55 -> 28052.27, ed interrotto in quanto la riaggiunta di una delle 2 variabili già escluse, o l’eliminazione di altre nel modello, porterebbe ad un nuovo incremento del parametro AIC non voluto.
2.6.3 BIC (Bayesian Information Criterion)
Il BIC (Bayesian Information Criterion) è un’altra misura utilizzata nella statistica per valutare e confrontare modelli di regressione o modelli probabilistici. È un criterio di selezione del modello che tiene conto della complessità del modello e della bontà di adattamento ai dati. Come per l’AIC, anche il per il criterio di informazione Bayesiano viene preferito il modello con il valore di BIC più basso:
stepwise_BIC <-step(mod1, direction = "both", k = log(nrow(dati)))## Start: AIC=28118.61
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36710 186779904 28111
## - Fumatrici 1 90677 186833870 28112
## - Ospedale 2 687555 187430749 28112
## - N.gravidanze 1 446244 187189438 28117
## - Tipo.parto 1 451073 187194266 28117
## <none> 186743194 28119
## - Sesso 1 3610705 190353899 28159
## - Gestazione 1 5458852 192202046 28183
## - Cranio 1 45318506 232061700 28654
## - Lunghezza 1 87861708 274604902 29074
##
## Step: AIC=28111.28
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 91599 186871503 28105
## - Ospedale 2 693914 187473818 28105
## - Tipo.parto 1 452049 187231953 28110
## <none> 186779904 28111
## - N.gravidanze 1 631082 187410986 28112
## + Anni.madre 1 36710 186743194 28119
## - Sesso 1 3617809 190397713 28151
## - Gestazione 1 5424800 192204704 28175
## - Cranio 1 45569477 232349381 28649
## - Lunghezza 1 87852027 274631931 29066
##
## Step: AIC=28104.68
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 702925 187574428 28098
## - Tipo.parto 1 444404 187315907 28103
## <none> 186871503 28105
## - N.gravidanze 1 608136 187479640 28105
## + Fumatrici 1 91599 186779904 28111
## + Anni.madre 1 37633 186833870 28112
## - Sesso 1 3601860 190473363 28145
## - Gestazione 1 5358199 192229702 28168
## - Cranio 1 45613331 232484834 28642
## - Lunghezza 1 88259386 275130889 29063
##
## Step: AIC=28098.41
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 467626 188042054 28097
## <none> 187574428 28098
## - N.gravidanze 1 648873 188223301 28099
## + Ospedale 2 702925 186871503 28105
## + Fumatrici 1 100610 187473818 28105
## + Anni.madre 1 44184 187530244 28106
## - Sesso 1 3644818 191219246 28139
## - Gestazione 1 5457887 193032315 28162
## - Cranio 1 45747094 233321522 28636
## - Lunghezza 1 87955701 275530129 29051
##
## Step: AIC=28096.81
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188042054 28097
## - N.gravidanze 1 621053 188663107 28097
## + Tipo.parto 1 467626 187574428 28098
## + Ospedale 2 726146 187315907 28103
## + Fumatrici 1 92548 187949505 28103
## + Anni.madre 1 45366 187996688 28104
## - Sesso 1 3650790 191692844 28137
## - Gestazione 1 5477493 193519547 28161
## - Cranio 1 46098547 234140601 28637
## - Lunghezza 1 87532691 275574744 29044summary(stepwise_BIC)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16Il criterio BIC conferma che il migliore compromesso è stato raggiunto dal mod3 con il valore di BIC=35207.48, il più basso fra i modelli creati manualmente (mod1, 2 e 3). Il criterio BIC, infatti , tende a premiare l’efficienza di modelli semplici e non troppo parametrizzati (a differenza del AIC che tende a preferire gli iperparametrizzati)
2.6.4 Considerazioni intermedie sul modello 3 (mod3)
Ricapitolando, dal punto di vista della solidità previsionale del modello 3, abbiamo:
\(R^2\) \(= 0.7264, con\) \(p-value <2.2e-16\)
VIF di tutte le variabili significative che lo compongono <5
Il modello mod3 può passare quindi alla fase di diagnostica dei residui, e fatte le doverose premesse, possibilmente si potranno iniziare a fare le prime previsioni se il modello 3 supererà la diagnostica della parte erratica.
Volendo, però, possiamo costruire anche altri modelli che tengono conto delle relazioni quadratiche con il target Peso, o altri in cui si tiene conto dell’influenza reciproca delle variabili indipendenti
2.7 MODELLO 4 (interazione quadratica fra Lunghezza e il Peso)
Ricordando la linea di tendenza (modello generalizzato) riportata nel
grafico del plot bi-dimensionale dei valori assunti della variabile peso
neonato in funzione della lunghezza dello stesso (sezione 2.2), la forma
di tale interpolazione ricorda quella di una funzione quadratica. Per
avvalorare tale intuizione dal punto di vista statistico numerico,
procedo con inserire una relazione quadratica nel modello lineare 3, e
valutare i risultati battezzando il nuovo modello che include la
modifica mod4:
mod4 <- update(mod3,~.+I(Lunghezza^2))
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1169.62 -181.77 -12.79 163.77 1786.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 212.288548 723.852095 0.293 0.769336
## N.gravidanze 14.085464 4.266175 3.302 0.000975 ***
## Gestazione 42.551398 3.876629 10.976 < 2e-16 ***
## Lunghezza -20.267001 3.162718 -6.408 1.76e-10 ***
## Cranio 10.651783 0.418894 25.428 < 2e-16 ***
## SessoM 69.968733 11.038797 6.338 2.75e-10 ***
## I(Lunghezza^2) 0.031655 0.003267 9.690 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7369, Adjusted R-squared: 0.7363
## F-statistic: 1163 on 6 and 2491 DF, p-value: < 2.2e-16kable_fun(round(vif(mod4),2), "Tab.20 - VIF DEL MODELLO 4")| x | |
|---|---|
| N.gravidanze | 1.03 |
| Gestazione | 1.80 |
| Lunghezza | 238.01 |
| Cranio | 1.63 |
| Sesso | 1.05 |
| I(Lunghezza^2) | 230.03 |
Interpretazione dei coefficienti mod4:
(Intercept): 212.29, ma non significativo, quindi non rilevante;
N.gravidanze: Ogni aumento di 1 nel numero di gravidanze aumenta la variabile dipendente in media di 14.09 grammi;
Gestazione: Ogni aumento di 1 settimana di gestazione, in media si incrementa la variabile dipendente di 42.55 grammi;
Lunghezza: Ogni aumento di 1 millimetro nella lunghezza del neonato si riduce in media la variabile dipendente di 20.27 grammi;
Cranio: Ogni aumento di 1 millimetro nella misura del cranio del neonato aumenta la variabile dipendente di 10.65 grammi.
SessoM: Il sesso maschile (rispetto al sesso di riferimento, femminile) aumenta in media la variabile dipendente di 69.97 grammi.
I(Lunghezza^2): L’effetto quadratico della lunghezza aggiunge 0.032 grammi per ogni incremento al quadrato della lunghezza, suggerendo una relazione non lineare.
Con l’aggiunta del termine quadratico \(Lunghezza^2\) nel modello, si osservano 3 fenomeni principali:
- L’aggiunta del termine quadratico nel modello riesce a migliorare l’R^2 aggiustato del modello rispetto al mod3 (0.7265 -> 0.7363) con un’ottima significatività della statistica F (p-value < 2.2e-16),con tutte le variabili esplicative che restano significative per il modello
- L’intercetta ha perso di significatività con un p-value elevato (0.769) ed un coefficiente di +212.29, facendo intendere che il modello non riesce a generalizzare correttamente quando tutte le variabili indipendenti tendono a zero. Questo fatto possiamo trascurarlo, in quanto nella pratica è difficile che ci possiamo ritrovare in una situazione del genere
- Il VIF del mod4 è variata rispetto al precedente modello più semplice, infatti ora le variabili Lunghezza e \(Lunghezza^2\) hanno un valore di VIF >>5 ( circa di 238 ognuno) e questo comporta instabilità al modello.
Per scopi di ricerca manteniamo il modello 4, un modello separato rispetto al modello 3 (preferito e più semplice), e continuiamo con con l’espansione alternativa della complessità del modello 3, provando effetti di interazioni fra le variabili Lunghezza e Cranio del neonato con nuovo modello 5 a seguire.
2.8 MODELLO 5 (modello con interazione fra le variabili Lunghezza e Cranio)
Come osservato in precedenza le 2 variabili Lunghezza neonato e Diametro del cranio possiedono fra loro un coefficiente di correlazione lineare significativa “r” pari a 0.60. Nel modello multivariato mod5 di questa sezione proviamo ad inserire tale effetto di interazione fra le 2 variabili (partendo dal modello 3) e capire come questa interazione influisce sulla varianza spiegata del modello e la relativa stabilità:
mod5 <- update(mod3,~.+Lunghezza*Cranio)
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Lunghezza:Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.65 -180.93 -13.48 165.99 2865.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.803e+03 1.018e+03 -1.771 0.0767 .
## N.gravidanze 1.293e+01 4.323e+00 2.991 0.0028 **
## Gestazione 3.815e+01 3.967e+00 9.616 < 2e-16 ***
## Lunghezza -3.060e-01 2.203e+00 -0.139 0.8895
## Cranio -4.755e+00 3.192e+00 -1.490 0.1365
## SessoM 7.324e+01 1.120e+01 6.537 7.59e-11 ***
## Lunghezza:Cranio 3.157e-02 6.531e-03 4.835 1.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7296, Adjusted R-squared: 0.7289
## F-statistic: 1120 on 6 and 2491 DF, p-value: < 2.2e-16kable_fun(round(vif(mod5),2), "Tab. 21 - VIF MODELLO 5")## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif| x | |
|---|---|
| N.gravidanze | 1.02 |
| Gestazione | 1.84 |
| Lunghezza | 112.32 |
| Cranio | 91.83 |
| Sesso | 1.05 |
| Lunghezza:Cranio | 313.72 |
- L’interazione fra lunghezza e Cranio migliora in misura marginale l’R2 del modello 3 portandolo nel modello 5 a 0.729, con un’ottima significatività del modello in generale
- Si registra una perdita di significatività nell’intercetta del modello, nella variabile singole di lunghezza e cranio, dove la loro interazione esplicativa nel modello sia risultata più importante (p-value di 1.41e-06) delle singole variabili Lunghezza e Cranio (p-value rispettivamente di 0.89 e 0.14)
- Il VIF delle variabili Lunghezza, Cranio e quella combinata Lunghezza*Cranio superano abbondamente il valore di 5, rispettivamente 112.32, 91.83 e 313.72 , causando instabilità al modello legato a fenomeni di multicollinearità fra le variabili nel modello
Da questo studio otteniamo che esiste una correlazione abbastanza significativa fra le variabili Lunghezza e Cranio che tale combinazione non può essere trascurata. E’ vero anche che l’introduzione del prodotto fra variabili nel modello altera molto la stabilità nel modello 5, preferendo ancora una volta il modello 3 (senza interazioni, ne effetti quadratici) più affidabile e robusto.
3. DIAGNOSTICA DEI MODELLI
La diagnostica dei residui in un modello lineare multivariato è fondamentale per valutare la validità e l’affidabilità del modello. Analizzando i residui, cioè le differenze tra i valori osservati e quelli previsti, è possibile verificare se il modello rispetta le ipotesi statistiche di base. Questo processo permette di controllare se vengono rispettate le assunzioni di linearità del modello:
- se i residui seguono una distribuzione normale (test di Shapiro-Wilk)
- se i residui hanno una varianza costante o omoschedasticità (test di Breush-Pagan)
- non sono affetti da autocorrelazione (test di Durbin-Watson)
- analisi della presenza di outliers nello spazio dei regressori (chiamati leverages o valori di leva) e nello spazio delle variabili risposte (detti outliers), e la loro gestione consapevole. Questi criteri alle volte vengono unificati con il calcolo della distanza di Cook.
3.1 Analisi diagnostica grafica dei modelli
3.1.1 Analisi grafiche e punti influenti
Come detto in precedenza il modello che ho preferito per le previsioni è il modello 3 che risulta essere quello più semplice, robusto ed affidabile. Al fine di promuovere o meno tale modello alla fase finale di previsione del peso dei neonati, dobbiamo prima analizzare la sua parte erratica con le regole generali della diagnostica:
par(mfrow=c(2,2))
plot(mod3)
- Le immagini sopra mostrate riguardano i quattro grafici diagnostici comunemente utilizzati per valutare la bontà di un modello di regressione lineare. Il primo grafico, Residuals vs Fitted, rappresenta i residui sull’asse verticale rispetto ai valori predetti dal modello sull’asse orizzontale. Questo grafico è utile per identificare eventuali schemi sistematici nei residui: idealmente, i punti dovrebbero distribuirsi in modo casuale attorno alla linea orizzontale (a zero). Nel caso mostrato, alcuni punti lontani dalla massa principale potrebbero indicare anomalie o una relazione non lineare tra le variabili.
- Il secondo grafico, il Q-Q Plot, confronta la distribuzione dei residui standardizzati con quella teorica di una normale. Se i residui seguono una distribuzione normale, i punti si dispongono lungo la linea diagonale. Tuttavia, nell’immagine si notano delle deviazioni significative per valori estremi, suggerendo che i residui non seguono perfettamente una distribuzione normale.
- Il terzo grafico, noto come Scale-Location, mostra la radice quadrata del valore assoluto dei residui standardizzati sull’asse verticale rispetto ai valori predetti sull’asse orizzontale. Questo grafico serve a verificare la presenza di eteroschedasticità, ossia variazioni non costanti nella varianza dei residui. Un andamento uniforme dei punti indicherebbe omoschedasticità, ma la presenza di schemi o variazioni potrebbe indicare problemi con la varianza.
- Infine, il quarto grafico, Residuals vs Leverage, evidenzia le osservazioni con alta leva e potenziale influenza sul modello. Qui si valuta l’effetto combinato del residuo standardizzato e della distanza di Cook per identificare osservazioni influenti. Nel grafico, viene identificato un punto influente con un’osservazione di Cook compresa fra 0.5 e 1, e tale punto è relativa all’osservazione “1551”.
Nel complesso, questi grafici suggeriscono la presenza di osservazioni potenzialmente problematiche (come outliers o punti con alta leva) e possibili problemi legati all’eteroschedasticità o alla non normalità dei residui.
- GRAFICI DEI VALORI LEVERAGES/OUTLIERS/DISTANZE DI COOK
Andiamo ad identificare i punti influenti sul modello lineare appena costruito, ed analizzarli caso per caso. Successivamente verranno eseguiti anche i diagnostici di tipo statistico per accompagnare l’analisi grafica ad una di tipo numerica:
print("Analisi dei punti influenti nel modello di regressione lineare multivariata mod3")## [1] "Analisi dei punti influenti nel modello di regressione lineare multivariata mod3"outlierTest(mod3)## rstudent unadjusted p-value Bonferroni p
## 1551 10.046230 2.6345e-23 6.5810e-20
## 155 5.025345 5.3818e-07 1.3444e-03
## 1306 4.824963 1.4848e-06 3.7092e-03#leverage
lev <-hatvalues(mod3)
plot(lev,
main="Leverages\n (Spazio dei regressori)",
ylab="Distanza di Cook"
)
p=sum(lev)
soglia = 2*p/n
#lev[lev>soglia]
abline(h=soglia,col="red")
#outliers
plot(rstudent(mod3),
main="Outliers\n (spazio delle risposte)"
)
abline(h=c(-2,2),col="red")
#Cook's distance
cook <- cooks.distance(mod3)
plot(cook,
main="Verifica punti influenti\nCook's distance (scatterplot)",
ylab="Cook's distance"
)
abline(h=0.5,col="red")
#Cook's distance 2
plot(mod3, which = 4, main = "Verifica Punti influenti")
abline(h=0.5,col="red")
#max di cook
max_distance <- as.data.frame(round(max(cook),2))
colnames(max_distance) <- c("Distanza di Cook")
kable_fun(max_distance, "Tab.21 - Osservazione influente\nn.1551")| Distanza di Cook |
|---|
| 0.83 |
##### osservazione 1551(punto influente)
kable_fun((dati%>%
filter(rownames(dati)=="1551")),
"Tabella 22 - Osservazione influente n.1551 del dataset neonati")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 35 | 1 | 0 | 38 | 4370 | 315 | 374 | Nat | osp3 | F |
#ricerca misure antropometriche
kable_fun(round(dati%>%
filter(Sesso=="F")%>%
filter(Gestazione==38)%>%
summarise(mediaLun=mean(Lunghezza),
sdLunghezza=sd(Lunghezza),
mediaCranio=mean(Cranio),
sdCranio=sd(Cranio)),2),
"Tab.23 - Media e Dev.Standard di misure antropometriche neonate nate in 38 settimane di gestazione (millimetri)")| mediaLun | sdLunghezza | mediaCranio | sdCranio |
|---|---|---|---|
| 485.47 | 22.67 | 337.84 | 12.65 |
#ricerca peso
kable_fun(dati%>%
filter(Sesso=="F")%>%
filter(Gestazione==38)%>%
summarise(MediaNeonate=mean(Peso),
sdNeonate=sd(Peso)),
"Tab.24 - Media e Dev.STD del peso di Neonate (in 38 settimane di Gestazione)")| MediaNeonate | sdNeonate |
|---|---|
| 3125.894 | 390.2284 |
##### osservazione 155 (outlier)
kable_fun((osservazione155 <-dati%>%
filter(rownames(dati)=="155")),
"Tabella 25 - osservazione outlier 155")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 30 | 0 | 0 | 36 | 3610 | 410 | 330 | Nat | osp1 | M |
#ricerca misura Peso (settimana 36)
kable_fun(dati%>%
filter(Gestazione==36)%>%
filter(Sesso=="M")%>%
summarise(mean(Peso), sd(Peso)),
"Tab.26 - Media e Dev.STD del peso di naeonati maschi (in 36 settimane di Gestazione)")| mean(Peso) | sd(Peso) |
|---|---|
| 2950.4 | 461.5958 |
##### osservazione 1306 (outlier)
kable_fun((dati%>%
filter(rownames(dati)=="1306")),
"Tab.27 - osservazione 1306 del dataset (outlier)")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 23 | 0 | 0 | 41 | 4900 | 510 | 352 | Nat | osp2 | F |
#ricerca misura Peso (settimana 41)
kable_fun((dati%>%
filter(Gestazione==41)%>%
filter(Sesso=="F")%>%
summarise(mean(Peso), sd(Peso))), "Tab.28 - Media e Dev.ST del peso di neonati maschi (in 41 settimane di gestazione)")| mean(Peso) | sd(Peso) |
|---|---|
| 3436.25 | 444.6357 |
ANALISI DEGLI OUTLIER/PUNTI INLFUENTI MOD3
Dal confronto dei risultati e delle analisi grafiche, emergono almeno tre outlier influenti sul modello, tra cui l’osservazione 1551, che supera la soglia di Cook (0.83) e risulta particolarmente critica per il modello 3. Di seguito si attenzionano gli outliers e i punti influenti sul modello 3:
Osservazione 155(Tab.26 -outlier) : registrazione di peso elevato di 3610g, superiore alla media (2950g ± 461g) dei neonati nati a 36 settimane come riportato in Tab.26;
Osservazione 1306 (Tab.27 - outlier): registrazione di peso elevato di 4900g, ben oltre la media 3436g ± 444g per neonati a 41 settimane, con una durata della gestazione eccezionalmente lunga, come riportato in Tab.28.
Osservazione 1551 (Tab.22 - outlier influente): peso anomalo registrato per la neonata di 4370g alquanto sospetto, rispetto alla media dei pesi delle neonate con lo stesso numero di settimane di gestazione (peso medio Neonate 3125.89+-390g, riportato in Tab.24) , con misure antropometriche del neonato “fuori dalla norma”. Infatti anche le altre caratteristiche antropometriche (lunghezza e diametro craniale) risultano incoerenti con il dato di peso-sesso della neonata 1551. Per una neonata nata con una gestazione 38 settimane rinveniamo in media le seguenti misure antropometriche (Tab.23):
Misure Neonate (gestazione 38) = 485.47+-22.67mm (Lunghezza) e 337.84+-12.65mm(Cranio)
Neonata 1551 = 315 mm (Lunghezza) e 374 mm (Cranio)
Nelle prossime sezioni proveremo di nuovo a costruire i modelli di regressione lineare multipla creati in precedenza (mod3,mod4 e mod5), eliminando l’osservazione 1551 dai modelli, considerato un punto influente secondo le osservazioni delle distanze di Cook (distanza di cook 0.83).
3.1.2 RICOSTRUZIONE MODELLO3_NEW (eliminazione osservazione 1551)
subset_dati2 <- dati%>%
filter(rownames(dati)!="1551")
modello3_new <- lm(data=subset_dati2,
Peso~N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso)
par(mfrow=c(2,2))
plot(modello3_new, which = 1, main = "(Verifica linearità/omoschedasticità)")
# Normal Q-Q
plot(modello3_new, which = 2, main = "(Verifica Normalità sui resisui)")
# Scale-Location
plot(modello3_new, which = 3, main = "(Verifica omoschedasticità)")
# Residuals vs Leverage
plot(modello3_new, which = 5, main = "(Osservazioni influenti sul modello)")
#leverage
lev <-hatvalues(modello3_new)
plot(lev,
main="Leverages\n (Spazio dei regressori)",
ylab="Distanza di Cook"
)
p=sum(lev)
soglia = 2*p/n
abline(h=soglia,col="red")
#outliers
plot(rstudent(modello3_new),
main="Outliers\n (spazio delle risposte)"
)
abline(h=c(-2,2),col="red")
#Cook's distance
cook <- cooks.distance(modello3_new)
plot(cook,
main="Verifica punti influenti\nCook's distance (scatterplot)",
ylab="Cook's distance"
)
abline(h=0.5,col="red")
#Cook's distance 2
plot(modello3_new, which = 4, main = "Verifica Punti influenti")
abline(h=0.5,col="red")
summary(modello3_new)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = subset_dati2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1165.68 -179.74 -12.42 162.92 1410.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6683.8326 133.1602 -50.194 < 2e-16 ***
## N.gravidanze 13.1448 4.2576 3.087 0.00204 **
## Gestazione 29.6341 3.7369 7.930 3.27e-15 ***
## Lunghezza 10.8899 0.3019 36.077 < 2e-16 ***
## Cranio 9.9192 0.4227 23.465 < 2e-16 ***
## SessoM 78.1376 10.9929 7.108 1.53e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7372, Adjusted R-squared: 0.7367
## F-statistic: 1398 on 5 and 2491 DF, p-value: < 2.2e-16Dall’eliminazione dell’osservazione influente 1551, il modello3_new (ex mod3) ha guadagnato in termini di \(R^2\) e di capacità di generalizzare i dati previsionali sul peso dei neonati, passando da un \(R^2 aggiustato\) di 0.7264 -> 0.7367. La modifica del modello è risultata soddisfacente, motivo per cui il modell3_new verrà utilizzato per eseguire previsioni del neonato.
kable_fun(round(vif(modello3_new),2), "Tab.29 - VIF DEL MODELLO3_NEW")| x | |
|---|---|
| N.gravidanze | 1.02 |
| Gestazione | 1.68 |
| Lunghezza | 2.13 |
| Cranio | 1.66 |
| Sesso | 1.04 |
3.1.3 RICOSTRUZIONE MODELLO4_NEW (eliminazione osservazione 1551)
Preso atto che l’osservazione 1551 risulta essere un punto influente per il modello, e ricalcando quanto fatto anche per il modello3_new, eliminiamo dal dataset di input al modello4_new l’osservazione in questione e ricostruiamo ex-novo il mod4 precedente, ovvero quello in cui si considerava il termine quadratico di \(Lunghezza^2\) all’interno della relazione lineare:
modello4_new <- lm(data=subset_dati2,
Peso~N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso+
I(Lunghezza^2))
par(mfrow=c(2,2))
plot(modello4_new, which = 1, main = "(Verifica linearità/omoschedasticità)")
# Normal Q-Q
plot(modello4_new, which = 2, main = "(Verifica Normalità sui resisui)")
# Scale-Location
plot(modello4_new, which = 3, main = "(Verifica omoschedasticità)")
# Residuals vs Leverage
plot(modello4_new, which = 5, main = "(Osservazioni influenti sul modello)")
#leverage
lev <-hatvalues(modello4_new)
plot(lev,
main="Leverages\n (Spazio dei regressori)",
ylab="Distanza di Cook"
)
p=sum(lev)
soglia = 2*p/n
abline(h=soglia,col="red")
#outliers
plot(rstudent(modello4_new),
main="Outliers\n (spazio delle risposte)"
)
abline(h=c(-2,2),col="red")
#Cook's distance
cook <- cooks.distance(modello4_new)
plot(cook,
main="Verifica punti influenti\nCook's distance (scatterplot)",
ylab="Cook's distance"
)
abline(h=0.5,col="red")
#Cook's distance 2
plot(modello4_new, which = 4, main = "Verifica Punti influenti")
abline(h=0.5,col="red")
eseguiamo il summary del modello4_new:
summary(modello4_new)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = subset_dati2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1176.71 -178.98 -11.61 164.25 1327.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.610e+03 7.589e+02 -2.121 0.033985 *
## N.gravidanze 1.418e+01 4.222e+00 3.358 0.000796 ***
## Gestazione 3.777e+01 3.893e+00 9.703 < 2e-16 ***
## Lunghezza -1.172e+01 3.343e+00 -3.505 0.000465 ***
## Cranio 1.015e+01 4.203e-01 24.146 < 2e-16 ***
## SessoM 7.220e+01 1.093e+01 6.606 4.8e-11 ***
## I(Lunghezza^2) 2.330e-02 3.431e-03 6.790 1.4e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.9 on 2490 degrees of freedom
## Multiple R-squared: 0.742, Adjusted R-squared: 0.7414
## F-statistic: 1193 on 6 and 2490 DF, p-value: < 2.2e-16kable_fun(round(vif(modello4_new),2), "Tab.30 - VIF DEL MODELLO4_NEW")| x | |
|---|---|
| N.gravidanze | 1.03 |
| Gestazione | 1.85 |
| Lunghezza | 266.44 |
| Cranio | 1.67 |
| Sesso | 1.05 |
| I(Lunghezza^2) | 255.52 |
Nel modello 4_new, eliminando l’osservazione 1551, riconosciuta in precedenza come influente, è possibile osservare che non vi sono più osservazioni di allarme nel grafico di cook, inoltre la distribuzione dei residui interno alla media di 0 è meglio descritta per i valori di pesi più inferiori, dove il modello basico 3 (privo di termini quadratici) non riusce bene a generalizzare le osservazione in quel range di peso.
Per concludere dall’eliminazione dell’osservazione influente 1551, il modello4_new (ex mod4) ha guadagnato in termini di \(R^2\) e di capacità di generalizzare i dati, passando da un \(R^2 aggiustato\) di 0.7289-> 0.7414. Rimane però un problema stabilità del modello, dato che il VIF per i termini di \(Lunghezza\) e \(Lunghezza^2\) rimane elevato e molto >> 5
3.1.3 RICOSTRUZIONE MODELLO5_NEW (eliminazione osservazione 1551)
Ripetiamo la stessa procedura di ricostruzione anche per il modello5_new, cioè il modello che considera l’effetto di interazione fra le variabili LUNGHEZZA x CRANIO del neonato, ed osserviamo i risultati di tale modifica:
modello5_new <- lm(data=subset_dati2,
Peso~N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso+
Lunghezza*Cranio)
par(mfrow=c(2,2))
plot(modello5_new, which = 1, main = "(Verifica linearità/omoschedasticità)")
# Normal Q-Q
plot(modello5_new, which = 2, main = "(Verifica Normalità sui resisui)")
# Scale-Location
plot(modello5_new, which = 3, main = "(Verifica omoschedasticità)")
# Residuals vs Leverage
plot(modello5_new, which = 5, main = "(Osservazioni influenti sul modello)")
#leverage
lev <-hatvalues(modello5_new)
plot(lev,
main="Leverages\n (Spazio dei regressori)",
ylab="Distanza di Cook"
)
p=sum(lev)
soglia = 2*p/n
abline(h=soglia,col="red")
#outliers
plot(rstudent(modello5_new),
main="Outliers\n (spazio delle risposte)"
)
abline(h=c(-2,2),col="red")
#Cook's distance
cook <- cooks.distance(modello5_new)
plot(cook,
main="Verifica punti influenti\nCook's distance (scatterplot)",
ylab="Cook's distance"
)
abline(h=0.5,col="red")
#Cook's distance 2
plot(modello5_new, which = 4, main = "Verifica Punti influenti")
abline(h=0.5,col="red")
eseguiamo il summary del modello5_new:
summary(modello5_new)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Lunghezza * Cranio, data = subset_dati2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1169.44 -179.23 -12.41 164.64 1444.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.584e+02 1.009e+03 0.157 0.875302
## N.gravidanze 1.390e+01 4.220e+00 3.293 0.001006 **
## Gestazione 3.739e+01 3.873e+00 9.655 < 2e-16 ***
## Lunghezza -3.831e+00 2.173e+00 -1.763 0.078012 .
## Cranio -1.161e+01 3.175e+00 -3.656 0.000262 ***
## SessoM 7.151e+01 1.094e+01 6.539 7.51e-11 ***
## Lunghezza:Cranio 4.428e-02 6.474e-03 6.840 9.96e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.9 on 2490 degrees of freedom
## Multiple R-squared: 0.742, Adjusted R-squared: 0.7414
## F-statistic: 1194 on 6 and 2490 DF, p-value: < 2.2e-16kable_fun(round(vif(modello5_new),2), "Tab.31 - VIF DEL MODELLO5_NEW")## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif| x | |
|---|---|
| N.gravidanze | 1.02 |
| Gestazione | 1.84 |
| Lunghezza | 112.59 |
| Cranio | 95.22 |
| Sesso | 1.05 |
| Lunghezza:Cranio | 322.19 |
Nel modello5_new, eliminando l’osservazione 1551, riconosciuta in precedenza come influente, è possibile osservare l’asseenza di punti di allarme nelle distanze di Cook. L’effetto di interazione fra le variabili Lunghezza e Cranio (assenti nei modelli base come modello3_new) migliora la capacità del modello cosi costruito di generalizzare verso i valori inferiori di Peso.
Concludendo, il modello5_new (ex mod5) ha guadagnato in termini di \(R^2\) e di capacità di generalizzare i dati, passando da un \(R^2 aggiustato\) di 0.7363 -> 0.7414. Rimane però un problema stabilità del modello, dato che l’effetto di interazione fra variabili introdotte nel modello porta il fattore VIF a valori molto >> 5 per i termini di \(Lunghezza e Cranio\).
3.3 Test statistici sulla parte residuale (modello3_new, modello4_new e modello5_new)
In questa sezione verrà trattata la parte diagnostica dei residui dei modelli costruiti nelle sezione precedenti, ma analizzati dal punto di ipotesi statistiche:
3.3.1 Test della normalità distributiva dei residui (Shapiro Wilk)
shapiro.test(residuals(modello3_new)) ##
## Shapiro-Wilk normality test
##
## data: residuals(modello3_new)
## W = 0.98891, p-value = 5.23e-13shapiro.test(residuals(modello4_new)) ##
## Shapiro-Wilk normality test
##
## data: residuals(modello4_new)
## W = 0.98994, p-value = 3.151e-12shapiro.test(residuals(modello5_new)) ##
## Shapiro-Wilk normality test
##
## data: residuals(modello5_new)
## W = 0.98881, p-value = 4.467e-13non è possibile accettare le condizioni di normalità per la distribuzione dei residui per i 3 modelli, in quanto il p-valore è prossimo allo zero, quindi rifiutiamo l’ipotesi H0 di distribuzione di tipo normale dei residui a favore dell’ipotesi alternativa H1.
dens3_new <- density(residuals(modello3_new))
plot(dens3_new,
main = "Kernel Density residuale \nmodello3_new vs Normale Teorica",
xlab = "Distribuzione dei residui",
ylab = "Densità", col = "blue", lwd = 2)
#dati della curva teorica
mean_resid <- mean(residuals(modello3_new))
sd_resid <- sd(residuals(modello3_new))
curve(dnorm(x,
mean = mean_resid,
sd = sd_resid),
col = "black",
lwd = 5,
add = TRUE)
legend("topright", legend = c("Modello 3new", "Normale teorica"),
col = c("blue", "black"), lwd = 4)
#modello4_new
dens4_new <- density(residuals(modello4_new))
plot(dens3_new,
main = "Kernel Density residuale \nmodello4_new vs Normale Teorica",
xlab = "Distribuzione dei residui",
ylab = "Densità", col = "green", lwd = 2)
mean_resid <- mean(residuals(modello4_new))
sd_resid <- sd(residuals(modello4_new))
curve(dnorm(x,
mean = mean_resid,
sd = sd_resid),
col = "black",
lwd = 5,
add = TRUE)
legend("topright", legend = c("Modello 4new", "Normale teorica"),
col = c("green", "black"), lwd = 4)
#modello5_new
dens5_new <- density(residuals(modello5_new))
plot(dens5_new,
main = "Kernel Density residuale \nmodello5_new vs Normale Teorica",
xlab = "Distribuzione dei residui",
ylab = "Densità", col = "red", lwd = 2)
mean_resid <- mean(residuals(modello5_new))
sd_resid <- sd(residuals(modello5_new))
curve(dnorm(x,
mean = mean_resid,
sd = sd_resid),
col = "black",
lwd = 5,
add = TRUE)
legend("topright", legend = c("Modello 5new", "Normale teorica"),
col = c("red", "black"), lwd = 4)
3.3 Test di omoschedasticità dei residui (test di Breush-Pagan)
bptest(modello3_new)##
## studentized Breusch-Pagan test
##
## data: modello3_new
## BP = 11.343, df = 5, p-value = 0.04499bptest(modello4_new)##
## studentized Breusch-Pagan test
##
## data: modello4_new
## BP = 15.429, df = 6, p-value = 0.01717bptest(modello5_new)##
## studentized Breusch-Pagan test
##
## data: modello5_new
## BP = 8.357, df = 6, p-value = 0.2131Nel modello3_new (senza interazione), il tramite test di Breush-Pagan restituisce un p-value di 0.045, approssimabile a 0.05 (pari ad un livello di significatività alpha impostato a 0.05). Non si rigetta l’ipotesi nulla H0 di varianza costanza dei residui, confermando omoschedasticità degli stessi nel modello lineare 3_new;
nel modello4_new (interazione quadratica Lunghezza), si prende atto che l’ipotesi nulla H0 di omoschedasticità (varianza costante) dei residui tramite test di Breush-Pagan viene rigettata a favore dell’ipotesi alternativa H1 di eteroschedasticità, perché il p-value è di 0.017, al di sotto del livello di significatività alpha di 0.05. Si conferma eteroschedasticità dei residui in questo modello.
Nel modello5_new non si rifiuta l’ipotesi di nulla H0 di omoschedasticità dei residui del modello per un p-valore di 0.2131. Ne consegue che i residui in questo modello presenza varianza costante.
3.4 Test di autocorrelazione dei residui (test di Durbin-Watson)
dwtest(modello3_new)##
## Durbin-Watson test
##
## data: modello3_new
## DW = 1.9536, p-value = 0.1227
## alternative hypothesis: true autocorrelation is greater than 0dwtest(modello4_new)##
## Durbin-Watson test
##
## data: modello4_new
## DW = 1.9495, p-value = 0.1032
## alternative hypothesis: true autocorrelation is greater than 0dwtest(modello5_new)##
## Durbin-Watson test
##
## data: modello5_new
## DW = 1.9539, p-value = 0.1246
## alternative hypothesis: true autocorrelation is greater than 0Tramite analisi Durbin-Watson, non viene rifiutata l’ipotesi H0 nulla di autocorrelazione dei residui pari 0, quindi il la distribuzione dei residui del modello 3_new, modello4_new e modello_new non riportano fenomeni di autocorrelazione fino a prova contraria. I modelli superano la diagnostica di Durbin-Watson, indicando che i resuidi dei relativi modelli non presentano fenomeni di autocorrelazione.
4. ANALISI DELLA QUALITA’ DEI MODELLI
deg.free <-as.data.frame(rbind(
round((mod1)$rank),
round((modello3_new)$rank),
round((modello4_new)$rank),
round((modello5_new)$rank)))
colnames(deg.free) <- c("Gradi di libertà")
R2.adj<- as.data.frame(rbind(round(summary(mod1)$adj.r.squared,3),
round(summary(modello3_new)$adj.r.squared,3),
round(summary(modello4_new)$adj.r.squared,3),
round(summary(modello5_new)$adj.r.squared,3)))
colnames(R2.adj) = c("R2 AGGIUSTATO")
RMSE<-as.data.frame(rbind(
round(summary(mod1)$sigma^2),
round(summary(modello3_new)$sigma^2),
round(summary(modello4_new)$sigma^2),
round(summary(modello5_new)$sigma^2)))
colnames(RMSE) = c("RMSE")
models_comparison <- as.data.frame(cbind(deg.free,
R2.adj,
RMSE))
rownames(models_comparison)<-c("modello1(base)","modello3_new","modello4_new","modello5_new")
kable_fun(models_comparison, "Tab.32 - Comparazione metriche previsionali dei modelli di regressione lineare multivariata (base e ottimizzati)")| Gradi di libertà | R2 AGGIUSTATO | RMSE | |
|---|---|---|---|
| modello1(base) | 11 | 0.728 | 75088 |
| modello3_new | 6 | 0.737 | 72549 |
| modello4_new | 7 | 0.741 | 71259 |
| modello5_new | 7 | 0.741 | 71240 |
I modelli ottimizzati, in particolare modello4_new e
modello5_new, evidenziano un buon equilibrio tra
semplificazione (riduzione dei gradi di libertà) e miglioramento delle
prestazioni, con un incremento del \(R^2
aggiustato\) e una riduzione dell’errore \(RMSE\). Tuttavia, entrambi i modelli
presentano problemi di multicollinearità, come indicato dall’analisi dei
\(VIF\). Al contrario, il
modello3_new rappresenta il miglior compromesso
tra performance e stabilità: pur avendo un \(R^2 aggiustato\) leggermente inferiore
rispetto agli altri modelli ottimizzati, mantiene un VIF generale sulle
variabili indipendenti < 5, dimostrando di non soffrire di
multicollinearità. Questo aspetto lo rende non solo potente, ma anche
affidabile per l’interpretazione e l’applicazione pratica.
5. PREVISIONE PESO DEL NEONATO
Una volta validati i modellI lI useremo per fare previsioni pratiche. Andiamo quindi a stimare stimare il peso di:
- una neonata (SESSO=“F”)
- nata da una madre alla terza gravidanza (N.gravidanze=2),
- che partorirà alla 39esima settimana(Gestazione =39).
#costruzione osservazione neonato
nuovo_neonato <- data.frame(N.gravidanze = 2,
Gestazione = 39,
Lunghezza = mean(Lunghezza),
Cranio = mean(Cranio),
Sesso = "F")
kable_fun(nuovo_neonato, "Tabella 33 - Dettagli sul nuovo neonato")| N.gravidanze | Gestazione | Lunghezza | Cranio | Sesso |
|---|---|---|---|---|
| 2 | 39 | 494.6958 | 340.0292 | F |
# Previsione del peso
peso_predetto_mod3_new<- as.data.frame(predict(modello3_new,
newdata = nuovo_neonato))
peso_predetto_mod4_new<- as.data.frame(predict(modello4_new,
newdata = nuovo_neonato))
peso_predetto_mod5_new<-as.data.frame(predict(modello5_new,
newdata = nuovo_neonato))
colnames(peso_predetto_mod3_new)=c("Modello3_new")
colnames(peso_predetto_mod4_new)=c("Modello4_new")
colnames(peso_predetto_mod5_new)=c("Modello5_new")
weights_comparison <-cbind(peso_predetto_mod3_new,
peso_predetto_mod4_new,
peso_predetto_mod5_new)
rownames(weights_comparison) <- c("Peso (grammi)")
kable_fun(weights_comparison, "Tabella 34 - Comparison previsionale modelli multivariati")| Modello3_new | Modello4_new | Modello5_new | |
|---|---|---|---|
| Peso (grammi) | 3258.184 | 3246.479 | 3250.72 |
6. VISUALIZZAZIONI
6.1 Impatto del fumo sul Peso del neonato fra settimane di gestazione e Peso del neonato (grammi)
ggplot(data=dati) +
geom_point(aes(x=Gestazione, y=Peso, col=as.factor(Fumatrici)), position = "jitter") +
geom_smooth(aes(x=Gestazione, y=Peso, col=as.factor(Fumatrici)), method="lm", se=FALSE) +
geom_smooth(aes(x=Gestazione, y=Peso), col="black", method="lm", se=FALSE) +
labs(x="Gestazione (settimane)", y="Peso neonato (grammi)",
color="Fumatrici") +
scale_color_manual(values=c("0"="blue", "1"="red"),
labels=c("Non fumatrici", "Fumatrici"))
6.2 - Relazione lineare fra settimane di gestazione e Peso del neonato (grammi)
La figura sottostante mostra la relazione lineare tra le settimane di gestazione e il peso alla nascita, distinti per sesso. Come atteso, il peso aumenta con l’avanzare delle settimane di gestazione. Le linee del modello evidenziano le differenze di crescita del peso tra maschi e femmine:
ggplot(data=dati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso), position = "jitter")+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso), se=F,method="lm")+
geom_smooth(aes(x=Gestazione,
y=Peso),
col="black", se=F, method="lm")+
labs(x="Gestazione (settimane)",
y="Peso neonato (grammi)")
6.3 - Relazione lineare fra la lunghezza (in millimetri) e il Peso del neonato
La figura sottostante mostra la relazione tra la lunghezza prenatale rilevata tramite ecografia e il peso alla nascita, distinti per sesso. Si osserva che all’aumentare della lunghezza prenatale aumenta il peso neonatale previsto. Le linee del modello evidenziano le differenze di crescita del peso tra maschi e femmine in funzione della lunghezza ecografica.
ggplot(data=dati)+
geom_point(aes(x=Lunghezza,
y=Peso,
col=Sesso), position = "jitter")+
geom_smooth(aes(x=Lunghezza,
y=Peso,
col=Sesso), se=F,method="lm")+
geom_smooth(aes(x=Lunghezza,
y=Peso),
col="black", se=F, method="lm")+
labs(x="Lunghezza neonato (millimetri)",
y="Peso neonato (grammi)")
6.4 - Relazione lineare fra la diametro craniale (in millimetri) e il Peso del neonato
La figura sottostante mostra la relazione tra il diametro craniale prenatale rilevato tramite ecografia e il peso alla nascita, distinti per sesso. Si osserva che all’aumentare del diametro craniale prenatale aumenta il peso neonatale previsto. Le linee del modello evidenziano le differenze di crescita del peso tra maschi e femmine in funzione del diametro craniale ecografico.
ggplot(data=dati)+
geom_point(aes(x=Cranio,
y=Peso,
col=Sesso), position = "jitter")+
geom_smooth(aes(x=Cranio,
y=Peso,
col=Sesso), se=F,method="lm")+
geom_smooth(aes(x=Cranio,
y=Peso),
col="black", se=F, method="lm")+
labs(x="Diametro craniale (millimetri)",
y="Peso neonato (grammi)")
library(plotly)
# Grafico interattivo con 4 dimensioni
plot_ly(
x = ~Lunghezza,
y = ~Cranio,
z = ~Peso,
color = ~Gestazione,
type = "scatter3d",
mode = "markers") %>%
layout(title="Modellazione 3D interattivo delle variabili prenatali - \n in funzione delle settimane di gestazione",
scene=list(
xaxis = list(title = "Lunghezza (mm)"),
yaxis = list(title = "Cranio (mm)"),
zaxis = list(title = "Peso (grammi)")),
colorbar(list(title="Gestazione della madre (settimane)",
side="top")))Il grafico 3D interattivo rappresenta la relazione tra tre variabili prenatali: il peso alla nascita in grammi( asse y), la lunghezza del neonato in millimetri (asse x) e il diametro craniale in millimetri (asse z). I punti sono colorati in base alle settimane di gestazione, dal viola (settimane minori) al giallo (settimane maggiori). Si osserva una relazione positiva tra peso, lunghezza e diametro craniale, con neonati nati da gestazione più lunghe (gialli) che tendono a presentare misure antropometriche più elevate, ed un peso maggiore. La distribuzione dei dati mostra una concentrazione centrale, indicando che la maggior parte dei neonati rientra in intervalli tipici per queste variabili.
7. CONSIDERAZIONI FINALI
Lo studio ha evidenziato che le principali variabili che influenzano la predizione del peso alla nascita del neonato sono il numero di settimane di gestazione, la lunghezza prenatale del neonato misurata tramite ecografia e il diametro craniale ecografico. In particolare, all’aumentare delle settimane di gestazione, aumenta anche il peso previsto del neonato. ll sesso dei neonati risulta anch’esso rilevante, con i neonati maschi che, in media, pesano di più rispetto alle femmine. Sebbene debole, il numero di gravidanze precedenti della madre si è dimostrato significativo nei modelli multivariati costruiti (modello3_new, modello4_new e modello5_new), pur non essendo emersa una relazione diretta con il peso (come ampiamente trattato nella sezione 1.4.1 di questo progetto). Sulla base del dataset fornitomi per questo progetto (neonati.csv), non è stata trovata alcuna associazione tra le variabili “Anni della madre” e “Fumatrici” con il peso del neonato, come inizialmente supposto nel testo della consegna di questo progetto. Infine, analizzando le correlazioni bi-variate, è emerso che la lunghezza del neonato presenta una correlazione quadratica con il peso. Inoltre, le variabili di controllo come lunghezza e diametro craniale sono correlate positivamente fra loro, cosi come ognuna di esse è correlata, a sua volta, linearmente e positivamente con il numero di settimane di gestazione.