1 Introduzione

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:

1.Miglioramento delle previsioni cliniche:

2.Ottimizzazione delle risorse ospedaliere:

3.Prevenzione e identificazione dei fattori di rischio:

4.Valutazione delle pratiche ospedaliere:

5.Supporto alla pianificazione strategica:

Librerie 
library(ggplot2)
library(moments)
library(corrplot)
library(car)
library(MASS)
library(stats)
library(dplyr)
library(lmtest)
library(knitr)
Importazione dei dati 
d <- read.csv("neonati.csv",
              stringsAsFactors = FALSE,
              fileEncoding = "UTF-8")

Controllo che il dataset sia stato caricato correttamente

head(d)
##   Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1         26            0         0         42 3380       490    325        Nat
## 2         21            2         0         39 3150       490    345        Nat
## 3         34            3         0         38 3640       500    375        Nat
## 4         28            1         0         41 3690       515    365        Nat
## 5         20            0         0         38 3700       480    335        Nat
## 6         32            0         0         40 3200       495    340        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
## 6     osp2     F

Verifico che il dataset non presenti valori mancanti

data_na <- sum(is.na(d))

Il dataseta presenta 0 valori mancanti

Converto i dati in fattori e variabili numeriche a seconda che si tratti di variabili qualitative o quantitative

to_factor <- function(x) if(!is.null(x)) d[[x]] <<- factor(d[[x]])
cols <- c("Tipo.parto", "Ospedale", "Sesso")
d[cols] <- lapply(d[cols], factor)
d$Fumatrici <- factor(ifelse(d$Fumatrici == 1, "si", "no"))

cols <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")
d[cols] <- lapply(d[cols], as.numeric)

2 Analisi descrittiva

Procediamo con un’analisi descrittiva delle variabili quantitative:

desc <- data.frame(
  media     = sapply(d[cols], mean, na.rm = TRUE),
  median    = sapply(d[cols], median, na.rm = TRUE),
  min       = sapply(d[cols], min, na.rm = TRUE),
  max       = sapply(d[cols], max, na.rm = TRUE),
  qtl.25    = sapply(d[cols], quantile, probs=0.25,na.rm = TRUE),
  qtl.75    = sapply(d[cols], quantile, probs=0.75,na.rm = TRUE),
  sd        = sapply(d[cols], sd, na.rm = TRUE),
  skewness  = sapply(d[cols], skewness, na.rm = TRUE),
  kurtosi   = sapply(d[cols], function(x) kurtosis(x, na.rm=TRUE)-3)
)
kable(desc)
media median min max qtl.25 qtl.75 sd skewness kurtosi
Anni.madre 28.1640 28 0 46 25 32 5.273578 0.0428115 0.3804165
N.gravidanze 0.9812 1 0 12 0 1 1.280587 2.5142541 10.9894064
Gestazione 38.9804 39 25 43 38 40 1.868639 -2.0653133 8.2581496
Peso 3284.0808 3300 830 4930 2990 3620 525.038744 -0.6470308 2.0315318
Lunghezza 494.6920 500 310 565 480 510 26.318644 -1.5146991 6.4871739
Cranio 340.0292 340 235 390 330 350 16.425330 -0.7850527 2.9462063

Il valore minimo della variabile Anni.madre suggerisce che all’interno del dataset potrebbero esserci dei valori anomali o inseriti in maniera errata. Si procede quindi con l’identificazione e gestione degli outliers.

2.1 Gestione degli outliers

Implementiamo una funzione che identifica i valori oltre 1.5×IQR dal 1° e 3° quartile

find_outliers <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower <- Q1 - 1.5 * IQR
  upper <- Q3 + 1.5 * IQR
  which(x < lower | x > upper)
}

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

outlier_indices <- lapply(d[cols], find_outliers)
names(outlier_indices) <- cols

outliers_all <- unique(unlist(outlier_indices))
head(d[outliers_all, ]) # Visualizzo degli esempi delle righe che contengono almeno un outlier
##     Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 138         13            0        no         38 2760       470    325
## 205         45            2        no         38 3850       505    384
## 230         43            1        no         35 2050       455    305
## 260         44            1        no         40 3500       480    346
## 335         44            0        si         38 3150       465    335
## 855         43            0        no         38 3600       510    336
##     Tipo.parto Ospedale Sesso
## 138        Nat     osp2     F
## 205        Nat     osp3     M
## 230        Nat     osp2     F
## 260        Nat     osp1     F
## 335        Nat     osp3     F
## 855        Nat     osp3     M
d$outlier_flag <- seq_len(nrow(d)) %in% unique(unlist(outlier_indices))
table(d$outlier_flag) # numero totale di outliers
## 
## FALSE  TRUE 
##  2144   356
sapply(outlier_indices, length) # outliers divisi per colonna (N.B.: potrebbero esserci dei valori incrociati)
##   Anni.madre N.gravidanze   Gestazione         Peso    Lunghezza       Cranio 
##           13          246           67           69           59           48

L’algoritmo ha identificato 356 outliers. Tuttavia, operando in ambito sanitario e grazie alla tabella degli indici statistici, si è deciso di rimuovere esclusivamente gli outliers associati alla variabile Anni.madre e di mantenere i restanti in quanto potrebbero essere legati a casi clinici rari che necessitano di un’analisi dedicata.

N.B.: La variabile Numero di gravidanze presenta unità statistiche con valore 0 pertanto è necessario indagare sul suo significato. Si è deciso di mantenerne i valori anomali nel dataset in quanto è stata interpretata come numero di gravidanze precedenti al parto corrente.

cols <- c("Anni.madre")

