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.

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

Analisi

Nella prima fase, esploriamo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie. Le variabili “Età della madre”, “numero di gravidanze”, “durata della gravidanza”, “Peso del neonato”, “lunghezza del cranio” e “diametro del cranio”, sono variabili quantitative per cui possono essere analizzate attraverso le statistiche descrittive.

variabili_numeriche <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")
dati_filtered <- dati[, variabili_numeriche]

cat("Summary delle variabili numeriche:\n")

## Summary delle variabili numeriche:

kable(summary(dati_filtered))

Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Min. : 0.00	Min. : 0.0000	Min. :25.00	Min. : 830	Min. :310.0	Min. :235
1st Qu.:25.00	1st Qu.: 0.0000	1st Qu.:38.00	1st Qu.:2990	1st Qu.:480.0	1st Qu.:330
Median :28.00	Median : 1.0000	Median :39.00	Median :3300	Median :500.0	Median :340
Mean :28.16	Mean : 0.9812	Mean :38.98	Mean :3284	Mean :494.7	Mean :340
3rd Qu.:32.00	3rd Qu.: 1.0000	3rd Qu.:40.00	3rd Qu.:3620	3rd Qu.:510.0	3rd Qu.:350
Max. :46.00	Max. :12.0000	Max. :43.00	Max. :4930	Max. :565.0	Max. :390

CV <- function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100

stats_df <- data.frame(
  Variabile = variabili_numeriche,
  Varianza = round(sapply(dati_filtered, var, na.rm = TRUE), 2),
  Deviazione_Standard = round(sapply(dati_filtered, sd, na.rm = TRUE), 2),
  Coefficiente_Variabilita = round(sapply(dati_filtered, CV), 2),
  Asimmetria = round(sapply(dati_filtered, skewness, na.rm = TRUE), 2)
)

cat("\nStatistiche descrittive avanzate:\n")

## 
## Statistiche descrittive avanzate:

kable(stats_df)

	Variabile	Varianza	Deviazione_Standard	Coefficiente_Variabilita	Asimmetria
Anni.madre	Anni.madre	27.81	5.27	18.72	0.04
N.gravidanze	N.gravidanze	1.64	1.28	130.51	2.51
Gestazione	Gestazione	3.49	1.87	4.79	-2.07
Peso	Peso	275665.68	525.04	15.99	-0.65
Lunghezza	Lunghezza	692.67	26.32	5.32	-1.51
Cranio	Cranio	269.79	16.43	4.83	-0.79

L’Età media della madre è di circa 28 anni, con una mediana molto simile, e varia da 0 (che potrebbe essere un errore nei dati) a 46 anni. La distribuzione sembra abbastanza simmetrica (asimmetria vicina a 0) In media le madri hanno avuto circa 1 gravidanza (0.98), il range va da 0 a 12 gravidanze. Ha un’alta asimmetria positiva (2.51) e un elevato coefficiente di variabilità (130.51%), indicando una distribuzione molto sbilanciata verso i valori bassi La durata media di gestazione è di circa 39 settimane, con una variabilità da 25 a 43 settimane e un’asimmetria negativa (-2.07), indicando una coda verso i valori più bassi. La media del peso del neonato è di 3284g, con una mediana simile (3300g). Varia da 830g a 4930g. Presenta una leggera asimmetria negativa (-0.65) La lunghezza e il diametro del cranio entrambe le misure mostrano un’asimmetria negativa, con coefficienti di variabilità relativamente bassi (5.32% e 4.83%) Per meglio rappresentare questi andamenti si riportano anche i boxplot, che permettono tra l’altro di individuare le variabili che hanno dei probabili outlier.

Dai boxplot si può notare la presenza di diversi outlier, particolarmente evidenti per il numero di gravidanze e il peso del neonato, la distribuzione asimmetrica del numero di gravidanze e la relativa simmetria dell’età della madre, che confermano quanto valutato attraverso le statistiche descrittive.

Distribuzioni di Frequenza Variabili Qualitative

variabili_cat <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
for (var in variabili_cat) {
  freq_ass <- table(dati[[var]])
  freq_rel <- prop.table(freq_ass)
  distr_freq <- cbind(Frequenza_Assoluta = freq_ass, Frequenza_Relativa = round(freq_rel, 2))
  print(kable(distr_freq, caption = paste("Distribuzione di Frequenza per", var)))
}

