1) Introduzione

Modello Statistico per la Previsione del Peso Neonatale Contesto Aziendale Azienda: Neonatal Health Solutions Obiettivo: Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario.

Raccolta dei Dati e Struttura del Dataset: import librerie, Set del wd ,Caricamento del dataset

2) Analisi esplorativa

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

2.1) Variabili Quantitative:

Anni.madre,N.gravidanze,Gestazione,Peso,Lunghezza,Cranio

#Funzione per calcolare la moda
moda_fun <- function(x) {
  moda_value <- names(table(x))[which.max(table(x))]
  return(moda_value)}

# Funzione principale per gli indici di posizione
 pos_index <- function(x) {
  
  quant <- quantile(x,na.rm=T)
  return(c(Media=round(mean(x, na.rm=T),2),
           Mediana=round(median(x, na.rm=T),2),
           Moda= moda_fun(x), #uso la funzione precedentemente creata
           Min=round(quant[1],2), #indicizzo dentro quant
           Q1= round(quant[2],2),
           Q3= round(quant[4],2), #ho tolto la moda- indice 4 corrisponde a Q3
           Max=round(quant[5],2)))} 


# Calcolo indici posizione usando la funzione pos_index e sapply
   #Creo un vettore c con i nomi e i valori da usare con sapply

var_quanti <- c('Anni.madre','N.gravidanze','Gestazione','Peso','Lunghezza', "Cranio")

indici_posizione_var <-sapply(dati_clean[var_quanti],pos_index) #da "dati_clean" prendo il vettore "c" e si applica la funz. pos_index

#Trasposizione delle righe con le colonne, posiziono le variabili sulle colonne
tab_index_pos <-as.data.frame(t(indici_posizione_var))
#uso kable
tab_index_pos %>%
  knitr::kable(caption="Tabella indici di posizione per le variabili quantitative di interesse", digits=2, align='c')

Tabella indici di posizione per le variabili quantitative di interesse
	Media	Mediana	Moda	Min.0%	Q1.25%	Q3.75%	Max.100%
Anni.madre	28.19	28	30	13	25	32	46
N.gravidanze	0.98	1	0	0	0	1	12
Gestazione	38.98	39	40	25	38	40	43
Peso	3284.18	3300	3300	830	2990	3620	4930
Lunghezza	494.7	500	500	310	480	510	565
Cranio	340.03	340	340	235	330	350	390

# variabilità per le variabili quantitative 
#creo le funzioni per calcolare gli indici di variabilità - deviazione standard e interquartile range (IQR) - coefficiente di variabilità

#funzione per il coefficiente di variazione CV
cv_fun <- function(x) {
  # Controlla per la divisione per zero e gestisce NA
  mu <- mean(x, na.rm = TRUE)
  sigma <- sd(x, na.rm = TRUE)
  if (is.na(mu) || mu == 0) return(NA)
  return(round((sigma / mu) * 100, 2))} # ritorna in percentuale


#funzione per indici di variabilità
var_index <-function(x){
  dev_st <- sd(x,na.rm = T)
  data_range <- max(x,na.rm=T) - min(x,na.rm = T)
  iqr_range <-IQR(x, na.rm=T)
  
  return(c( D_stand =round(sd(x,na.rm=T),2),
            Range=round(max(x,na.rm=T) - min(x,na.rm=T),2),
            IQR=round(IQR(x,na.rm=T),2),
            CV=cv_fun(x) ))} 



### CALCOLO DELLA VARIABILITA' e TABELLE

tab_Variability <- as.data.frame(t(sapply(dati_clean[var_quanti],var_index)))
knitr::kable( tab_Variability,caption= " INDICI DI VARIABILITA'", digits = 2, align='c')

INDICI DI VARIABILITA’
	D_stand	Range	IQR	CV
Anni.madre	5.22	33	7	18.51
N.gravidanze	1.28	12	1	130.50
Gestazione	1.87	18	2	4.79
Peso	525.23	4100	630	15.99
Lunghezza	26.33	255	30	5.32
Cranio	16.43	155	20	4.83