remove_outliers <- function(df, col) {
  Q1 <- quantile(df[[col]], 0.25, na.rm = TRUE)
  Q3 <- quantile(df[[col]], 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower <- Q1 - 1.5 * IQR
  upper <- Q3 + 1.5 * IQR
  df[df[[col]] >= lower & df[[col]] <= upper, ]
}

for (c in cols) {
  d <- remove_outliers(d, c)
}

Ripeto l’analisi descrittiva sul dataset pulito

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

desc_new <- data.frame(
  media     = sapply(d[cols], mean, na.rm = TRUE),
  median    = sapply(d[cols], median, na.rm = TRUE),
  min       = sapply(d[cols], min, na.rm = TRUE),
  max       = sapply(d[cols], max, na.rm = TRUE),
  qtl.25    = sapply(d[cols], quantile, probs=0.25,na.rm = TRUE),
  qtl.75    = sapply(d[cols], quantile, probs=0.75,na.rm = TRUE),
  sd        = sapply(d[cols], sd, na.rm = TRUE),
  skewness  = sapply(d[cols], skewness, na.rm = TRUE),
  kurtosi   = sapply(d[cols], function(x) kurtosis(x, na.rm=TRUE)-3)
)
kable(desc_new)
media median min max qtl.25 qtl.75 sd skewness kurtosi
Anni.madre 28.1523924 28 15 42 25 32 5.124800 0.1095618 -0.2605367
N.gravidanze 0.9798955 1 0 12 0 1 1.277788 2.5229692 11.1167651
Gestazione 38.9835143 39 25 43 38 40 1.869681 -2.0724822 8.2940113
Peso 3284.6328910 3300 830 4930 2990 3620 525.323328 -0.6470005 2.0357109
Lunghezza 494.7426618 500 310 565 480 510 26.350546 -1.5201400 6.4988027
Cranio 340.0096502 340 235 390 330 350 16.411686 -0.7928552 2.9614758

  • L’età delle madri presenta una distribuzione pressoché simmetrica (skewness prossima a zero) e una moderata dispersione, come indicato dalla deviazione standard.

  • Il numero medio di gravidanze mostra una lieve asimmetria positiva, segnalando la presenza di alcune madri con un numero di gravidanze superiore alla media.

  • La durata della gestazione presenta una media intorno alle 39 settimane, con una bassa variabilità: ciò suggerisce che la maggior parte dei parti avviene a termine. I valori minimi e massimi permettono di identificare i casi di nascita prematura e post-termine.

  • Il peso dei neonati ha una distribuzione leggermente asimmetrica a sinistra, tipica di variabili biologiche soggette a limiti inferiori fisiologici. La deviazione standard indica una moderata dispersione, coerente con la variabilità osservata nei neonati sani.

  • Le misure antropometriche (lunghezza e diametro cranico) mostrano valori medi compatibili con la letteratura clinica. La bassa kurtosi e skewness indicano una distribuzione vicina alla normalità, confermando che i dati non presentano valori estremi o anomali evidenti.

2.2 Variabili qualitative

Valuto le frequenze relative di ciascuna variabile qualitativa

perc <- function(var) {
  tab <- prop.table(table(var)) * 100
  round(tab, 2)
}

# Distribuzione per Ospedale
kable(perc(d$Ospedale))
var Freq
osp1 32.65
osp2 33.94
osp3 33.41 
# Distribuzione per Sesso
kable(perc(d$Sesso))
var Freq
F 50.26
M 49.74 
# Distribuzione per Fumo
kable(perc(d$Fumatrici))
var Freq
no 95.86
si 4.14 
# Distribuzione per tipo di parto
kable(perc(d$Tipo.parto))
var Freq
Ces 29.19
Nat 70.81

2.3 Grafici esplorativi

ggplot(d, aes(x = Peso)) +
  geom_histogram(aes(y=..density..), bins = 30) +
  geom_density(color="red") +
  labs(title="Distribuzione del peso alla nascita", x="Peso (grammi)", y="Densità")

Il grafico mostra la distribuzione del peso dei neonati alla nascita.
L’istogramma rappresenta la frequenza relativa dei valori di peso, mentre la curva rossa sovrapposta indica la densità stimata.

Osservando la forma della distribuzione, si nota che essa è approssimativamente simmetrica, con un picco intorno ai 3200–3300 grammi, valore coerente con la media calcolata in precedenza.
La coda destra leggermente più lunga suggerisce una lieve asimmetria negativa (skewness < 0), coerente con la presenza di alcuni neonati con peso inferiore alla media, probabilmente dovuto a nascite premature o a condizioni materne particolari.

ggplot(d, aes(x = Sesso, y = Peso)) +
  geom_boxplot() +
  labs(title="Peso per sesso", x="Sesso", y="Peso (grammi)")

Il boxplot mostra la distribuzione del peso alla nascita suddivisa per sesso. Osservando il grafico, si nota che i neonati di sesso maschile (M) presentano in media un peso leggermente superiore rispetto alle neonate (F).
La mediana e la posizione della scatola risultano infatti spostate verso valori più alti nei maschi, confermando una differenza di peso tipica e coerente con la letteratura clinica.

La variabilità interna (indicata dall’ampiezza del box) appare simile tra i due gruppi, ma si osserva la presenza di outlier in entrambe le categorie.

ggplot(d, aes(x = Ospedale, fill = Tipo.parto)) +
  geom_bar(position="fill") +
  labs(title="Composizione % tipo di parto per ospedale",
       x="Ospedale", y="Proporzione") +
  guides(fill = guide_legend(title="Tipo di parto"))

Il grafico mostra la distribuzione percentuale dei tipi di parto (naturale e cesareo) nei tre ospedali analizzati. Osservando il grafico, si nota che la composizione è piuttosto simile tra le tre strutture, con una prevalenza di parti naturali in tutti gli ospedali.

Queste ed altre considerazioni verranno saggiate nella sezione successiva tramite opportuni test statistici.

2.4 Test statistici

2.4.1 Tipo di Parto x Ospedale

L’obiettivo è verificare se c’è associazione tra le due variabili categoriche oppure, in altri termini, esaminare se la distribuzione del tipo di parto (naturale o cesareo) cambia al variare dell’ospedale. Per condurre quest’analisi il test statistico più indicato è il test χ² di indipendenza

Ipotesi

  • \(H_0\): tipo di parto e ospedale sono indipendenti (stesse proporzioni di cesarei in tutti gli ospedali).

  • \(H_1\): esiste associazione (almeno un ospedale ha proporzione diversa).

tab_parto_osp <- table(Ospedale = d$Ospedale, Parto = d$Tipo.parto)
kable(tab_parto_osp)
Ces Nat
osp1 240 572
osp2 254 590
osp3 232 599 
test_chi <- chisq.test(tab_parto_osp, correct = FALSE)
print(test_chi)
## 
##  Pearson's Chi-squared test
## 
## data:  tab_parto_osp
## X-squared = 1.0374, df = 2, p-value = 0.5953

Siccome p-value = 0,5953 > 0,05 non rifiuto \(H_0\): i dati non forniscono evidenza di differenze nelle percentuali di parti cesarei tra i tre ospedali.

Concludiamo che il campione non fornisce alcuna differenza statisticamente rilevante nella distribuzione del tipo di parto tra i tre ospedali.

2.4.2 Media del peso e della lunghezza VS popolazione

L’obiettivo è quello di confrontare la media di due variabili quantitative continue (peso e lunghezza) di un campione con dei valori noti di riferimento nel caso in cui la varianza della popolazione non è nota. Si è scelto pertanto di condurre il t-test per un campione per ciascuna variabile.

Ipotesi

  • \(H_0\): \(\mu = \mu_0\) il campione rappresenta la popolazione di riferimento

  • \(H_1\): \(\mu \not= \mu_0\) la media del campione è diversa da quella di riferimento

I valori medi della popolazione sono 3300 g e 500 mm (fonte: Ospedale Bambino Gesù)

mu_peso_pop <- 3300
mu_lung_pop <- 500      

t_peso <- t.test(d$Peso, mu = mu_peso_pop)
t_lung <- t.test(d$Lunghezza, mu = mu_lung_pop)

print(t_peso)
## 
##  One Sample t-test
## 
## data:  d$Peso
## t = -1.4588, df = 2486, p-value = 0.1447
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.977 3305.289
## sample estimates:
## mean of x 
##  3284.633
print(t_lung)
## 
##  One Sample t-test
## 
## data:  d$Lunghezza
## t = -9.9498, df = 2486, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.7065 495.7788
## sample estimates:
## mean of x 
##  494.7427

Peso

Siccome p-value = 0,1447 > 0,05 non rifiuto \(H_0\): il peso medio dei neonati nel campione (3284,6 g) non differisce in modo significativo dal valore medio della popolazione (3300 g) quindi il campione può essere considerato rappresentativo della popolazione per quanto riguarda il peso.

Lunghezza

Siccome p-value = 2,2e-16 < 0,05 rifiuto \(H_0\): la lunghezza media dei neonati nel campione (494,7 mm) è significativamente inferiore rispetto alla popolazione di riferimento (500 mm, circa 5 mm in meno). Siccome si tratta di una differenza statisticamente significativa ma piccola in termini pratici si è ritenuto opportuno considerare questo risultato marginale nelle valutazioni successive.

2.4.3 Differenze tra sessi nelle misure antropometriche

L’obiettivo è quello di verificare se le misure antropometriche (descritte dalle variabili quantitative continue Peso, Lunghezza e Cranio) sono significativamente diverse tra i due sessi. Poichè la variabile di confronto ha due categorie e i campioni sono indipendenti (ogni neonato appartiene ad un solo gruppo, maschi o femmine) si è scelto di utilizzare il t-test per campioni indipendenti

Ipotesi

  • \(H_0\): \(\mu_\text{maschi} = \mu_\text{femmine}\) nessuna differenza di media tra i sessi

  • \(H_1\): \(\mu_\text{maschi} \not= \mu_\text{femmine}\) esiste una differenza di media tra i sessi

t_peso_sesso <- t.test(Peso ~ Sesso, data = d)
t_lung_sesso <- t.test(Lunghezza ~ Sesso, data = d)
t_cranio_sesso <- t.test(Cranio ~ Sesso, data = d)

print(t_peso_sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.05, df = 2478.2, 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:
##  -286.8803 -206.5763
## sample estimates:
## mean in group F mean in group M 
##        3161.914        3408.642
print(t_lung_sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.4851, df = 2447.5, 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:
##  -11.877170  -7.807573
## sample estimates:
## mean in group F mean in group M 
##        489.8472        499.6896
print(t_cranio_sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.2915, df = 2478.2, p-value = 4.102e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.025036 -3.471203
## sample estimates:
## mean in group F mean in group M 
##        337.6480        342.3961

Siccome per tutte le misure antropometrice vale p-value < 0,05 rifiuto \(H_0\): tutte le differenze tra i sessi sono statisticamente significative. In media, i neonati maschi presantano valori maggiori in tutte le misure antropometriche considerate: peso, lunghezza e diametro craniale.

3 Modello di regressione lineare multipla

Prima di procedere con la costruzione del modello di regressione lineare multipla, è fondamentale esplorare le relazioni tra le variabili quantitative presenti nel dataset.

A tal fine, è stata realizzata una matrice di dispersione che consente di osservare contemporaneamente:

Questa analisi preliminare è utile per:

d <- d[,1:10] #rimuovo eventuali colonne create precedentemente

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  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)
}

pairs(d, upper.panel = panel.smooth, lower.panel = panel.cor)

Dalla matrice di dispersione si osservano alcune relazioni interessanti:

Si evidenzia quindi un insieme di forti correlazioni tra variabili antropometriche, che sarà importante considerare nel modello di regressione per evitare ridondanze o collinearità.

Si procede con la costruzione di un modello di regressione lineare multipla con l’obiettivo di stimare il peso alla nascita del neonato.

Tale approccio consente di quantificare l’effetto simultaneo di più fattori sulla variabile risposta, tenendo conto delle possibili relazioni e interazioni tra di essi.

In particolare, il modello adottato segue la forma generale:

\[ \textit{Peso} = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2} + \dots + \beta_k X_{k} + \varepsilon \]

dove:

mod1 <- lm(Peso ~ ., data = d)

summary(mod1)
## 
## Call:
## lm(formula = Peso ~ ., data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1125.81  -181.68   -15.02   161.00  2613.60 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6732.0237   141.8849 -47.447  < 2e-16 ***
## Anni.madre        0.8754     1.1708   0.748   0.4547    
## N.gravidanze     11.5221     4.6946   2.454   0.0142 *  
## Fumatricisi     -33.3922    27.6949  -1.206   0.2280    
## Gestazione       32.3695     3.8274   8.457  < 2e-16 ***
## Lunghezza        10.3003     0.3014  34.173  < 2e-16 ***
## Cranio           10.4639     0.4281  24.443  < 2e-16 ***
## Tipo.partoNat    30.3139    12.1145   2.502   0.0124 *  
## Ospedaleosp2    -10.1785    13.4854  -0.755   0.4505    
## Ospedaleosp3     27.7962    13.5394   2.053   0.0402 *  
## SessoM           78.5901    11.2117   7.010 3.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.1 on 2476 degrees of freedom
## Multiple R-squared:  0.7288, Adjusted R-squared:  0.7277 
## F-statistic: 665.4 on 10 and 2476 DF,  p-value: < 2.2e-16

Verifichiamo ora le principali assunzioni teoriche del modello lineare classico, al fine di garantire la validità delle inferenze statistiche e la correttezza delle stime dei coefficienti.

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

shapiro.test(residuals(mod1))  # normalità dei residui
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod1)
## W = 0.97397, p-value < 2.2e-16
lmtest::bptest(mod1)           # omoschedasticità
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 93.666, df = 10, p-value = 1.001e-15
lmtest::dwtest(mod1)           # indipendenza degli errori
## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 1.9583, p-value = 0.1494
## alternative hypothesis: true autocorrelation is greater than 0
car::vif(mod1)                 # multicollinearità 
##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.191168  1        1.091406
## N.gravidanze 1.190548  1        1.091122
## Fumatrici    1.007852  1        1.003918
## Gestazione   1.694200  1        1.301614
## Lunghezza    2.087152  1        1.444698
## Cranio       1.633153  1        1.277949
## Tipo.parto   1.004065  1        1.002030
## Ospedale     1.004509  2        1.001125
## Sesso        1.040108  1        1.019857