## 
## 
## Table: Distribuzione di Frequenza per Fumatrici
## 
## |   | Frequenza_Assoluta| Frequenza_Relativa|
## |:--|------------------:|------------------:|
## |0  |               2396|               0.96|
## |1  |                104|               0.04|
## 
## 
## Table: Distribuzione di Frequenza per Tipo.parto
## 
## |    | Frequenza_Assoluta| Frequenza_Relativa|
## |:---|------------------:|------------------:|
## |Ces |                728|               0.29|
## |Nat |               1772|               0.71|
## 
## 
## Table: Distribuzione di Frequenza per Ospedale
## 
## |     | Frequenza_Assoluta| Frequenza_Relativa|
## |:----|------------------:|------------------:|
## |osp1 |                816|               0.33|
## |osp2 |                849|               0.34|
## |osp3 |                835|               0.33|
## 
## 
## Table: Distribuzione di Frequenza per Sesso
## 
## |   | Frequenza_Assoluta| Frequenza_Relativa|
## |:--|------------------:|------------------:|
## |F  |               1256|                0.5|
## |M  |               1244|                0.5|

Da queste analisi emerge che solo il 4% delle madri sono fumatrici (104 su 2500), per quanto riguarda il tipo di parto, Il 71% dei parti sono naturali e Il 29% sono cesarei. La distribuzione è quasi uniforme tra i tre ospedali (33-34% ciascuno) I grafici a barre confermano visivamente queste distribuzioni di frequenza. La variabile dipendente è il peso, in prima analisi facciamo un test di Shapiro-Wilk per saggiare l’ipotesi di normalità. Il risultato è di 0.97, con un p-value prossimo a 0, quindi si rifiuta l’ipotesi nulla di normalità della variabile Peso.

Per avere un riscontro, visualizziamo la distribuzione di densità della variabile peso.

Analisi variabile Peso

shapiro.test(Peso)

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

plot(density(Peso), main="Distribuzione del Peso")

La distribuzione del peso mostrata nel grafico appare infatti approssimativamente normale (a campana), anche se si può notare una leggera asimmetria negativa (coda più lunga a sinistra), che è coerente con il coefficiente di asimmetria -0.65 riportato nelle statistiche descrittive. Inoltre c’è da considerare che la dimensione del campione è grande (N = 2500). Con campioni così grandi, il test di Shapiro-Wilk diventa estremamente sensibile anche a piccole deviazioni dalla normalità. Questo può portare al rifiuto dell’ipotesi nulla anche quando la distribuzione è “abbastanza” normale per scopi pratici.

Matrice di Correlazione

dati <- dati[, c("Peso", setdiff(names(dati), "Peso"))]
panel.cor <- function(x, y) {
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- cor(x, y)
  text(0.5, 0.5, format(r, digits=2), cex = 1.5)
}

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

Dall’analisi della matrice si mostra come il peso abbia una correlazione lineare positiva con le variabili Gestazioni, quindi il numero di settimane della gestazione, con un coefficiente di 0.59, con la lunghezza del Cranio (R2 = 0.8) e con il diametro del cranio con un R2 di 0.7. Una leggera correlazione positiva con un R2 di 0.24 è valutata anche con il sesso. Per le restanti variabili i coefficienti di correlazioni sono molto prossimi a zero. Per le variabili di tipo qualitativo per avere un informazione più di dettaglio costruiamo dei boxplot.

variabili_categoriche <- c("Sesso", "Fumatrici", "Tipo.parto", "Ospedale")
par(mfrow = c(2, 2))

for (var in variabili_categoriche) {
  formula <- as.formula(paste("Peso ~", var))
  
  boxplot(formula, data = dati, 
          main = paste("Distribuzione Peso per", var),
          col = rainbow(length(unique(dati[[var]]))),
          xlab = var, ylab = "Peso")
  
  cat("\n", rep("=", 50), "\n")
  cat("ANALISI PER:", var, "\n")
  cat(rep("=", 50), "\n")
  
  if (var == "Ospedale") {
    anova_result <- aov(formula, data = dati)
    cat("ANOVA - Confronto tra ospedali:\n")
    print(summary(anova_result))
    
    if (summary(anova_result)[[1]][["Pr(>F)"]][1] < 0.05) {
      cat("\nTest post-hoc (Tukey HSD):\n")
      print(TukeyHSD(anova_result))
    }
  } else {
    test_result <- t.test(formula, data = dati)
    cat("T-test - Confronto tra gruppi:\n")
    print(test_result)
  }
}