#FUNZIONI PER CALCOLARE LA FORMA DELLA DISTRIBUZIONE E RELATIVI CALCOLI E PLOT 
# variabili utilizzate: qualitative

#creo funzioni per calcolo della forma:(funziona solo con pacchetto Moments! runnato all'inizio dello script)
forma_index <-function(x) {
  asimmetria <- skewness(x,na.rm=T)
  curtosi <- kurtosis(x,na.rm=T) -3 # -3 per calcolare la cursoti in eccesso e calcolare il valore partendo dallo zero. 
  
  return(c( Asimmetria=round(asimmetria,3),
            Curtosi=round(curtosi,3)))}


#CALCOLO DELLA FORMA E TABELLE

tab_forma <- as.data.frame(t(sapply(dati_clean[var_quanti], forma_index)))
knitr::kable( tab_forma, caption= "Forma di distribuzione- asimmetria e curtosi",align="c")

Forma di distribuzione- asimmetria e curtosi
	Asimmetria	Curtosi
Anni.madre	0.151	-0.106
N.gravidanze	2.513	10.982
Gestazione	-2.065	8.256
Peso	-0.647	2.029
Lunghezza	-1.515	6.481
Cranio	-0.785	2.945

#PLOT CON CICLO FOR DELLE TABELLE PER LE VARIABILI QUANTITATIVE DI INTERESSE. 

# Ciclo for per il plot densità per ogni variabile
for (var in var_quanti) {
  
  #Calcolo della media e mediana per il plot
media_val <-mean(dati_clean[[var]],na.rm=T)
mediana_val <-median(dati_clean[[var]],na.rm=T)
asim_val    <- round(skewness(dati_clean[[var]], na.rm = TRUE), 2)

  
  # crea il grafico
  p<-ggplot(data=dati_clean,aes(x=.data[[var]])) +
    # Curva densità
    geom_density(fill = "lightblue", alpha = 0.5,color = NA) + 
    
    # Linea della media rossa tratteggiata
    geom_vline(aes(xintercept = media_val, 
               color = 'Media'), 
               linetype = 'dashed', 
               linewidth= 1) + 
    geom_vline(aes(xintercept = mediana_val, 
               color = "Mediana"), 
               linetype = "dotted", linewidth = 1)+
    
    #assegno colori per la legenda
    scale_color_manual(name = "Index", 
                       values = c("Media" = "red", "Mediana" = "darkgreen")) +
    
    # Titoli dei grafici
    labs(
      title = paste("Distribuzione:", var),
      subtitle = paste("asimmetria:", asim_val), #contiene calcolo asimmetria
      x = var, 
      y = "densità"
    ) +
    theme_minimal()+
    theme(legend.position = "bottom")
  
  # Stampa i grafici
  print(p)}

2.2) Variabili Qualitative

SESSO, OSPEDALE, FUMATRICI, TIPO.PARTO

#Funzione per calcolare indice di Gini normalizzato

gini_normalized_fun <- function(x) {
  
  frequenze_relative <- prop.table(table(x))
  xx <- length(unique(x))
  # Gestione se k <= 1 (Gini non è definito)
  if (xx <= 1) return(NA)
  # Calcolo
  gini_non_norm <- 1 - sum(frequenze_relative^2)
  gini_norm <- gini_non_norm / ((xx - 1) / xx)
  
  return(round(gini_norm, 2))
}


#Vettore contenente i nomi delle variabili qualitative
var_quali <- c("Sesso","Ospedale","Tipo.parto","Fumatrici")

# Calcolo dell'indice di Gini solo per le variabili qualitative
result <- sapply(dati_clean[var_quali], gini_normalized_fun)

tab_gini <- data.frame(Gini = result)

knitr::kable(tab_gini, caption = "Indice di Gini Normalizzato per Variabile",align = 'c')

Indice di Gini Normalizzato per Variabile
	Gini
Sesso	1.00
Ospedale	1.00
Tipo.parto	0.83
Fumatrici	0.16