Nel complesso, le assunzioni fondamentali del modello lineare risultano soddisfatte in misura accettabile. Le lievi violazioni della normalità e dell’omoschedasticità dei residui non compromettono la validità inferenziale del modello, ma suggeriscono che modelli alternativi potrebbero migliorare ulteriormente l’adattamento ai dati.

Si procede costruendo un modello completo che includa anche termini di secondo ordine e interazioni tra variabili suggerite dall’osservazione della matrice di correlazione.

L’obiettivo è verificare se la relazione tra il peso neonatale e le variabili esplicative segua effettivamente un andamento lineare o se siano presenti effetti non lineari o combinati che migliorano la capacità predittiva del modello.

d <- d[,1:10]

mod_full <- lm(Peso ~ 
               Gestazione + Anni.madre + N.gravidanze + Lunghezza + Cranio 
               + Fumatrici + Sesso + Tipo.parto 
               + I(Gestazione^2) + I(Anni.madre^2) + I(N.gravidanze^2) 
               + I(Lunghezza^2) + I(Cranio^2) 
               + Gestazione:Lunghezza + Gestazione:Cranio
               + Lunghezza:Cranio + Gestazione:Fumatrici
               + Gestazione:Sesso, 
               data=d)

summary(mod_full)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Anni.madre + N.gravidanze + 
##     Lunghezza + Cranio + Fumatrici + Sesso + Tipo.parto + I(Gestazione^2) + 
##     I(Anni.madre^2) + I(N.gravidanze^2) + I(Lunghezza^2) + I(Cranio^2) + 
##     Gestazione:Lunghezza + Gestazione:Cranio + Lunghezza:Cranio + 
##     Gestazione:Fumatrici + Gestazione:Sesso, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1137.91  -178.44   -11.94   163.32  1325.39 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -158.13507 1210.25632  -0.131  0.89605    
## Gestazione              152.71131   76.07563   2.007  0.04482 *  
## Anni.madre               11.77320    8.91444   1.321  0.18673    
## N.gravidanze             33.27562    8.32437   3.997 6.59e-05 ***
## Lunghezza                -3.21543    5.93314  -0.542  0.58791    
## Cranio                  -24.44760   10.23996  -2.387  0.01704 *  
## Fumatricisi             722.47916  742.40880   0.973  0.33057    
## SessoM                 -145.38883  234.17649  -0.621  0.53475    
## Tipo.partoNat            30.37478   11.78236   2.578  0.01000 ** 
## I(Gestazione^2)          -4.48544    1.63072  -2.751  0.00599 ** 
## I(Anni.madre^2)          -0.20610    0.15619  -1.320  0.18710    
## I(N.gravidanze^2)        -3.34706    1.28674  -2.601  0.00935 ** 
## I(Lunghezza^2)            0.05596    0.00651   8.596  < 2e-16 ***
## I(Cranio^2)               0.05100    0.01886   2.705  0.00689 ** 
## Gestazione:Lunghezza     -0.29341    0.19719  -1.488  0.13688    
## Gestazione:Cranio         1.10252    0.20933   5.267 1.51e-07 ***
## Lunghezza:Cranio         -0.08682    0.01672  -5.193 2.23e-07 ***
## Gestazione:Fumatricisi  -19.21502   18.89188  -1.017  0.30920    
## Gestazione:SessoM         5.66238    5.99860   0.944  0.34529    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266 on 2468 degrees of freedom
## Multiple R-squared:  0.7455, Adjusted R-squared:  0.7436 
## F-statistic: 401.5 on 18 and 2468 DF,  p-value: < 2.2e-16