## 
##  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## ANALISI PER: Sesso 
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## T-test - Confronto tra gruppi:
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.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.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

## 
##  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## ANALISI PER: Fumatrici 
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## T-test - Confronto tra gruppi:
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

## 
##  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## ANALISI PER: Tipo.parto 
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## T-test - Confronto tra gruppi:
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -46.27992  40.54037
## sample estimates:
## mean in group Ces mean in group Nat 
##          3282.047          3284.916

## 
##  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## ANALISI PER: Ospedale 
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
## ANOVA - Confronto tra ospedali:
##               Df    Sum Sq Mean Sq F value Pr(>F)
## Ospedale       2    936237  468118   1.699  0.183
## Residuals   2497 687952305  275512

par(mfrow = c(1, 1))

Il peso risulta essere maggiore per i bambini di sesso maschile rispetto al peso delle femmine. Valutando un t.test delle medie che permette di saggiare l’ipotesi di uguaglianza delle medie tra gruppi indipendenti e in questo caso con un p-value prossimo a zero si rifiuta l’ipotesi nulla, quindi le due medie sono significativamente diverse, per questo motivo mi potrò aspettare un beta di regressione per questa variabile signficativo, per cui può aver senso inserire questa variabile come variabile di controllo all’interno del modello. Procedendo allo stesso modo con la variabile fumatrici. Dal grafico si può vedere come i due boxplot non sono molto diversi. Analizzando anche le medie si vede come il peso dei bambini nati da donne fumatrici (3286.153 g) e quello dei bambini nati da donne non fumatrici (3236.346 g) sia molto simile. Infatti anche il t-test con un p-value di 0.303 non si rifiuta l’ipotesi nulla per cui le due medie non sono significativamente diverse. Quindi per questa variabile possiamo aspettarci una scarsa influenza nel nostro modello di regressione. Allo stesso modo anche per il tipo di parto, come si vede dal grafico i due boxplot sono identici e ciò è confermato anche con un p-value del test t-student di 0.90. Per quanto riguarda l’ospedale anche i box plot di seguito riportato mostrano una distribuzione molto simile del peso per i tre ospedali. In questo caso non si può utilizzare il test t di student a causa di più di due livelli della variabile ospedale, per questo motivo si è proceduto con l’uso dell’ANOVA. Dai risultati dell’ANOVA, il valore di F (1.699) e il p-value (0.183) suggeriscono che non ci sono differenze significative nel peso dei neonati tra i diversi ospedali, a un livello di significatività del 5%. In altre parole, non ci sono prove sufficienti per dire che l’ospedale in cui nasce un neonato abbia un effetto significativo sul peso del neonato.

Modello Lineare Multiplo - Iterazioni

cat("=== MODELLO INIZIALE COMPLETO ===\n")

## === MODELLO INIZIALE COMPLETO ===

mod_full <- lm(Peso ~ ., data = dati)
print(summary(mod_full))

## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1124.40  -181.66   -14.42   160.91  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6738.4762   141.3087 -47.686  < 2e-16 ***
## Anni.madre        0.8921     1.1323   0.788   0.4308    
## N.gravidanze     11.2665     4.6608   2.417   0.0157 *  
## Fumatrici       -30.1631    27.5386  -1.095   0.2735    
## Gestazione       32.5696     3.8187   8.529  < 2e-16 ***
## Lunghezza        10.2945     0.3007  34.236  < 2e-16 ***
## Cranio           10.4707     0.4260  24.578  < 2e-16 ***
## Tipo.partoNat    29.5254    12.0844   2.443   0.0146 *  
## Ospedaleosp2    -11.2095    13.4379  -0.834   0.4043    
## Ospedaleosp3     28.0958    13.4957   2.082   0.0375 *  
## SessoM           77.5409    11.1776   6.937 5.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 669.2 on 10 and 2489 DF,  p-value: < 2.2e-16

cat("\n=== SELEZIONE AUTOMATICA DEL MODELLO ===\n")

## 
## === SELEZIONE AUTOMATICA DEL MODELLO ===