for (var in var_quali) {
  pp <- ggplot(data=dati_clean,aes(x=.data[[var]]))+
    geom_bar(fill="lightblue")+
    labs(title = paste("barplot variabili categoriche:",var),
x= var,
y="conteggio")+
    theme_minimal()
  print(pp)
}

3) Ipotesi

Inoltre si saggeranno le seguenti ipotesi con i test adatti:

In alcuni ospedali si fanno più parti cesarei
La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione
Le misure antropometriche sono significativamente diverse tra i due sessi

3.1) Ipotesi 1

Test- 1 - Chi-quadro

H₀: la proporzione di parti cesarei è uguale nei tre ospedali

H₁: la proporzione di parti cesarei differisce in almeno un ospedale

tab_contingenza <-table(dati_clean$Ospedale,dati_clean$Tipo.parto)
X=chisq.test(tab_contingenza)

#tabelle
knitr::kable(tab_contingenza,caption="Tabella di contingenza: Ospedale vs tipo di parto", align = 'c')

Tabella di contingenza: Ospedale vs tipo di parto
	Ces	Nat
osp1	242	574
osp2	254	594
osp3	232	602

#frequenze attese
knitr::kable(X$expected,caption = "frequenze attese",
             digits=2,
             align = 'c')

frequenze attese
	Ces	Nat
osp1	237.81	578.19
osp2	247.14	600.86
osp3	243.06	590.94

#risultati chi-quadro
tab_risultato <- data.frame(
  Statistica = round(X$statistic, 3),
  Df = X$parameter,
  P_value = round(X$p.value, 4)
)

knitr::kable(tab_risultato,
             caption = "Risultato test chi-quadro",
             align = 'c')

Risultato test chi-quadro
	Statistica	Df	P_value
X-squared	1.083	2	0.5819

Conclusione: Con p-value = 0.58 > 0.05, non rifiutiamo H₀. Non c’è evidenza statistica che la proporzione di parti cesarei differisca tra i tre ospedali.

3.2) Ipotesi 2

###Test- 2 - t-test 1 group

Per le ipotesi sono stati ricavati i valori di peso e lunghezza media da letteratura medica.

Riferimenti: WHO Child Growth Standards: Methods and development: Length/height-for-age, weight-for-age, weight-for-length, weight-for-height and body mass index-for-age. (2006).

H₀: Il peso/ la lunghzezza medio/a del campione è uguale a quello della popolazione media popolazione di riferimento: μ_peso = 3300g , μ_lunghezza = 500mm

H₁: Il peso/la lunghezza medio/a del campione è diversa da quella della popolazione media popolazione di riferimento: μ_peso ≠ 3300g μ_lunghezza ≠ 500mm

#t-test PESO
t_peso=t.test(dati_clean$Peso,mu=3300, alternative="two.sided")
t_peso

## 
##  One Sample t-test
## 
## data:  dati_clean$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

#t-test Lunghezza
t_lunghezza=t.test(dati_clean$Lunghezza,mu=500, alternative = "two.sided")
t_lunghezza

## 
##  One Sample t-test
## 
## data:  dati_clean$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6628 495.7287
## sample estimates:
## mean of x 
##  494.6958