Il modello presenta un \(R^2_\text{adj}=\) 0.744 indicando un ottimo livello di adattamento ai dati.

Per identificare il modello più parsimonioso applichiamo due criteri di selezione automatica:

mod_aic <- stepAIC(mod_full, direction = "both", k = 2,  trace =F)

summary(mod_aic)
## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + Tipo.parto + I(Gestazione^2) + I(N.gravidanze^2) + 
##     I(Lunghezza^2) + I(Cranio^2) + Gestazione:Lunghezza + Gestazione:Cranio + 
##     Lunghezza:Cranio, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1171.90  -177.21   -11.02   163.71  1318.86 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.989e+01  1.195e+03   0.042  0.96671    
## Gestazione            1.486e+02  7.529e+01   1.973  0.04856 *  
## N.gravidanze          3.249e+01  7.876e+00   4.126 3.82e-05 ***
## Lunghezza            -3.537e+00  5.924e+00  -0.597  0.55048    
## Cranio               -2.408e+01  1.023e+01  -2.354  0.01867 *  
## SessoM                7.421e+01  1.093e+01   6.787 1.43e-11 ***
## Tipo.partoNat         3.016e+01  1.177e+01   2.564  0.01042 *  
## I(Gestazione^2)      -4.492e+00  1.613e+00  -2.784  0.00541 ** 
## I(N.gravidanze^2)    -3.311e+00  1.271e+00  -2.605  0.00925 ** 
## I(Lunghezza^2)        5.604e-02  6.505e-03   8.614  < 2e-16 ***
## I(Cranio^2)           5.083e-02  1.884e-02   2.698  0.00703 ** 
## Gestazione:Lunghezza -2.811e-01  1.969e-01  -1.428  0.15345    
## Gestazione:Cranio     1.104e+00  2.091e-01   5.282 1.39e-07 ***
## Lunghezza:Cranio     -8.742e-02  1.671e-02  -5.231 1.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266 on 2473 degrees of freedom
## Multiple R-squared:  0.7449, Adjusted R-squared:  0.7436 
## F-statistic: 555.5 on 13 and 2473 DF,  p-value: < 2.2e-16
n <- nrow(d)
mod_bic <- stepAIC(mod_full, direction = "both", k = log(n), trace =F)