mod_step <- step(mod_full, direction = "both", trace = FALSE)
cat("Variabili selezionate automaticamente:", paste(names(coef(mod_step))[-1], collapse = ", "), "\n")

## Variabili selezionate automaticamente: N.gravidanze, Gestazione, Lunghezza, Cranio, Tipo.partoNat, Ospedaleosp2, Ospedaleosp3, SessoM

cat("\n=== MODELLO MANUALE OTTIMIZZATO ===\n")

## 
## === MODELLO MANUALE OTTIMIZZATO ===

modelli <- list(
  mod_full = mod_full,
  mod1 = update(mod_full, ~.-Fumatrici),
  mod2 = update(mod_full, ~.-Fumatrici-Anni.madre),
  mod3 = update(mod_full, ~.-Fumatrici-Anni.madre-Ospedale),
  mod4 = update(mod_full, ~.-Fumatrici-Anni.madre-Ospedale-Tipo.parto),
  mod5 = update(mod_full, ~.-Fumatrici-Anni.madre-Ospedale-Tipo.parto+Lunghezza*Cranio),
  mod6 = update(mod_full, ~.-Fumatrici-Anni.madre-Ospedale-Tipo.parto+Lunghezza*Cranio-Lunghezza),
  mod7 = update(mod_full, ~.-Fumatrici-Anni.madre-Ospedale-Tipo.parto+Lunghezza*Cranio-Lunghezza-Cranio)
)

cat("\n=== CONFRONTO MODELLI ===\n")

## 
## === CONFRONTO MODELLI ===

model_comparison <- data.frame(
  Modello = names(modelli),
  R2 = sapply(modelli, function(x) round(summary(x)$r.squared, 4)),
  R2_adj = sapply(modelli, function(x) round(summary(x)$adj.r.squared, 4)),
  AIC = sapply(modelli, AIC),
  BIC = sapply(modelli, BIC),
  Variabili = sapply(modelli, function(x) length(coef(x)) - 1)
)
print(kable(model_comparison))

## 
## 
## |         |Modello  |     R2| R2_adj|      AIC|      BIC| Variabili|
## |:--------|:--------|------:|------:|--------:|--------:|---------:|
## |mod_full |mod_full | 0.7289| 0.7278| 35171.95| 35241.84|        10|
## |mod1     |mod1     | 0.7288| 0.7278| 35171.15| 35235.22|         9|
## |mod2     |mod2     | 0.7287| 0.7278| 35169.79| 35228.03|         8|
## |mod3     |mod3     | 0.7277| 0.7270| 35175.16| 35221.75|         6|
## |mod4     |mod4     | 0.7270| 0.7265| 35179.33| 35220.10|         5|
## |mod5     |mod5     | 0.7295| 0.7289| 35157.97| 35204.57|         6|
## |mod6     |mod6     | 0.7295| 0.7290| 35155.99| 35196.76|         5|
## |mod7     |mod7     | 0.7257| 0.7252| 35189.37| 35224.31|         4|

mod_finale <- modelli[[which.max(model_comparison$R2_adj)]]
cat("\n=== MODELLO FINALE OTTIMALE ===\n")

## 
## === MODELLO FINALE OTTIMALE ===

print(summary(mod_finale))

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Sesso + 
##     Lunghezza:Cranio, data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1150.9  -180.8   -13.8   165.8  2860.8 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.940e+03  1.828e+02 -10.611  < 2e-16 ***
## N.gravidanze      1.294e+01  4.320e+00   2.995  0.00277 ** 
## Gestazione        3.790e+01  3.700e+00  10.244  < 2e-16 ***
## Cranio           -4.329e+00  7.261e-01  -5.962 2.85e-09 ***
## SessoM            7.338e+01  1.116e+01   6.573 5.97e-11 ***
## Cranio:Lunghezza  3.068e-02  8.871e-04  34.589  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.3 on 2494 degrees of freedom
## Multiple R-squared:  0.7295, Adjusted R-squared:  0.729 
## F-statistic:  1345 on 5 and 2494 DF,  p-value: < 2.2e-16