#assunsioni di normalitò del t-test
#i dati sono realmente normali e indipendenti? 
#Shapiro-Wilk
shapiro.test(dati_clean$Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  dati_clean$Peso
## W = 0.97068, p-value < 2.2e-16

shapiro.test(dati_clean$Lunghezza)

## 
##  Shapiro-Wilk normality test
## 
## data:  dati_clean$Lunghezza
## W = 0.90944, p-value < 2.2e-16

#tab t-test PESO
tab_ttest_peso <- data.frame(
  Stat_t = round(t_peso$statistic, 3),
  Df = t_peso$parameter,
  #P_val = round(t_peso$p.value, 4),
  P_val = ifelse(t_peso$p.value < 0.001, "< 0.001", 
               round(t_peso$p.value, 4)),
  IC_low = round(t_peso$conf.int[1], 2),
  IC_up = round(t_peso$conf.int[2], 2)
)

knitr::kable(tab_ttest_peso, 
             caption = "Risultato t-test Peso",
             align = 'c')

Risultato t-test Peso
	Stat_t	Df	P_val	IC_low	IC_up
t	-1.505	2497	0.1324	3263.58	3304.79

#tab t-test Lunghezza
tab_ttest_lunghezza <- data.frame(
  Stat_t =round(t_lunghezza$statistic,3),
  Df = t_lunghezza$parameter,
  #P_val= round(t_lunghezza$p.value,4),
  P_val = ifelse(t_lunghezza$p.value < 0.001, "< 0.001", 
               round(t_lunghezza$p.value, 4)),
  IC_low= round(t_lunghezza$conf.int[1],2),
  IC_up=round(t_lunghezza$conf.int[2],2))

knitr::kable(tab_ttest_lunghezza, 
             caption = "Risultato t-test Lunghezza",
             align = 'c')

Risultato t-test Lunghezza
	Stat_t	Df	P_val	IC_low	IC_up
t	-10.069	2497	< 0.001	493.66	495.73

Conclusione: i neonati del campione hanno una lunghezza media significativamente inferiore alla media della popolazione (p-value < 0.05), mentre il peso è in linea, non differisce significativamente dalla popolazione (p-value = 0.1324).

Le assunzioni dei t-test sono state controllate con test della normalità e indipendenza:

Il test di Shapiro-Wilk rifiuta l’ipotesi di normalità (p < 0.05). Tuttavia, con un campione di 2498 osservazioni, il test è eccessivamente sensibile a piccole deviazioni dalla normalità. Visivamente i density plot mostrano distribuzioni approssimativamente normali, e grazie al Teorema del Limite Centrale, con n > 30 la distribuzione campionaria della media converge alla normalità. Il t-test rimane quindi valido e robusto. Per un supporto visivo consulare i grafici di distrubizione delle variabili Peso e Lunghezza.

3.3) Ipotesi 3

H₀: : Le misure antropometriche non sono significativamente diverse tra i due sessi H₁: : Le misure antropometriche sono significativamente diverse tra i due sessi

###Test- 3 - t-test 2 sample

Si eseguono t-test a due gruppi per confrontare le differenze delle variabili antropometriche (Peso, Lunghezza, Cranio) tra i due gruppi (Maschi e Femmine).

Si esegue test dell’omogeneità per verificare se cambia la varianza tra i due gruppi e scegliere il t-test adatto tra t-student e Welch. t-student verrà eseguito per le variabili che mostrano varianza uguale (p >0.05 al test di Levene).

#libreria car per usare il test omogeneità della varianza Levene
#sulla base di questo scego il t-test adatto (t student o Welch)
#t student se varianza è uguale altrimenti Welch se varianza è diversa

leveneTest(Peso ~ Sesso, data = dati_clean)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  0.8222 0.3646
##       2496

leveneTest(Lunghezza ~ Sesso, data = dati_clean)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value   Pr(>F)   
## group    1  10.571 0.001164 **
##       2496                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

leveneTest(Cranio~ Sesso, data = dati_clean)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  1.2098 0.2715
##       2496

#t-test a due gruppi per variabili antropometriche
#t-test Welch variabile Lunghezza,  t-student per variabili Peso e Cranio

t_2_Peso=t.test(Peso ~ Sesso, data=dati_clean, var.equal=TRUE) #t-student
t_2_Lunghezza=t.test(Lunghezza ~ Sesso, data=dati_clean, var.equal=FALSE) #Welch
t_2_Cranio=t.test(Cranio ~ Sesso, data=dati_clean, var.equal=TRUE) #t-student

#tabelle:
#Peso
tab_t_test2_Peso <- data.frame(
  Stat_t =round(t_2_Peso$statistic,3),
  Df = t_2_Peso$parameter,
  #P_val= round(t_2_Peso$p.value,4),
  P_val = ifelse(t_2_Peso$p.value < 0.001, "< 0.001", 
               round(t_2_Peso$p.value, 4)),
  IC_low= round(t_2_Peso$conf.int[1],2),
  IC_up=round(t_2_Peso$conf.int[2],2))