summary(mod_bic)
## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + 
##     Gestazione:Cranio + Lunghezza:Cranio, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1168.08  -182.32   -11.74   168.64  1319.67 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.651e+02  1.193e+03  -0.138  0.88993    
## Gestazione         1.787e+02  7.449e+01   2.399  0.01652 *  
## N.gravidanze       1.497e+01  4.237e+00   3.532  0.00042 ***
## Lunghezza         -6.958e+00  5.669e+00  -1.227  0.21982    
## Cranio            -2.102e+01  1.012e+01  -2.077  0.03790 *  
## SessoM             7.366e+01  1.095e+01   6.724 2.19e-11 ***
## I(Gestazione^2)   -6.272e+00  1.034e+00  -6.068 1.49e-09 ***
## I(Lunghezza^2)     5.028e-02  5.049e-03   9.960  < 2e-16 ***
## I(Cranio^2)        5.575e-02  1.854e-02   3.006  0.00267 ** 
## Gestazione:Cranio  1.012e+00  2.062e-01   4.907 9.86e-07 ***
## Lunghezza:Cranio  -9.292e-02  1.581e-02  -5.877 4.75e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.6 on 2476 degrees of freedom
## Multiple R-squared:  0.7434, Adjusted R-squared:  0.7424 
## F-statistic: 717.3 on 10 and 2476 DF,  p-value: < 2.2e-16
AIC(mod_full, mod_aic, mod_bic)
##          df      AIC
## mod_full 20 34851.12
## mod_aic  15 34846.30
## mod_bic  12 34855.20
BIC(mod_full, mod_aic, mod_bic)
##          df      BIC
## mod_full 20 34967.50
## mod_aic  15 34933.58
## mod_bic  12 34925.03

Il confronto tra i criteri indica che:

mod_star <- mod_bic

In conclusione, il modello selezionato secondo BIC rappresenta un compromesso ottimale tra accuratezza statistica e semplicità interpretativa.

È interessante osservare che nel modello selezionato da BIC compaiono solo 5 dei 9 possibili regressori, in particolare sono state rimosse le variabili Anni.madre e Fumatrici suggerendo che il loro impatto sul peso neonatale è molto meno significativo di quanto si potrebbe ipotizzare.

3.1 Analisi della qualità del modello

Sottoponiamo il modello così selezionato alla verifica delle assunzioni a cui abbiamo sottoposto i modelli precedenti con l’obiettivo di confermare la bontà delle stime e l’adeguatezza del modello finale in termini inferenziali e predittivi.

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

I grafici mostrano:

  • una distribuzione dei residui centrata attorno a zero, senza pattern evidenti che indichino relazioni non lineari;

  • una varianza dei residui quasi costante, con un leggero aumento nella parte destra del grafico;

  • il Q–Q plot mostra solo piccole deviazioni nelle code, indice di una buona approssimazione alla normalità;

  • nessun punto di leverage particolarmente influente, confermando l’assenza di osservazioni anomale che distorcano il modello.

shapiro.test(mod_star$residuals)  
## 
##  Shapiro-Wilk normality test
## 
## data:  mod_star$residuals
## W = 0.99095, p-value = 2.181e-11
plot(density(residuals(mod_star)))