Analizzando la tabella si vede che il valore del R2 aggiustato sia di 0.73 un valore già abbastanza alto ma che potrebbe sicuramente migliorare in quanto sono presenti nel modello comunque variabili non significative o comunque che possono avere una bassa significatività sul peso. Tali variabili sono: gli anni della madre, madre fumatrici, tipo di parto e l’ospedale. Per quanto riguarda le variabili quantitative, Gestazione, Lunghezza, Cranio e N. Gravidanze mostrano un coefficiente positivo e un p-value fortemente significativo. Analizziamo i risultati del modello lineare multiplo per comprendere l’effetto di ciascuna variabile sul peso del neonato: 1. Variabili altamente significative: Gestazione: Per ogni settimana in più, il peso aumenta in media di 32.57g Lunghezza: Per ogni unità in più, il peso aumenta di 10.29g Cranio: Per ogni unità di diametro in più, il peso aumenta di 10.47g Sesso M: I maschi pesano in media 77.54g più delle femmine 2. Variabili moderatamente significative (p < 0.05, ): N.gravidanze: Ogni gravidanza precedente aggiunge in media 11.27g al peso Tipo.partoNat: I parti naturali sono associati a un aumento di 29.53g rispetto ai cesarei Ospedaleosp3: I neonati nell’ospedale 3 pesano in media 28.10g in più rispetto all’ospedale di riferimento (presumibilmente osp1) 3. Variabili non significative: Età della madre: Effetto non significativo (p = 0.4308) Fumatrici: Effetto negativo (-30.16g) ma non significativo (p = 0.2735) *Ospedaleosp2: Differenza non significativa rispetto all’ospedale di riferimento

Inoltre l’intercetta presenta un valore -6738.48 rappresenta il peso “base” teorico quando tutte le altre variabili sono zero.

Analisi dei Residui

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

plot(density(residuals(mod_finale)), main="Distribuzione dei Residui")
shapiro.test(residuals(mod_finale))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_finale)
## W = 0.96979, p-value < 2.2e-16

bptest(mod_finale)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_finale
## BP = 107.65, df = 5, p-value < 2.2e-16

dwtest(mod_finale)

## 
##  Durbin-Watson test
## 
## data:  mod_finale
## DW = 1.9534, p-value = 0.1219
## alternative hypothesis: true autocorrelation is greater than 0

Il grafico Residui vs Valori Fitted (in alto a sinistra) mostra un pattern a forma di ventaglio che si allarga leggermente all’aumentare dei valori fitted, si evidenziano alcuni outlier, in particolare il punto 1551. La linea rossa (loess) non è perfettamente orizzontale, suggerendo una possibile non-linearità. Il Q-Q Plot dei Residui (in alto a destra), mostra che i punti deviano leggermente dalla linea diagonale ideale, soprattutto alle code suggerendo che i residui non sono perfettamente normali. Ci sono alcuni punti estremi (1551, 1306) che si discostano notevolmente Il terzo grafico “Scale-Location” (in basso a sinistra) mostra che la variabilità dei residui standardizzati non è costante notando un leggero pattern crescente, indicando possibile eteroschedasticità Infine il grafico Residui vs Leverage (in basso a destra), indica che la maggior parte dei punti ha un leverage basso (<0.02). Ci sono alcuni punti ad alto leverage e residui standardizzati elevati e il punto 1551 ha un’influenza particolarmente alta sul modello. Si evince che il test di Shapiro-Wilk per la normalità mostra un p-value < 2.2e-16, è altamente significativo, confermando quanto osservato nel Q-Q plot, i residui non seguono una distribuzione normale. Il Test di Breusch-Pagan è non significativo, conferma l’omoschedasticità che avevamo notato nei grafici. Infine il Test di Durbin-Watson per l’autocorrelazione, è l’unico test non significativo che suggerisce che non c’è una forte autocorrelazione nei residui, il valore di 1.9535 è vicino a 2, che è l’ideale per questo test. Riportiamo la distribuzione dei residui, che mostra un andamento normale.

Leverage, Outlier, Distanza di Cook

lev <- hatvalues(mod_finale)
plot(lev, main="Leverage")
abline(h=2*mean(lev), col="red")

plot(rstudent(mod_finale), main="Outliers")
abline(h=c(-2,2), col="blue")

cook <- cooks.distance(mod_finale)
plot(cook, main="Distanza di Cook")