#Lunghezza
tab_t_test2_Lunghezza <- data.frame(
  Stat_t =round(t_2_Lunghezza$statistic,3),
  Df = t_2_Lunghezza$parameter,
  #P_val= round(t_2_Lunghezza$p.value,4),
  P_val = ifelse(t_2_Lunghezza$p.value < 0.001, "< 0.001", 
               round(t_2_Lunghezza$p.value, 4)),
  IC_low= round(t_2_Lunghezza$conf.int[1],2),
  IC_up=round(t_2_Lunghezza$conf.int[2],2))
#Cranio
tab_t_test2_Cranio <- data.frame(
  Stat_t =round(t_2_Cranio$statistic,3),
  Df = t_2_Cranio$parameter,
  #P_val= round(t_2_Cranio$p.value,4),
  P_val = ifelse(t_2_Cranio$p.value < 0.001, "< 0.001", 
               round(t_2_Cranio$p.value, 4)),
  IC_low= round(t_2_Cranio$conf.int[1],2),
  IC_up=round(t_2_Cranio$conf.int[2],2))

#colonne
tab_t_test2_Peso$Variabile <- "Peso"
tab_t_test2_Lunghezza$Variabile <-"Lunghezza"
tab_t_test2_Cranio$Variabile <- "Cranio"

tab_completa_t_test_2 <- rbind(tab_t_test2_Peso,tab_t_test2_Lunghezza,tab_t_test2_Cranio)

knitr::kable(tab_completa_t_test_2, 
             caption = "Risultato t-test 2 sample Variabili antropometrice vs Sesso",
             align = 'c')

Risultato t-test 2 sample Variabili antropometrice vs Sesso
	Stat_t	Df	P_val	IC_low	IC_up	Variabile
t	-12.111	2496.000	< 0.001	-287.50	-207.37	Peso
t1	-9.582	2457.301	< 0.001	-11.94	-7.88	Lunghezza
t2	-7.434	2496.000	< 0.001	-6.11	-3.56	Cranio

Conclusione: Per tutte e tre le misure antropometriche, i maschi hanno valori significativamente più alti delle femmine (p-value < 2.2e-16). Questo è coerente con la letteratura medica perchè i neonati maschi tendono ad essere più grandi delle femmine alla nascita.

4) 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 ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

variabile dipendente Y = Peso variabili indipendenti X = Anni.madre, N.gravidanze, Gestazione, Lunghezza e Cranio ( valori ottenuto da ecografia prima del parto) Fumatrici, Sesso

Le variabili: Tipo.parto non inclusa per casualità temporale. Ospedale: verrà utilizzata come effetto fisso, variabile di controllo per controllare se esistono differenze tra ospedali come valutazione delle pratiche ospedaliere.

4.1) Matrice di correlazione

In questa sessione verrà svolta una matrice di correlazione tra la variabile peso (y) e possibili varaibili indipendenti (x). Osserviamo quali variabili sono più correlate con il Peso e identifico la possibile multicollinearità tra i predittori (le x).

matrice_cor <- cor(dati_clean[var_quanti])
matrice_cor

##               Anni.madre N.gravidanze Gestazione         Peso   Lunghezza
## Anni.madre    1.00000000  0.383317212 -0.1349423 -0.023781381 -0.06496158
## N.gravidanze  0.38331721  1.000000000 -0.1015306  0.002277118 -0.06046524
## Gestazione   -0.13494228 -0.101530623  1.0000000  0.591959182  0.61915103
## Peso         -0.02378138  0.002277118  0.5919592  1.000000000  0.79604149
## Lunghezza    -0.06496158 -0.060465242  0.6191510  0.796041486  1.00000000
## Cranio        0.01620505  0.038826618  0.4607450  0.704843839  0.60343791
##                  Cranio
## Anni.madre   0.01620505
## N.gravidanze 0.03882662
## Gestazione   0.46074505
## Peso         0.70484384
## Lunghezza    0.60343791
## Cranio       1.00000000