lmtest::bptest(mod_star)           
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_star
## BP = 41.138, df = 10, p-value = 1.067e-05
lmtest::dwtest(mod_star)           
## 
##  Durbin-Watson test
## 
## data:  mod_star
## DW = 1.9644, p-value = 0.187
## alternative hypothesis: true autocorrelation is greater than 0
car::vif(mod_star, type="predictor")
## GVIFs computed for predictors
##                     GVIF Df GVIF^(1/(2*Df))                     Interacts With
## Gestazione   1681.088593  5        2.101644            I(Gestazione^2), Cranio
## N.gravidanze    1.024914  1        1.012380                               --  
## Lunghezza    1735.975208  5        2.108407             I(Lunghezza^2), Cranio
## Cranio          1.073263  8        1.004429 I(Cranio^2), Gestazione, Lunghezza
## Sesso           1.049274  1        1.024341                               --  
##                                         Other Predictors
## Gestazione                N.gravidanze, Lunghezza, Sesso
## N.gravidanze        Gestazione, Lunghezza, Cranio, Sesso
## Lunghezza                Gestazione, N.gravidanze, Sesso
## Cranio                               N.gravidanze, Sesso
## Sesso        Gestazione, N.gravidanze, Lunghezza, Cranio

  • Il test di Shapiro–Wilk (\(W =\) 0.9909533, \(\textit{p-value} =\) 2.1807516^{-11}) evidenzia una lieve deviazione dalla normalità teorica, ma il grafico di densità mostra una distribuzione simmetrica e centrata attorno a zero. Considerata la numerosità campionaria elevata (\(n \approx 2500\)), tale deviazione non compromette la validità del modello.

  • Il test di Breusch–Pagan (\(BP =\) 41.1379173, \(\textit{p-value} =\) 1.0667069^{-5}) segnala una lieve eterogeneità della varianza, ma molto ridotta rispetto ai modelli precedenti. La distribuzione dei residui nei grafici Residuals vs Fitted e Scale–Location appare più omogenea e stabile.

  • Il test di Durbin–Watson (\(DW =\) 1.9643581, \(\textit{p-value} =\) 0.1870408) non fornisce evidenze di autocorrelazione, quindi l’ipotesi di indipendenza può ritenersi soddisfatta.

  • I valori del GVIF, corretti per i gradi di libertà, sono tutti prossimi a 1, ben al di sotto della soglia critica di 5. Ciò conferma che le variabili incluse non risultano eccessivamente correlate tra loro.

Per completare la valutazione del modello finale, si conduce un’analisi dei punti influenti e degli outlier, essenziale per verificare che la stima dei coefficienti non sia eccessivamente condizionata da poche osservazioni anomale.

lev <- hatvalues(mod_star)
plot(lev, main = "Leverage dei punti", ylab = "Leverage")
abline(h = 2*mean(lev), col = "red")

Il grafico dei valori di leverage mostra la maggior parte delle osservazioni al di sotto della soglia 2h (linea rossa). Solo pochi punti si distinguono come potenzialmente influenti, ma la loro incidenza sul modello risulta marginale. In generale, non emergono casi con leverage eccessivo tali da richiedere esclusione.

stud_res <- rstudent(mod_star)
plot(stud_res, main = "Residuals Studentizzati", ylab = "Valore")
abline(h = c(-2,2), col="red", lty=2)

outlierTest(mod_star)
##       rstudent unadjusted p-value Bonferroni p
## 1306  4.978656         6.8434e-07    0.0017019
## 155   4.604442         4.3449e-06    0.0108060
## 1399 -4.403448         1.1107e-05    0.0276230
## 1694  4.396472         1.1467e-05    0.0285190

La maggior parte dei residui studentizzati si distribuisce nell’intervallo [−2,2], mentre solo un numero limitato di punti supera le soglie critiche. Questo suggerisce che la quasi totalità delle osservazioni è ben rappresentata dal modello, con pochi casi anomali isolati.

L’analisi numerica tramite outlierTest() conferma la presenza di quattro osservazioni potenzialmente influenti che, pur risultando statisticamente rilevanti dal test di Bonferroni, il loro effetto complessivo è contenuto, come confermato dalle metriche di influenza.

cook <- cooks.distance(mod_star)
plot(cook, main = "Distanza di Cook", ylab = "Cook's distance")
abline(h = 4/length(cook), col="red", lty=2)

Il grafico della distanza di Cook mostra che la quasi totalità delle osservazioni presenta valori inferiori alla soglia di riferimento. Solo poche unità (coincidenti con quelle già identificate come outlier) presentano valori leggermente più alti e solo una ha una distanza di Cook maggiore di 1, valore che in letteratura è indice di un osservazione altamente influente da indagare attentamente.

out <- outlierTest(mod_star)
outlier_ids <- names(out$rstudent)
outliers_df <- d[outlier_ids, ]
print(outliers_df)
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1306         23            0        no         41 4900       510    352
## 155          30            0        no         36 3610       410    330
## 1399         42            2        no         38 2560       525    349
## 1694         23            1        no         36 3850       460    334
##      Tipo.parto Ospedale Sesso
## 1306        Nat     osp2     F
## 155         Nat     osp1     M
## 1399        Ces     osp2     M
## 1694        Ces     osp3     F

Osservando gli outliers possiamo concludere che non si tratta di veri e propri valori anomali bensì, trovandosi in un contesto biologico, di eventi clinici da segnalare come casi speciali e in quanto tali possiamo scegliere di mantenerli all’interno del dataset.

Modello senza outliers

Per dimostrare questa scelta costruiamo un nuovo dataset rimuovendo gli outliers e ripetiamo il procedimento di creazione del modello e di selezione delle variabili con BIC.

outlier_ids <- as.numeric(names(out$rstudent))

d_no_out <- d[-outlier_ids, ]

mod_out <- lm(Peso ~ 
               Gestazione + Anni.madre + N.gravidanze + Lunghezza + Cranio 
               + Fumatrici + Sesso + Tipo.parto 
               + I(Gestazione^2) + I(Anni.madre^2) + I(N.gravidanze^2) 
               + I(Lunghezza^2) + I(Cranio^2) 
               + Gestazione:Lunghezza + Gestazione:Cranio
               + Lunghezza:Cranio + Gestazione:Fumatrici
               + Gestazione:Sesso, 
               data=d_no_out)