Si procede ora con la rimozione dei leverage e degli outliers, rimuovendo prima il solo punto 1551 e poi successivamente i punti 155, 155 e 1306 e confrontato i modelli che si ottengono con il modello 4. Otteniamo un miglioramento del R2 aggiustato da 0.73 a 0.74, un leggero miglioramento che però può essere significativo. Confrontiamo anche il BIC, che si riduce a 35055 con il modello senza outlier. Analizziamo nuovamente i residui.

Rimozione Outlier e Confronto Modelli

model_data <- model.frame(mod_finale)
model_data_no_outliers <- model_data[-c(1551, 155, 1306),]
mod_finale_no_outliers <- update(mod_finale, data = model_data_no_outliers)
summary(mod_finale_no_outliers)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Sesso + 
##     Lunghezza:Cranio, data = model_data_no_outliers)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1173.17  -177.57   -12.21   166.89  1162.83 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.547e+03  1.796e+02  -8.615  < 2e-16 ***
## N.gravidanze      1.459e+01  4.177e+00   3.493 0.000486 ***
## Gestazione        3.464e+01  3.586e+00   9.658  < 2e-16 ***
## Cranio           -6.439e+00  7.223e-01  -8.914  < 2e-16 ***
## SessoM            7.220e+01  1.080e+01   6.687 2.81e-11 ***
## Cranio:Lunghezza  3.334e-02  8.841e-04  37.714  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 264.2 on 2491 degrees of freedom
## Multiple R-squared:  0.7462, Adjusted R-squared:  0.7457 
## F-statistic:  1465 on 5 and 2491 DF,  p-value: < 2.2e-16

BIC(mod_finale, mod_finale_no_outliers)

##                        df      BIC
## mod_finale              7 35196.76
## mod_finale_no_outliers  7 34984.35

Previsione Peso Neonato

new_born <- data.frame(N.gravidanze = 3, Gestazione = 39, Sesso = "F")
mod_to_predict <- lm(Peso ~ N.gravidanze + Gestazione + Sesso, data = dati)
predict(mod_to_predict, newdata = new_born)

##        1 
## 3252.311

Visualizzazioni Finali

dati$Fumatrici <- factor(dati$Fumatrici, levels = c(0, 1), labels = c("Non Fumatrice", "Fumatrice"))

ggplot(dati) +
  geom_point(aes(x = Lunghezza, y = Peso, col = Sesso), position = "jitter") +
  geom_smooth(aes(x = Lunghezza, y = Peso, col = Sesso), se = FALSE, method = "lm") +
  labs(title = "Peso vs. Lunghezza per Sesso") +
  theme_minimal()

ggplot(dati) +
  geom_point(aes(x = N.gravidanze, y = Peso, col = Fumatrici), position = "jitter") +
  geom_smooth(aes(x = N.gravidanze, y = Peso, col = Fumatrici), se = FALSE, method = "lm") +
  labs(title = "Peso vs. N. Gravidanze per Fumatrici") +
  theme_minimal()

Grafico - Peso vs. Lunghezza per Sesso Lo scatter plot mostra una forte correlazione positiva tra lunghezza e peso del neonato, con alcune osservazioni interessanti: Relazione generale: Entrambi i sessi mostrano una correlazione lineare molto chiara: all’aumentare della lunghezza aumenta proporzionalmente il peso Le linee di regressione sono quasi parallele, indicando che la relazione lunghezza-peso è simile tra maschi e femmine

Differenze tra sessi: I maschi (blu) tendono ad avere un peso leggermente superiore alle femmine a parità di lunghezza La differenza è consistente lungo tutto il range di lunghezze, confermando il coefficiente positivo del sesso maschile nel modello (+77g circa) *La distribuzione è abbastanza omogenea senza particolari outlier evidenti

Grafico 2 - Peso vs. N. Gravidanze per Fumatrici Questo grafico presenta un pattern molto diverso e più complesso.La maggior parte dei dati si concentra nelle prime gravidanze (0-3), con una netta diminuzione per gravidanze successive Ci sono pochissimi casi oltre la 4ª gravidanza, rendendo meno affidabili le stime per valori alti.

L’effetto del fumo: Le non fumatrici (rosso) dominano numericamente il dataset Le fumatrici (blu) sono molto meno rappresentate *Le linee di regressione mostrano trend opposti: leggermente crescente per non fumatrici, decrescente per fumatrici

Modello Statistico per la Previsione del Peso Neonatale

Domenico Toscano

2025-06-12

Raccolta dei Dati