corrplot(matrice_cor,
         method = "color", #usa i colori
         type="upper", #mostra solo parte superiore simmetrica
         addCoef.col="black", #aggiunge valori numerici
         tl.col="black") #crea etichette nere

Conclusione: Le variabili che correlano meglio con Peso sono: Gestazione (0.59 moderata/forte, correlazione positiva), Lunghezza (0.80 forte correlazione positiva), Cranio (0.70 forte correlazione positiva), Anni.madre e N.graviande: assenza di correlazione (-0.02 , 0.00).

4.2) Modello Completo

modello predittivo Creo un modello preditto di regressione lineare multipla inserendo prima tutte le variabili: Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio, Fumatrici, Sesso) Verranno esplorati i p_value , icoefficienti e la bontà del modello R(Adjusted R-squared). Tenendo fisse tutte le altre variabili, un aumento di 1 unità in X produce una variazione di β grammi nel Peso

#modello completo con tutte le variabili

main_model <-lm(Peso ~ Anni.madre + N.gravidanze + Gestazione + Lunghezza+Cranio + Fumatrici+Sesso,data=dati_clean )
summary(main_model)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
##     Lunghezza + Cranio + Fumatrici + Sesso, data = dati_clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.6  -181.3   -15.7   163.6  2630.7 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6712.2405   141.3339 -47.492  < 2e-16 ***
## Anni.madre       0.8803     1.1491   0.766    0.444    
## N.gravidanze    11.3789     4.6767   2.433    0.015 *  
## Gestazione      32.9472     3.8288   8.605  < 2e-16 ***
## Lunghezza       10.2316     0.3011  33.979  < 2e-16 ***
## Cranio          10.5198     0.4271  24.633  < 2e-16 ***
## Fumatrici1     -30.3958    27.6080  -1.101    0.271    
## SessoM          78.0787    11.2132   6.963 4.24e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic: 948.3 on 7 and 2490 DF,  p-value: < 2.2e-16

Conclusione: Le variabili che correlano in modo significativo (p_value <0.05) sono N.gravidanze, Gestazione, Lunghezza, Cranio e Sesso maschile. Il modello ha un buon fit e spiega il 0.72 di variabilità di Y. Tenendo fisse tutte le altre variabili, un aumento di 1 unità in X produce una variazione di β grammi nel Peso: Ogni settimana di gestazione in più il peso aumenta di 32.95g Ad ogni mm in più di lunghezza il peso aumenta di 10.23g Ad ogni mm in più di diametro del cranio il peso aumenta di 10.52g I maschi pesano in media 78.08g in più delle femmine

4.3) Selezione modello (AIC)

Multicollinearità Controllo della multicollinearità tra variabili tramite il VIF. I valori sono < 5 , escludo che nel modello ci sia multicollinearità tra le variabili, i coefficienti del modello sono stabili e affidabili.

vif(main_model)

##   Anni.madre N.gravidanze   Gestazione    Lunghezza       Cranio    Fumatrici 
##     1.189264     1.187447     1.694331     2.079749     1.628987     1.006692 
##        Sesso 
##     1.040493

4.4) Modello ottimale

Trade off del modello con AIC e BIC

Applico la funzione stepAIC per verificare le penalità delle varie combinazioni del modello. Individuo la combinazione delle variabili che mi restituisce il modello ottimale. Il modello ottimale deve essere parsimonioso e con AIC più basso.

Dall’analisi risulta che l’optimal model comprende 5 variabili anzichè 7 1.N.Gravidanze, 2.Gestazione, 3.Lunghezza, 4.Cranio, 5.Sesso

Le variabili eliminate dal modello che facevano aumentare il valore AIC erano Anni.madre e Fumatrici, le stesse variabili che non mostravano significatività, p-value 0.444 e 0.271. Pertanto vengono eliminate dal modello perchè non aggiungono informazioni utili.

AIC_model<-stepAIC(main_model, direction = "both")