summary(mod_out)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Anni.madre + N.gravidanze + 
##     Lunghezza + Cranio + Fumatrici + Sesso + Tipo.parto + I(Gestazione^2) + 
##     I(Anni.madre^2) + I(N.gravidanze^2) + I(Lunghezza^2) + I(Cranio^2) + 
##     Gestazione:Lunghezza + Gestazione:Cranio + Lunghezza:Cranio + 
##     Gestazione:Fumatrici + Gestazione:Sesso, data = d_no_out)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1137.25  -178.12   -11.91   163.15  1326.01 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -163.35867 1210.24598  -0.135  0.89264    
## Gestazione              155.30684   76.10172   2.041  0.04138 *  
## Anni.madre               11.52420    8.91578   1.293  0.19628    
## N.gravidanze             33.22873    8.32635   3.991 6.78e-05 ***
## Lunghezza                -3.26337    5.93451  -0.550  0.58244    
## Cranio                  -24.58858   10.24091  -2.401  0.01642 *  
## Fumatricisi             716.30393  742.40228   0.965  0.33472    
## SessoM                 -152.44417  234.29200  -0.651  0.51533    
## Tipo.partoNat            30.26116   11.79745   2.565  0.01037 *  
## I(Gestazione^2)          -4.62552    1.63538  -2.828  0.00472 ** 
## I(Anni.madre^2)          -0.20239    0.15621  -1.296  0.19522    
## I(N.gravidanze^2)        -3.33284    1.28688  -2.590  0.00966 ** 
## I(Lunghezza^2)            0.05543    0.00652   8.502  < 2e-16 ***
## I(Cranio^2)               0.05188    0.01887   2.750  0.00600 ** 
## Gestazione:Lunghezza     -0.27375    0.19787  -1.383  0.16664    
## Gestazione:Cranio         1.09757    0.20933   5.243 1.71e-07 ***
## Lunghezza:Cranio         -0.08737    0.01672  -5.225 1.89e-07 ***
## Gestazione:Fumatricisi  -19.05430   18.89175  -1.009  0.31326    
## Gestazione:SessoM         5.85013    6.00189   0.975  0.32980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266 on 2464 degrees of freedom
## Multiple R-squared:  0.7456, Adjusted R-squared:  0.7438 
## F-statistic: 401.2 on 18 and 2464 DF,  p-value: < 2.2e-16
n <- nrow(d_no_out)
mod_bic_out <- stepAIC(mod_out, direction = "both", k = log(n), trace =F)

summary(mod_bic_out)
## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + 
##     Gestazione:Cranio + Lunghezza:Cranio, data = d_no_out)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1168.31  -182.24   -11.35   168.51  1320.27 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.664e+02  1.193e+03  -0.140  0.88904    
## Gestazione         1.797e+02  7.449e+01   2.413  0.01590 *  
## N.gravidanze       1.494e+01  4.238e+00   3.526  0.00043 ***
## Lunghezza         -6.814e+00  5.669e+00  -1.202  0.22948    
## Cranio            -2.132e+01  1.012e+01  -2.107  0.03525 *  
## SessoM             7.394e+01  1.096e+01   6.745 1.90e-11 ***
## I(Gestazione^2)   -6.279e+00  1.034e+00  -6.075 1.43e-09 ***
## I(Lunghezza^2)     5.016e-02  5.049e-03   9.935  < 2e-16 ***
## I(Cranio^2)        5.630e-02  1.855e-02   3.036  0.00242 ** 
## Gestazione:Cranio  1.010e+00  2.062e-01   4.898 1.03e-06 ***
## Lunghezza:Cranio  -9.295e-02  1.581e-02  -5.880 4.66e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.6 on 2472 degrees of freedom
## Multiple R-squared:  0.7436, Adjusted R-squared:  0.7426 
## F-statistic: 716.9 on 10 and 2472 DF,  p-value: < 2.2e-16
Confrontando i risultati così ottenuti con i modelli precedenti concludiamo che la rimozione degli outliers non produce risultati significativi in quanto i coefficienti e le metriche sono rimasti pressocchè invariati.


Nel complesso:

  • i valori di leverage rientrano entro limiti accettabili;

  • i residui studentizzati evidenziano solo un numero limitato di outlier;

  • la distanza di Cook conferma l’assenza di osservazioni eccessivamente influenti.

Si può quindi concludere che il modello finale è stabile e robusto, in quanto le stime dei coefficienti non sono condizionate in modo rilevante da singole osservazioni anomale. Il comportamento dei residui e la distribuzione dei punti influenti confermano la buona qualità complessiva del modello selezionato secondo il criterio BIC.

3.2 Esempio di previsione

Stima del peso di una neonata con Gestazione e N.gravidanze note

newdata <- data.frame(
  Gestazione   = 39,
  Lunghezza    = mean(d$Lunghezza),
  Cranio       = mean(d$Cranio), 
  Anni.madre   = mean(d$Anni.madre),
  N.gravidanze = 2,
  Sesso        = factor("F",   levels = levels(d$Sesso))
)
pr <- predict(mod_star, newdata = newdata, se.fit = TRUE)
pr
## $fit
##        1 
## 3245.774 
## 
## $se.fit
## [1] 9.522485
## 
## $df
## [1] 2476
## 
## $residual.scale
## [1] 266.6487

3.3 Grafici

ggplot(data = d) +
  geom_point(aes(x = Gestazione,
                 y = Peso,
                 col = Sesso),
             position = "jitter", alpha = 0.2) +
  geom_smooth(aes(x = Gestazione,
                  y = Peso,
                  col = Sesso),
              se = FALSE, method = "lm", linewidth = 1.1) +
  labs(title = "Relazione tra settimane di gestazione e peso alla nascita",
       x = "Settimane di gestazione",
       y = "Peso del neonato (grammi)",
       color = "Sesso")+
  geom_smooth(aes(x=Gestazione, y=Peso), col="black", se=FALSE, method="lm",linewidth = 1) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

ggplot(data = d) +
  geom_point(aes(x = Gestazione, y = Peso, color = Sesso),
             position = "jitter", alpha = 0.2) +
  geom_smooth(aes(x = Gestazione, y = Peso, color = Sesso),
              se = F, method = "lm", linewidth = 1.1) +
  facet_wrap(~ Sesso) +
  labs(title = "Relazione tra gestazione e peso per sesso",
       x = "Settimane di gestazione",
       y = "Peso (g)",
       color = "Sesso") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Il grafico mostra una relazione fortemente positiva tra il numero di settimane di gestazione e il peso del neonato. Si osserva come, all’aumentare delle settimane di gestazione, il peso alla nascita cresca in modo quasi lineare. Le rette di regressione separate per sesso evidenziano che, a parità di settimane, i maschi tendono a pesare leggermente di più rispetto alle femmine, sebbene la differenza non sia marcata. La retta nera rappresenta la tendenza media generale del campione e conferma che la durata della gravidanza è una delle principali determinanti del peso neonatale, coerentemente con quanto atteso dal punto di vista biologico.

ggplot(data=d) +
  geom_boxplot(aes(x=Fumatrici, y=Peso, fill=Fumatrici)) +
  labs(title="Distribuzione del peso per fumatrici e non fumatrici",
       x="Fumo materno", y="Peso (g)") +
  theme_minimal()

Dal boxplot emerge che le madri fumatrici tendono ad avere neonati con peso medio leggermente inferiore rispetto alle non fumatrici. Tuttavia, la sovrapposizione delle due distribuzioni e la presenza di numerosi valori estremi (outlier) suggeriscono che l’effetto del fumo, pur potenzialmente presente, non è particolarmente marcato nel dataset considerato. Per questo motivo, oltre alla scarsa significativà del regressore, nella fase di modellizzazione si è deciso di non includere la variabile “Fumo” nel modello finale, poiché il suo contributo esplicativo al peso neonatale risultava statisticamente debole rispetto ad altre variabili.

ggplot(d, aes(x=Gestazione, y=Peso, col=Fumatrici)) +
  geom_point(position="jitter") +
  geom_smooth(aes(x=Gestazione, y=Peso, col=Fumatrici), method="lm", se=FALSE) +
  geom_smooth(aes(x=Gestazione, y=Peso), method="lm", se=FALSE, col="black", linewidth = 0.5) +
  labs(title="Interazione tra settimane di gestazione e fumo sul peso",
       x="Settimane di gestazione", y="Peso (g)")

Il grafico mostra l’interazione tra settimane di gestazione e fumo materno rispetto al peso del neonato. Le due rette (per fumatrici e non fumatrici) suggeriscono che i neonati di madri fumatrici hanno in media un peso lievemente inferiore rispetto alla controparte dei non fumatori. Tuttavia, essendo il dataset fortemente sbilanciato a favore delle non fumatrici, la relazione principale (indicata dalla retta nera) si confonda facilmente con quella delle non fumatrici. Questa vicinanza suggerisce che la variabile “Fumatrici” non modifica la pendenza della relazione principale e non aggiunge informazione rilevante al modello confermando la scelta metodologica di escludere il fumo come regressore nel modello finale.

dati <- d
dati$Peso_pred <- predict(mod_star)

ggplot(dati, aes(x=Peso, y=Peso_pred)) +
  geom_point(alpha=0.6) +
  geom_abline(intercept=0, slope=1, col="red") +
  labs(title="Peso osservato vs Peso previsto dal modello",
       x="Peso osservato (g)", y="Peso previsto (g)") +
  theme_minimal()

ggplot(dati, aes(x=fitted(mod_star), y=residuals(mod_star))) +
  geom_point(alpha=0.6) +
  geom_hline(yintercept=0, col="red") +
  labs(title="Analisi dei residui", x="Valori previsti", y="Residui") +
  theme_minimal()

Lo scatterplot confronta i valori di peso osservati nel dataset con quelli previsti dal modello finale. La disposizione dei punti lungo la bisettrice rossa indica un buon allineamento tra osservazioni e previsioni: il modello è in grado di riprodurre accuratamente la variabilità del peso all’interno del campione. Le lievi dispersioni intorno alla linea sono fisiologiche e rappresentano la componente di errore non spiegata dal modello, confermando una buona capacità predittiva complessiva.

Il grafico dei residui in funzione dei valori previsti permette di verificare la validità delle assunzioni del modello lineare. La distribuzione simmetrica e casuale dei residui intorno alla linea orizzontale rossa suggerisce che non vi sono evidenze di eteroschedasticità o di andamenti non lineari. Inoltre, l’assenza di pattern sistematici indica che le assunzioni di linearità e omoschedasticità sono rispettate, e che il modello stimato è statisticamente coerente con i dati osservati.

4 Conclusioni

Il progetto sviluppato per Neonatal Health Solutions ha permesso di costruire un modello statistico affidabile per la previsione del peso neonatale a partire da variabili cliniche materne e perinatali. L’analisi descrittiva preliminare ha fornito un quadro completo delle caratteristiche del campione, evidenziando differenze significative tra i sessi e tra i diversi ospedali.

Il modello di regressione lineare multipla ha confermato l’importanza delle settimane di gestazione come principale fattore predittivo del peso alla nascita, seguita da variabili come il numero di gravidanze e il sesso del neonato. La variabile relativa al fumo materno, pur mostrando una tendenza negativa sul peso, non è risultata statisticamente significativa e per questo non è stata inclusa nel modello finale.

Le metriche di valutazione, come \(R^2\) e RMSE, hanno evidenziato una buona capacità esplicativa del modello, mentre l’analisi dei residui e delle osservazioni influenti ha confermato l’adeguatezza dell’approccio adottato. Le visualizzazioni grafiche hanno inoltre permesso di interpretare in modo chiaro le relazioni tra le variabili, rendendo i risultati facilmente comunicabili anche in contesti clinici.

In conclusione, il modello proposto rappresenta uno strumento predittivo efficace a supporto delle decisioni cliniche, utile per l’individuazione precoce dei neonati a rischio e per la pianificazione ottimale delle risorse ospedaliere.