## Start:  AIC=28064.05
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Fumatrici + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     44292 187949505 28063
## - Fumatrici     1     91474 187996688 28063
## <none>                      187905214 28064
## - N.gravidanze  1    446756 188351970 28068
## - Sesso         1   3658879 191564093 28110
## - Gestazione    1   5587942 193493156 28135
## - Cranio        1  45789523 233694736 28607
## - Lunghezza     1  87128339 275033553 29014
## 
## Step:  AIC=28062.64
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Fumatrici + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     92548 188042054 28062
## <none>                      187949505 28063
## + Anni.madre    1     44292 187905214 28064
## - N.gravidanze  1    643981 188593487 28069
## - Sesso         1   3666800 191616305 28109
## - Gestazione    1   5544825 193494331 28133
## - Cranio        1  46056754 234006260 28608
## - Lunghezza     1  87116561 275066067 29012
## 
## Step:  AIC=28061.87
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188042054 28062
## + Fumatrici     1     92548 187949505 28063
## + Anni.madre    1     45366 187996688 28063
## - N.gravidanze  1    621053 188663107 28068
## - Sesso         1   3650790 191692844 28108
## - Gestazione    1   5477493 193519547 28132
## - Cranio        1  46098547 234140601 28608
## - Lunghezza     1  87532691 275574744 29015

#Optimal Model La differenza di R² è minima (0.7264 → 0.7265). Eliminando le variabili non significative, il modello è diventato più parsimonioso, la bontà del modello non è cambiata ma a parità di risultato e potere esplicativo viene preferito questo modello più parsimonioso.

Optimal_model <- lm(Peso ~ N.gravidanze + Gestazione + 
                        Lunghezza + Cranio + Sesso, 
                        data = dati_clean)
summary(Optimal_model)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati_clean)
## 
## 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-16

4.5 Analisi dei Residui

Residuals vs Fitted: La linea rossa è orizzontale e i residui sono dispersi in modo uniforme suggerendo l’omoschedasticità del modello.

Q-Q Plot: I residui seguono la diagonale nella parte centrale, le code tendono a variare.

Scale-Location: La linea mostra una lieve pendenza decrescente, suggerendo una leggera eteroschedasticità per valori bassi di peso

Residuals vs Leverage — Valori Influenti: L’analisi della Cook’s distance ha identificato l’osservazione 1551 come particolarmente influente (Cook’s D ≈ 0.7, soglia = 4/n = 0.0016).

Esaminando i valori di questa osservazione sembra che per una gestazione di39 settimane il peso 2740g è basso per l’età gestazionale

è stato svolto un confronto tra i due modelli ( con e senza questa osservazione) e mostrava una differenza minima in R² (0.7270 vs 0.7276). Si è quindi scelto di mantenere l’osservazione nel modello finale. l’informazione è documentata ma non riportato il codice.

par(mfrow = c(2,2))
plot(Optimal_model)

#plotting generale Cook's distance
plot(Optimal_model, which = 4)

# Calcola la Cook's distance per tutte le osservazioni
cooks <- cooks.distance(Optimal_model)

5) Previsioni:

Stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

Dati:

Sesso: Femmina, N.gravidanze: 2, Gestazione: 39

per la lunghezza e circonferenza del cranio in assenza di informazioni specifiche, il valore più rappresentativo è la media.

Per la previsione uso la funzione “predict”: questa usa i coefficienti del modello già stimato per calcolare il valore di Y per un nuovo caso.

Conclusione:

il modello prevede un peso di 3258g per questa neonata. Con il 95% di probabilità, il peso effettivo sarà compreso tra 2719g e 3797g. Il valore puntuale è in linea con la media del campione (3284g).

#Ricavo i valori medi di Lunghezza e Cranio

#creo data frame con i valori delle variabili del modello 
case_study <- data.frame(
  N.gravidanze = 2,
  Gestazione = 39,
  Lunghezza = round(mean(dati_clean$Lunghezza), 1),
  Cranio = round(mean(dati_clean$Cranio), 1),
  Sesso = "F"
)

#Funzione predict per fare previsione del modello
prevision_model <-predict(Optimal_model, #modello già stimato
        newdata = case_study, #nuovocaso
        interval = "prediction") #intervalli da calcolare


#tabella con kaable
tab_previsione <- data.frame(
  Caso = "Neonata, 3a gravidanza, 39 settimane",
  Peso_previsto = round(prevision_model[1], 0),
  IC_lower = round(prevision_model[2], 0),
  IC_upper = round(prevision_model[3], 0)
)

knitr::kable(tab_previsione,
             caption = "Previsione peso neonatale",
             align = 'c')

Previsione peso neonatale
Caso	Peso_previsto	IC_lower	IC_upper
Neonata, 3a gravidanza, 39 settimane	3258	2720	3797

Funzione predici peso

predict_peso <- function(n_gravidanze, gestazione, 
                          lunghezza, cranio, sesso) {
  case_study <- data.frame(
    N.gravidanze = n_gravidanze,
    Gestazione = gestazione,
    Lunghezza = lunghezza,
    Cranio = cranio,
    Sesso = sesso)
  
  predict(Optimal_model, 
          newdata = case_study,
          interval = "prediction")}

predict_peso(2,39,495,340,"F") # inserire i valori a disposizione del modello: N.parti, mese gestazione, lunghezza, circonferenza cranio e sesso.

##        fit      lwr      upr
## 1 3261.536 2722.583 3800.489

6) Visualizzazioni e conclusione

1)Plot n.1 mostra la relazione tra il peso osservato e il peso predetto dal modello

2)Plot n.2 mostra la relazione positiva tra settimane di gestazione e peso alla nascita. All’aumentare delle settimane di gestazione il peso aumenta in modo lineare.

3)Plot n.3 mostra le differenze di peso dei neonati delle madri fumatrici vs non fumatrici.
C’è una lieve differenza nel peso tra figli di madri fumatrici e non fumatrici. Però le fumatrici rappresentano solo una piccola parte del campione (4% ) e limita il potere statistico dell’analisi. Inoltre, il coefficiente nel modello non risulta statisticamente significativo (p = 0.271). Per questo motivo il predittore è stato tolto dal modello finale. In futuro, un campione più bilanciato permetterebbe di stimare l’effetto del fumo con maggiore precisione

4)Plot n.4 mostra la differenza di peso alla nascita tra neonati maschi e femmina. Come visto anche dal t-test e dal modello di regressione: i neonati maschi hanno un peso medio significativamente superiore alle femmine (+78g, p < 2.2e-16)

#1 PESO OSSERVATO VS IL PESO PREDETTO DAL MODELLO
ggplot(dati_clean, aes(x = fitted(Optimal_model), 
                        y = Peso)) +
  geom_point(alpha = 0.3, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, 
              color = "red", linewidth = 1) +
  labs(title = "Valori predetti vs osservati",
       x = "Peso predetto (g)",
       y = "Peso osservato (g)") +
  theme_minimal()

#2 GESTAZIONE-PESO
ggplot(dati_clean, aes(x = Gestazione, y = Peso)) +
  geom_point(alpha = 0.3, color = "steelblue") +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Effetto della gestazione sul peso",
       x = "Settimane di gestazione",
       y = "Peso (g)") +
  theme_minimal()

#3 EFFETTO DEL FUMO SUL PESO DEL NEUNATO
ggplot(dati_clean, aes(x = Fumatrici, y = Peso, 
                        fill = Fumatrici)) +
  geom_boxplot() +
  labs(title = "Effetto del fumo sul peso",
       x = "Fumatrici",
       y = "Peso (g)") +
  scale_x_discrete(labels = c("0" = "Non fumatrice",
                               "1" = "Fumatrice")) +
  theme_minimal()

#4 Effetto del sesso M o F sul peso
ggplot(dati_clean, aes(x = Sesso, y = Peso, 
                        fill = Sesso)) +
  geom_boxplot() +
  labs(title = "Effetto del sesso sul peso",
       x = "Sesso",
       y = "Peso (g)") +
  theme_minimal()

PROF_AI_Statistica_inferenziale_GL

Giovanni Lazzaro

2026-03-26