Introduzione e Contesto

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.

Per attuare ciò si utilizza un dataset raccolto con dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

Età della madre: Misura dell’età in anni.
Numero di gravidanze: Quante gravidanze ha avuto la madre.
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).
Durata della gravidanza: Numero di settimane di gestazione.
Peso del neonato: Peso alla nascita in grammi.
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
Tipo di parto: Naturale o cesareo.
Ospedale di nascita: Ospedale 1, 2 o 3.
Sesso del neonato: Maschio (M) o femmina (F).

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

Caricamento del Dataset

# Importazione del dataset e trasformazione delle variabili categoriche in fattori

dati <- read.csv("neonati.csv", 
                 sep=",", 
                 stringsAsFactors = TRUE, 
                 fileEncoding = "Latin1")

# Collegamento al dataset
attach(dati)
n <- nrow(dati)

# Statistiche descrittive
summary(dati)

##    Anni.madre     N.gravidanze       Fumatrici        Gestazione   
##  Min.   : 0.00   Min.   : 0.0000   Min.   :0.0000   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000   Median :0.0000   Median :39.00  
##  Mean   :28.16   Mean   : 0.9812   Mean   :0.0416   Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.0000   3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000   Max.   :1.0000   Max.   :43.00  
##       Peso        Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   : 830   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :3300   Median :500.0   Median :340              osp3:835           
##  Mean   :3284   Mean   :494.7   Mean   :340                                 
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :4930   Max.   :565.0   Max.   :390

Analisi della Variabile Risposta

Per costruire un modello di previsione del peso neonatale, consideriamo il peso alla nascita (Peso) come variabile dipendente (Y), mentre tutte le altre variabili saranno considerate variabili esplicative (X).

L’obiettivo è individuare quali fattori influenzano significativamente il peso del neonato da un punto di vista statistico.

La prima analisi si concentra sulla distribuzione della variabile Peso, per verificare se segue una distribuzione normale, in modo da facilitare le analisi di regressione.

Utilizziamo la libreria moments per calcolare la asimmetria (skewness) e la curtosi (kurtosis) della distribuzione, oltre ad eseguire il test di Shapiro-Wilk per verificare la normalità.

# Calcolo dell'asimmetria (skewness)
library(moments)
moments::skewness(Peso)

## [1] -0.6470308

# Calcolo della curtosi (kurtosis) centrata sulla distribuzione normale
moments::kurtosis(Peso) - 3

## [1] 2.031532

# Test di Shapiro-Wilk per la normalità
test_shapiro <- shapiro.test(Peso)
test_shapiro

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

# Visualizzazione della densità della distribuzione del peso
plot(density(Peso), 
     main = "Distribuzione della Variabile Peso", 
     xlab = "Peso (g)", 
     col = "blue", 
     lwd = 2)

L’analisi degli indici di forma evidenzia che la variabile Peso presenta un’asimmetria leggermente negativa (-0.647), indicando una lieve tendenza verso valori inferiori rispetto alla distribuzione normale. La curtosi di circa 2.03 suggerisce una distribuzione leptocurtica, ovvero caratterizzata da code più pronunciate rispetto a una normale.

Il grafico della densità conferma quanto emerso dalle misure di skewness e kurtosis, mostrando una distribuzione non perfettamente normale. Inoltre, il test di Shapiro-Wilk fornisce un p-value molto basso, portandoci a rifiutare l’ipotesi nulla di normalità.

Nonostante queste evidenze, la distribuzione osservata non si discosta in modo eccessivo da una normale, risultando solo leggermente leptocurtica. Pertanto, procediamo con l’analisi senza applicare trasformazioni ai dati.

Analisi delle Correlazioni

Dopo aver studiato la distribuzione della variabile risposta, esaminiamo le relazioni tra le variabili presenti nel dataset.

Utilizziamo la funzione pairs() per rappresentare graficamente le correlazioni tra le variabili numeriche, fornendo sia una rappresentazione quantitativa che visiva.

# Funzione personalizzata per visualizzare le correlazioni nei grafici a coppie
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)
}

# Grafico delle correlazioni
pairs(dati, upper.panel=panel.smooth, lower.panel = panel.cor)

Poiché la funzione pairs() non è particolarmente utile per l’interpretazione delle variabili categoriche, utilizzeremo successivamente dei boxplot per analizzare la relazione tra queste variabili e il peso neonatale.

Dall’output del pairs(), osserviamo che le correlazioni più elevate si riscontrano tra:

Peso e Lunghezza (0.80)
Peso e Cranio (0.70)
Peso e Gestazione (0.59)
Gestazione e Lunghezza (0.62)
Lunghezza e Cranio (0.60)
Gestazione e Cranio (0.46)

Inoltre, la parte superiore del grafico evidenzia la presenza di relazioni non strettamente lineari tra alcune variabili. In particolare, si nota una curvatura nelle relazioni tra: - Peso e Gestazione - Gestazione e Lunghezza - Gestazione e Cranio - Peso e Cranio - Peso e Lunghezza

Queste evidenze suggeriscono la possibilità di introdurre termini quadratici nelle successive analisi di regressione per catturare meglio la relazione tra queste variabili.

A parte ciò, l’analisi delle correlazioni per quanto alcune possano essere maggiormente evidenti non sembrano esserci particolari problemi di multicollinearità che comunque verranno approfonditi in seguito.

Analisi delle Variabili Categoriali

Per le variabili categoriche utilizziamo i boxplot per visualizzare la distribuzione del peso in funzione di: - Sesso - Fumo materno - Tipo di parto

par(mfrow = c(2, 2))
boxplot(Peso, main = "Distribuzione generale del Peso")
boxplot(Peso ~ Sesso, main = "Peso per Sesso")
boxplot(Peso ~ Fumatrici, main = "Peso per Fumatrici")
boxplot(Peso ~ Tipo.parto, main = "Peso per Tipo di Parto")

Dai boxplot osserviamo che: - Il peso dei neonati maschi è mediamente più alto rispetto alle femmine. - Il fumo materno non sembra avere un impatto significativo sulla distribuzione del peso. - Non emergono differenze rilevanti tra i neonati nati con parto cesareo e quelli nati con parto naturale.

Per confermare queste osservazioni, applichiamo un t-test per confrontare le medie.

# Test t per confronto tra gruppi
t.test(Peso ~ Sesso)

## 
##  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

t.test(Peso ~ Fumatrici)

## 
##  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

t.test(Peso ~ Tipo.parto)

## 
##  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

I risultati confermano: - Una differenza statisticamente significativa tra il peso dei maschi e delle femmine (p-value < 2.2e-16). - Nessuna differenza significativa tra neonati di madri fumatrici e non (p-value = 0.3033). - Nessuna differenza significativa tra parto cesareo e naturale (p-value = 0.8968).

Questo suggerisce che il sesso del neonato potrebbe essere incluso nel modello predittivo, mentre il fumo materno e il tipo di parto potrebbero non essere rilevanti.

Costruzione del Modello di Regressione Multipla

Adotteremo un approccio stepwise backward, partendo da un modello completo con tutte le variabili e rimuovendo progressivamente quelle non significative.

# Modello iniziale con tutte le variabili esplicative
mod1 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, 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

L’output del primo modello mostra che l’R² Adjusted è pari a 0.7278, indicando già una buona capacità predittiva. Tuttavia, alcune variabili come Anni.madre non risultano statisticamente significative.

Procediamo quindi rimuovendo la variabile Anni.madre e rivalutando il modello.

# Rimozione della variabile Anni.madre
mod2 <- update(mod1, ~ . - Anni.madre)
summary(mod2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.93  -180.11   -16.36   160.58  2616.96 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6708.1065   135.9394 -49.346  < 2e-16 ***
## N.gravidanze     12.6085     4.3381   2.906  0.00369 ** 
## Fumatrici       -30.3092    27.5359  -1.101  0.27113    
## Gestazione       32.2501     3.7968   8.494  < 2e-16 ***
## Lunghezza        10.2944     0.3007  34.239  < 2e-16 ***
## Cranio           10.4876     0.4255  24.651  < 2e-16 ***
## Tipo.partoNat    29.5351    12.0834   2.444  0.01458 *  
## Ospedaleosp2    -11.0816    13.4359  -0.825  0.40957    
## Ospedaleosp3     28.3660    13.4903   2.103  0.03559 *  
## SessoM           77.6205    11.1763   6.945 4.81e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2490 degrees of freedom
## Multiple R-squared:  0.7288, Adjusted R-squared:  0.7278 
## F-statistic: 743.6 on 9 and 2490 DF,  p-value: < 2.2e-16

L’output del modello 2 mostra che l’R² Adjusted rimane sostanzialmente invariato (0.7278), confermando che l’esclusione della variabile Anni.madre non ha influenzato negativamente la capacità predittiva del modello. Questo supporta la decisione di eliminare la variabile per migliorare la parsimonia del modello.

Dopo aver rimosso Anni.madre, passiamo alla valutazione della variabile Fumatrici, che non è risultata significativa nel modello.

# Rimozione della variabile Fumatrici
mod3 <- update(mod2, ~ . - Fumatrici)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## N.gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16

Dall’output osserviamo che: - L’R² Adjusted rimane invariato (0.7278), indicando che la rimozione della variabile non ha avuto impatti negativi. - I coefficienti delle altre variabili restano stabili, confermando che Fumatrici non aveva un’influenza significativa sul peso neonatale.

A questo punto, ci concentriamo sulla variabile Ospedale. Anche se la differenza tra Ospedale 1 e Ospedale 3 è statisticamente significativa (p = 0.037), il passaggio da Ospedale 1 a Ospedale 2 non lo è (p = 0.412). Tuttavia, il contesto suggerisce che il peso del neonato non dovrebbe dipendere significativamente dalla struttura ospedaliera.

Decidiamo quindi di rimuovere completamente la variabile Ospedale e verificare l’impatto sul modello.

# Rimozione della variabile Ospedale
mod4 <- update(mod3, ~ . - Ospedale)
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1129.31  -181.70   -16.31   161.07  2638.85 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.2971   135.9911 -49.322  < 2e-16 ***
## N.gravidanze     12.7558     4.3366   2.941   0.0033 ** 
## Gestazione       32.2713     3.7941   8.506  < 2e-16 ***
## Lunghezza        10.2864     0.3007  34.207  < 2e-16 ***
## Cranio           10.5057     0.4260  24.659  < 2e-16 ***
## Tipo.partoNat    30.0342    12.0969   2.483   0.0131 *  
## SessoM           77.9285    11.1905   6.964 4.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared:  0.7277, Adjusted R-squared:  0.727 
## F-statistic:  1110 on 6 and 2493 DF,  p-value: < 2.2e-16

Come è possibile osservare, l’R^2 Adjusted rimane sostanzialmente invariato rispetto ai modelli precedenti, attestandosi a 0.727. Titte le variabili, inoltre, rimaste nel modello risultano significative, indicando che il modello mantiene una buona capacità predittiva. Come detto in precedenza, la decisione di eliminare la variabile Ospedale si basa più su considerazioni teoriche che statistiche, in quanto non ci si aspetta che l’ospedale di nascita possa influenzare il peso del neonato in modo determinante.

Per lo stesso ragionamento fatto con la variabile Ospedale, proviamo a rimuovere Tipo.Parto, in quanto è poco plausibile che il tipo di parto (naturale o cesareo) abbia un impatto significativo sul peso del neonato.

# Rimozione della variabile Tipo.parto
mod5 <- update(mod4,~.-Tipo.parto) 
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Il Multiple R² resta 0.727, mentre l’Adjusted R² è 0.7265, indicando che la capacità predittiva del modello non ha subito perdite significative, mantenendo un ottimo livello di spiegazione della variabilità del peso neonatale.

Infine, per lo stesso ragionamento tenitamo di rimuovere la variabile N.gravidanze anche se significativa in quanto realisticamente potrebbe avere un impatto minimo.

# Rimozione della variabile N.gravidanze
mod6 <- update(mod5,~.-N.gravidanze) 
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16

Rimuovendo la variabile N.gravidanze, il modello finale conserva le variabili Gestazione, Lunghezza, Cranio e Sesso, tutte altamente significative. Il valore di R² Adjusted passa da 0.7265 a 0.7257, una variazione minima che suggerisce come N.gravidanze non aggiungesse un contributo rilevante alla spiegazione del peso neonatale. Anche l’errore standard dei residui rimane stabile, confermando che la qualità della stima non ha subito peggioramenti significativi. Considerando questi elementi, possiamo adottare questo modello finale, che risulta più parsimonioso senza compromettere la capacità predittiva.

Il modello di regressione multipla finale può essere scritto come:

\[ \hat{Peso} = \beta_0 + \beta_1 Gestazione + \beta_2 Lunghezza + \beta_3 Cranio + \beta_4 Sesso + \epsilon \]

Dove: - \(\beta_0\) è l’intercetta, - \(\beta_1, \beta_2, \beta_3, \beta_4\) sono i coefficienti stimati per le variabili esplicative, - \(\epsilon\) rappresenta l’errore residuo.

I coefficienti stimati dal modello finale sono:

\[ \hat{Peso} = -6651.1188 + 31.2737 \cdot Gestazione + 10.2054 \cdot Lunghezza + 10.6704 \cdot Cranio + 79.1049 \cdot SessoM + \epsilon \]

Dove: - Gestazione è il numero di settimane di gravidanza, - Lunghezza è la lunghezza del neonato alla nascita, - Cranio è la circonferenza cranica del neonato, - SessoM assume valore 1 se il neonato è maschio, 0 se femmina.

Tutti i coefficienti sono statisticamente significativi (\(p < 0.001\)), confermando l’importanza di queste variabili nella previsione del peso neonatale.

Effetto quadratico

Tuttavia, all’inizio del progetto avevamo sottolineato come alcune relazioni tra variabili potevano presupporre una relazione quadratica, quindi una convessità o concavità. Per tale motivo, abbiamo proseguito con l’aggiornamento del modello andando a verificare se l’aspetto quadratico della variabile Gestazione, Lunghezza e Cranio può migliorare il modello.

# Aggiunta di Gestazione^2
mod7 <- update(mod6,~.+I(Gestazione^2)) 
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1132.84  -181.45   -15.99   162.62  2649.78 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4643.6094   899.4989  -5.162 2.63e-07 ***
## Gestazione        -80.7957    49.7863  -1.623   0.1047    
## Lunghezza          10.3087     0.3039  33.920  < 2e-16 ***
## Cranio             10.7663     0.4262  25.259  < 2e-16 ***
## SessoM             76.9359    11.2436   6.843 9.75e-12 ***
## I(Gestazione^2)     1.4960     0.6627   2.258   0.0241 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7267, Adjusted R-squared:  0.7261 
## F-statistic:  1326 on 5 and 2494 DF,  p-value: < 2.2e-16

Con l’introduzione del termine quadratico per la Gestazione, notiamo un cambiamento interessante: la variabile Gestazione perde significatività, mentre il suo termine al quadrato diventa significativo. Questo suggerisce che l’effetto della gestazione sul peso non è strettamente lineare, ma segue una curva (come avevamo intravisto nel plot delle correlazioni)

A livello interpretativo, ciò significa che all’aumentare delle settimane di gestazione, il peso cresce, ma il ritmo di crescita potrebbe non essere costante. Inizialmente potrebbe aumentare in modo più marcato, per poi rallentare nelle ultime settimane. Questa dinamica ha senso: nei primi mesi di gravidanza il feto cresce rapidamente, # ma verso la fine della gestazione il tasso di crescita può stabilizzarsi. Il modello sembra quindi catturare una relazione più realistica tra le settimane di gestazione e il peso, nonostante la variabile lineare perde di significatività. Un altro problema, inotlre, è che per ora il modello mostra una relazione inversa rispetto a quella presupposta dalla realtà che abbiamo appena spiegato, in quanto la variabile lineare assume un valor enegativo (-80) e leggermente positiva quella quadratica (+1.496). Tuttavia prima di trarre conclusioni, anche considerando il fatto che l’R² Adjusted rimane pressochè uguale) verifichiamo se anche la variabile lunghezza segue un andamento simile e vediamo come variano le significatività.

# Aggiunta di Lunghezza^2
mod8 <- update(mod7,~.+I(Lunghezza^2)) 
summary(mod8)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2) + I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1177.46  -182.44    -9.66   164.95  1406.21 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.384e+03  9.072e+02  -2.628  0.00865 ** 
## Gestazione       3.310e+02  6.281e+01   5.270 1.48e-07 ***
## Lunghezza       -3.161e+01  4.043e+00  -7.820 7.76e-15 ***
## Cranio           1.060e+01  4.177e-01  25.376  < 2e-16 ***
## SessoM           7.397e+01  1.101e+01   6.716 2.30e-11 ***
## I(Gestazione^2) -3.819e+00  8.261e-01  -4.624 3.97e-06 ***
## I(Lunghezza^2)   4.310e-02  4.145e-03  10.398  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.1 on 2493 degrees of freedom
## Multiple R-squared:  0.738,  Adjusted R-squared:  0.7374 
## F-statistic:  1170 on 6 and 2493 DF,  p-value: < 2.2e-16

L’inserimento della variabile Lunghezza^2 riequilibria da un putno di vista cocnettuale quanto abbiamo detto nella precedente sezione. La variabile Gestazione torna significativa e inverte la relazione con la variabile quadratica (leggermente positivo il coefficiente per la variabile lineare e leggermente negativo per la variabile quadratica). Inoltre, la variabile Lunghezza è significativa sia nella forma lineare che quadratica con un’interpretazione opposta rispetto a Gestazione, in quanto inizialmente all’aumentare dei mm di lunghezza il peso è negativo (-31) mentre risulta positivo nella forma quadratica (presupponendo il fatto che se i mm aumentano tanto questo avrà sicuramente un impatto sul peso). Realisticamente parlando può essere considerato non realistica una interpretazione del genere, in quanto concettualmente all’aumentare della lunghezza del bambino il peso dovrebbe aumentare. Al netto di ciò, potrebbe anche essere che inizialmente magari il bambino ancora non è nemmeno formato e quindi la lunghezza può essere un avariabile che “sballa” leggermente, mentre nel suo aspetto quadratico torna la posiitività. Non mi sento di toglierla in quanto la rimozione di questa variabile crea disagi nelle altre variabili, quindi per una questione di trade-off tra statistica e realtà proseguo lasciandola dentro. Considerando, infine, anche che l’R² Adjusted che migliora superando il 73% rispetto al precedente modello, crediamo sia un buon modello in temirni di trade off tra realtà e statistica, come detto poc’anzi.

Un ultimo tentativo è stato eseguito con l’inserimento anche della quadraticità della variabile Cranio

# Aggiunta di Cranio^2
mod9 <- update(mod8,~.+I(Cranio^2))
summary(mod9)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1173.39  -183.29   -12.16   163.59  1473.80 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.214e+03  1.163e+03  -1.044    0.296    
## Gestazione       3.709e+02  6.751e+01   5.494 4.33e-08 ***
## Lunghezza       -2.865e+01  4.442e+00  -6.448 1.35e-10 ***
## Cranio          -5.181e+00  9.825e+00  -0.527    0.598    
## SessoM           7.366e+01  1.101e+01   6.690 2.75e-11 ***
## I(Gestazione^2) -4.332e+00  8.852e-01  -4.894 1.05e-06 ***
## I(Lunghezza^2)   4.008e-02  4.550e-03   8.810  < 2e-16 ***
## I(Cranio^2)      2.328e-02  1.448e-02   1.608    0.108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269 on 2492 degrees of freedom
## Multiple R-squared:  0.7383, Adjusted R-squared:  0.7376 
## F-statistic:  1004 on 7 and 2492 DF,  p-value: < 2.2e-16

L’aggiunta della variabile Cranio^2 non ha prodotto risultati statistici validi in quanto perdono la significatività entrambe le variabili, nonostante un R² Adjusted che rimane in linea.

Per togliere ogni dubbio sul fatto che la variabile Cranio possa essere non considerata all’interno del modello dopo che comunque abbiamo visto che la correlazione è evidente con la variabile risposta, è stato testato il modello senza quest’ultima.

# Rimozione di Cranio e Cranio^2
mod9 <- update(mod8,~.-Cranio)
summary(mod9)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Sesso + I(Gestazione^2) + 
##     I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1187.44  -191.62   -18.77   188.10  2238.03 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.970e+03  1.017e+03  -2.920  0.00353 ** 
## Gestazione       5.056e+02  7.002e+01   7.221 6.80e-13 ***
## Lunghezza       -3.234e+01  4.534e+00  -7.133 1.29e-12 ***
## SessoM           8.705e+01  1.234e+01   7.055 2.23e-12 ***
## I(Gestazione^2) -5.960e+00  9.216e-01  -6.467 1.20e-10 ***
## I(Lunghezza^2)   4.716e-02  4.645e-03  10.152  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 301.8 on 2494 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.6697 
## F-statistic:  1014 on 5 and 2494 DF,  p-value: < 2.2e-16

Come è possibile osservare dall’output notiamo che nonostante tutte le variabili risultano essere significative l’R² Adjusted scende di molto (passando da circa 73% a circa 67%). Per tale motivo scegliamo come modello il numero 8.

Riepilogando in notazione statistica:

Notazione statistica del modello

\[ Peso_i = \beta_0 + \beta_1 \cdot Gestazione_i + \beta_2 \cdot Lunghezza_i + \beta_3 \cdot Cranio_i + \beta_4 \cdot SessoM_i + \beta_5 \cdot Gestazione_i^2 + \beta_6 \cdot Lunghezza_i^2 + \beta_7 \cdot Cranio_i^2 + \epsilon_i \]

Sostituendo i coefficienti stimati dal modello:

\[ Peso_i = -1214 + (370.9 \cdot Gestazione_i) + (-28.65 \cdot Lunghezza_i) + (-5.18 \cdot Cranio_i) + (73.66 \cdot SessoM_i) + (-4.332 \cdot Gestazione_i^2) + (0.402 \cdot Lunghezza_i^2) + (2.328 \cdot Cranio_i^2) + \epsilon_i \]

A questo punto procediamo come delle metodologie statistiche che ci aiutano a confermare o meno quanto abbiamo deciso sul modello 8 che ragionevolmente è il migliore. L’ANOVA test fornisce un’idea su quale modello riuslta essere statisticamente migliore considerando l’aggiunta di nuove variabili e conducendo il test di Fisher che considera simultaneamente tutte le variabili esplicative.

# ANOVA
anova(mod1,mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9)

## Analysis of Variance Table
## 
## Model 1: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## Model 3: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## Model 4: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## Model 5: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 6: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 7: Peso ~ Gestazione + Lunghezza + Cranio + Sesso + I(Gestazione^2)
## Model 8: Peso ~ Gestazione + Lunghezza + Cranio + Sesso + I(Gestazione^2) + 
##     I(Lunghezza^2)
## Model 9: Peso ~ Gestazione + Lunghezza + Sesso + I(Gestazione^2) + I(Lunghezza^2)
##   Res.Df       RSS Df Sum of Sq        F    Pr(>F)    
## 1   2489 186762521                                    
## 2   2490 186809099 -1    -46578   0.6207  0.430846    
## 3   2491 186899996 -1    -90897   1.2114  0.271162    
## 4   2493 187601677 -2   -701680   4.6757  0.009401 ** 
## 5   2494 188065546 -1   -463870   6.1820  0.012971 *  
## 6   2495 188688687 -1   -623141   8.3047  0.003988 ** 
## 7   2494 188303891  1    384797   5.1282  0.023626 *  
## 8   2493 180477124  1   7826767 104.3080 < 2.2e-16 ***
## 9   2494 227092921 -1 -46615797 621.2527 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dall’output otteiamo conferma sul modello 8, in quanto nonostante il modello 9 ha un RSS leggermente minore non supera il test di Fisher.

Procediamo anche con l’AIC e il BIC test che sono anch’essi criteri di selezione del modello usati per confrontare modelli diversi bilanciando adattamento e complessità. Mentre l’AIC criteria va a prefeerire modelli più complessi (con più regressori), il BIC tende a preferire modelli più semplici.

# AIC e BIC test
AIC(mod1, mod2, mod3,mod4, mod5, mod6, mod7, mod8, mod9)

##      df      AIC
## mod1 12 35171.95
## mod2 11 35170.57
## mod3 10 35169.79
## mod4  8 35175.16
## mod5  7 35179.33
## mod6  6 35185.60
## mod7  7 35182.50
## mod8  8 35078.36
## mod9  7 35650.75

BIC(mod1, mod2, mod3,mod4, mod5,mod6, mod7, mod8, mod9)

##      df      BIC
## mod1 12 35241.84
## mod2 11 35234.64
## mod3 10 35228.03
## mod4  8 35221.75
## mod5  7 35220.10
## mod6  6 35220.54
## mod7  7 35223.27
## mod8  8 35124.96
## mod9  7 35691.52

Considerando che l’AIC più o meno fornisce un valore uguale trA il modello 8 e il 9 che sono quelli ritenuti migliori, il BIC mostra una differenza più sostanziale tra i due. I due criteri confermano quanto detto in precedenza sul modello 8.

MULTICOLLINEARITA’

Per verificare se ci sono variabili affetti da multicollinearità ci rifacciamo alla misura del VIF (Variance Inflation Factor). Per non avere multicollinearità si preferisce che il valore sia inferiore a 5.

# VIF (Variance Inflation Factor)
library(car)
vif(mod8)

##      Gestazione       Lunghezza          Cranio           Sesso I(Gestazione^2) 
##      475.567702      390.760206        1.624664        1.047157      453.369759 
##  I(Lunghezza^2) 
##      372.185476

L’output mostra valori alti (ordine di 300 e 400) per le variabili Gestazione, Lunghezza, Gestazione^2 e Lunghezza^2, mentre valori molto bassi e accettabili per Cranio e Lunghezza. Tuttavia, il fatto che ci siano valori alti per le variabili sopra citate risulta esser enormale in quanto sono presenti le variabili originali e le stesse al quadrato. Per tale motivo, possiamo procedere.

Per levare ogni dubbio è stato applicato il VIF al modello 6 che è lo stesso del modello 8 senza le variabili quadratiche.

# VIF (Variance Inflation Factor)
library(car)
vif(mod6)

## Gestazione  Lunghezza     Cranio      Sesso 
##   1.653502   2.069517   1.606131   1.038813

L’Output conferma quanto descitto di sopra, in quanto tutte le variabili mostrano valori al di sotto di 5.

Diagnosi dei residui - Leverage e Outliers e Test sui residui (Omoschedasticità, Autocorrelazione, Normalità)

#Plot dei residui
plot(mod8)

L’output mostra i grafici diagnostici dei residui del modello selezionato (modello 8) e ogni grafico mostrato serve per osservare graficamente le ipotesi fondamentali della regressione:

Media dei residui uguale a 0
Omoschedaticità (varianza dei residui costante)
No autocorrelazione
Distribuzione Normale dei residui

Il primo grafico (Residuals vs Fitted) aiuta a verificare l’ipotesi di linearità e l’omoscedasticità. Quello che osserviamo dal grafico è che i residui sembrano distribuiti in modo relativamente casuale anche se ci sono alcuni outlier identificati (come 1551, 1306 e 1399). Il fatto, inoltre, che alcuni punti si discostano molto dalla nuvola centrale suggerisce possibili violazioni dell’omoschedasticità.

Il secondo grafico (Q-Q Residuals) controlla se i residui seguono una distribuzione normale, dove i vari punti dovrebbero seguire la bisettrice diagonale.I residui si allineano abbastanza bene alla diagonale, ma ci sono deviazioni nelle code (es. osservazioni come 1551 e 1306). Questo indica che ci sono outlier o che la distribuzione dei residui potrebbe avere code più pesanti della normale.

Il terzo grafico (Scale-Location) serve per verificare l’omoschedasticità dove per rispettare l’ipotesi i punti dovrebbero essere distribuiti in modo omogeneo senza un trend evidente. Tuttavia, il grafico mostra una leggera tendenza crescente nella dispersione (cioè il valore dei residui tende ad aumentare per valori predetti più alti).

Il quarto grafico (Residuals vs Leverage) aiuta ad identificare punti influenti nel modello. Per non avere punti influenti, tutti i valori dovrebbero più o meno giaccere sulla linea orizzontale ed esser elontani dalle soglie identificate dalla distanza di Cook (0.5 e 1). Come si può osservare dal grafico, l’osservazione 1551 sembrerebbe essere un punto influete (alto leverage e residuo standardizzato elevato).

ANALISI DEI VALORI DI LEVA E OUTLIERS

Quando si analizzano i dati in un modello di regressione, è importante distinguere tra valori leverage e outliers, poichè indicano problemi diversi nel modello.

Il leverage misura quanto un’osservazione è distante dal centro dei dati (media) rispetto alle variabili esplicative (X), dove un’osservazione con alto leverage ha un valore di X molto distante dalla media, il che significa che potrebbe avere un’influenza sproporzionata sulla regressione. Tramite la funzione hatvalues riesco a calcolare il leverage per ogni osservazione. Si ricorda, tuttavia, che un leverage alto, non significa necessariamente che il punto è problematico, ma solo che si trova in una regione con pochi altri punti.

L’Outlier è un’osservazione, invece, che ha un valore della variabile dipendente (Y) molto distante rispetto a quanto previsto dal modello e sono identificabili con i residui standardizzati o studentizzati (rstudent). Quindi con la funzione rstudent si calcolano i residui standardizzati del modello di ogni osservazione e con la funzione OutlierTest si confronta il valore del residuo con un nuovo valore soglia p che tiene conto del numero di test effettuati (chiamasi anche Aggiustamento di Bonferroni). Fondamentalmente, avendo nel dataset 2500 osservazioni viene effettuato un t-test sulla significatività del reiduo standardizzato ma il valore soglia al 5% (alpha) viene diviso per 2500 (in questo caso) così che risultano veramente outliers quelli che hanno un p-value al di sotto di 0.05/2500.

La Cook’s Distance misura quanto il modello cambierebbe se rimuovessimo un’osservazione. E’ una metrica che misura quanto la rimozione di una singola osservazione cambierebbe i coefficienti della regressione. Un’ooservazione può essere un outlier ma ininfluente sui coefficienti di regressione e quindi potrebbe essere tranquillamente tenuta. Le soglie tipiche della Cook’s Distance sono 0.5 e 1.0.

par(mfrow=c(2,2))
#valori di leva e outlier
lev <- hatvalues(mod8)
plot(lev)
p = sum(lev)
soglia = 2*p/n
abline(h=soglia, col=2)

lev[lev>soglia]

##          15          34          36          61          67         101 
## 0.007884257 0.006286413 0.007367715 0.006973015 0.007988541 0.017041152 
##         106         119         131         151         155         165 
## 0.022328004 0.005982011 0.007581362 0.010927167 0.012328598 0.005633158 
##         206         220         304         305         310         312 
## 0.013573772 0.007681754 0.008515946 0.008062414 0.092606575 0.038557719 
##         317         322         328         378         383         428 
## 0.006569069 0.005829405 0.006988231 0.057856611 0.007630319 0.005633158 
##         445         471         486         492         587         592 
## 0.008505526 0.010015281 0.006778847 0.020294394 0.020571069 0.011798253 
##         629         638         666         684         697         702 
## 0.006160189 0.006402403 0.008760601 0.008786380 0.005954875 0.007428351 
##         729         748         750         765         805         895 
## 0.008376656 0.014633423 0.012722154 0.007920693 0.034721336 0.007290926 
##         928         947         956         958         964         988 
## 0.115239342 0.012042213 0.009005388 0.006513797 0.007817221 0.006554396 
##         991        1014        1067        1091        1130        1181 
## 0.007143099 0.027842764 0.015040313 0.009637166 0.011482376 0.005965660 
##        1188        1200        1215        1222        1238        1248 
## 0.009482560 0.005761662 0.005692630 0.005982011 0.006872665 0.025840065 
##        1267        1273        1275        1283        1293        1323 
## 0.006194008 0.012029692 0.006807527 0.007178150 0.005746074 0.008376656 
##        1356        1357        1385        1400        1428        1429 
## 0.005703204 0.012504263 0.039739259 0.007791495 0.029494489 0.043907337 
##        1513        1551        1553        1556        1593        1610 
## 0.009303269 0.247420148 0.007120529 0.006084760 0.006650003 0.011972178 
##        1619        1628        1634        1639        1686        1701 
## 0.063579635 0.008560232 0.007140466 0.006554396 0.018192107 0.018297776 
##        1712        1767        1780        1809        1827        1893 
## 0.007251875 0.006033315 0.129882566 0.010155854 0.006236266 0.006026162 
##        1920        1977        2016        2040        2087        2089 
## 0.007433404 0.007853849 0.006340711 0.018839571 0.008143447 0.008933365 
##        2114        2115        2120        2140        2149        2175 
## 0.038720587 0.022263050 0.068987431 0.006934197 0.021747942 0.050343017 
##        2200        2216        2224        2225        2257        2307 
## 0.023849200 0.017866605 0.006115518 0.005961769 0.010689685 0.025300710 
##        2391        2408        2437        2452        2458        2478 
## 0.011874776 0.013697537 0.093290977 0.093247923 0.017586843 0.010755526

#Outliers
plot(rstudent(mod8))
abline(h=c(-2,2), col=2)
outlierTest(mod8)

##       rstudent unadjusted p-value Bonferroni p
## 1551  6.067677         1.4956e-09   3.7390e-06
## 1306  4.902591         1.0066e-06   2.5165e-03
## 1399 -4.397823         1.1393e-05   2.8483e-02
## 155   4.391160         1.1745e-05   2.9363e-02
## 1694  4.282947         1.9142e-05   4.7856e-02

#Cook
cook <- cooks.distance(mod8)
which(cook > 1)

## 1551 
## 1551

plot(cook)
max(cook)

## [1] 1.704646

Il primo plot in alto a sinistra mostra le osservazioni (nello spazio dei regressori) che hanno un leverage superiore alla soglia calcolata. lev[lev>soglia] elenca le osservazioni che superano la soglia e come si può notare risultano essere tante. Potenzialmente queste osservazioni sono critiche perchè si trovano lontano nel dominio delle variabili esplicative.

Il secondo plot in alto a destra mostra i residui standardizzati del modello di regressione rispetto all’indice delle osservazioni. Serve, appunto, per individuare outlier o osservazioni influenti nel modello. Come si può notare, La maggior parte dei residui si trova tra -2 e +2, il che suggerisce che non ci sono gravi problemi di outlier nei residui. Alcuni punti più estremi si avvicinano a ±3 o oltre, ma non sembrano numerosi. Questo indica che ci potrebbero essere pochi outlier, ma non in una quantità tale da compromettere il modello. Nessun pattern evidente: Se ci fosse una struttura chiara nei residui (ad esempio una forma a onda), significherebbe che il modello potrebbe essere mal specificato. Qui, invece, i punti sembrano distribuiti in modo casuale. L’OutlierTest mostra le osservazioni che si discostano significativamente dal modello. L’osservazione 1551 ha un rstudent di 6.07, il che significa che è un residuo molto estremo. Il p-value non aggiustato è 1.4956e-09, estremamente basso. Anche il p-value corretto con Bonferroni (3.7390e-06) è molto piccolo, indicando che l’osservazione è un outlier significativo. Anche le osservazioni 1306, 1399, 155 e 1694 hanno valori di rstudent superiori a ±4, che sono soglie critiche per identificare outlier. Il p-value Bonferroni è inferiore a 0.05 in tutti i casi, quindi possiamo considerarli outlier significativi.

Possiamo concludere che le osservazioni da attenzionare sono la numero 1551 e la 155 che oltre ad essere significativi in temrini di outliers erano stati screennati anche come valori di leva.

Tuttavia, per verificare se incidono sui coefficienti è necessario la verifica della Cook’s Distance (grafico in basso a sinistra) dove possiamo intravedere che solamente una osservazione supera il valore soglia di 1.0, ed è la 1551 (che ha un valore di Cook pari a 1.70)

#Analisi dell'osservazione screennata da Cook
dati[1551, ]

##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551         35            1         0         38 4370       315    374
##      Tipo.parto Ospedale Sesso
## 1551        Nat     osp3     F

summary(dati$Peso)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2990    3300    3284    3620    4930

summary(dati$Lunghezza)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   500.0   494.7   510.0   565.0

summary(dati$Cranio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     235     330     340     340     350     390

Analizzando l’osservazione 1551 si evidenziano alcune caratteristiche che la rendono particolarmente anomala rispetto alle distribuzioni delle variabili nel dataset. In particolare:

Peso del neonato: 4370 g → Molto alto rispetto alla media 3284 g e al terzo quartile 3620 g. Si avvicina molto al valore massimo 4930 g.
Lunghezza: 315 mm → Inferiore alla media 494.7 mm e al terzo quartile 510 mm. È molto più vicina al valore minimo 310 mm.
Cranio: 374 mm → Sopra il terzo quartile 350 mm e vicino al massimo 390 mm.
Gestazione: 38 settimane → Perfettamente nella norma.

Anche realisticamente tale osservazioni può risultare un outlier in quanto quello che salta all’occhio è che il neonato ha un peso molto alto (4370 g), ma una lunghezza molto bassa (315 mm) rispetto alla distribuzione della variabile.

Per tale motivo si è deciso di rimuovere tale osservazioni e vedere se ci sono variazioni statistiche in miglioramento.

# Modello senza l'osservazione 1551
dati_senza_1551 <- dati[-1551, ]  # Rimuove la riga 1551
mod10 <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
                    I(Gestazione^2) + I(Lunghezza^2), data = dati_senza_1551)
summary(mod10)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2) + I(Lunghezza^2), data = dati_senza_1551)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1174.18  -184.34   -12.62   163.26  1307.99 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.836e+03  9.038e+02  -3.138  0.00172 ** 
## Gestazione       1.947e+02  6.629e+01   2.936  0.00335 ** 
## Lunghezza       -1.872e+01  4.541e+00  -4.123 3.87e-05 ***
## Cranio           1.025e+01  4.188e-01  24.465  < 2e-16 ***
## SessoM           7.475e+01  1.094e+01   6.835 1.02e-11 ***
## I(Gestazione^2) -2.078e+00  8.689e-01  -2.391  0.01686 *  
## I(Lunghezza^2)   3.030e-02  4.625e-03   6.553 6.83e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.1 on 2492 degrees of freedom
## Multiple R-squared:  0.7414, Adjusted R-squared:  0.7408 
## F-statistic:  1191 on 6 and 2492 DF,  p-value: < 2.2e-16

Come possiamo notare dall’output, l’R² Adjusted supera il 74% (rispetto al 73% del modello precedente con l’osservazione 1551) e il coefficiente della variabile Lunghezza passa da -31 a -18, smussando quella problematica che avevamo cercato di spiegare in fase di definizione del modello e dell’incoerenza del fatto che all’aumentare della lunghezza del bambino il coefficiente risultava troppo negativo. I regressori, inoltre, rimangono tutti statisticamente significativi.

Altre visualizzazioni grafiche

In questa sezione mostriamo due plot che ci aiutano a confermare l’andamento concavo/convesso della variabile Gestazione e Lunghezza

# Effetto della durata della gestazione sul peso
ggplot(dati, aes(x=Gestazione, y=Peso)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="lm", formula=y ~ poly(x, 2), col="red") +
  labs(title="Effetto della durata della gestazione sul peso",
       x="Settimane di Gestazione",
       y="Peso del neonato (g)")

# Effetto della lunghezza del neonato sul peso
ggplot(dati, aes(x=Lunghezza, y=Peso)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="lm", formula=y ~ poly(x, 2), col="red") +
  labs(title="Effetto della lunghezza del neonato sul peso",
       x="Lunghezza del neonato (mm)",
       y="Peso del neonato (g)")

Guardando il grafico, specialmente, sulla variabile Lunghezza alcuni dubbi rimangono sul segno dei coefficienti tra quello lineare e quadratico, però la scelta del modello è stata basata su un miglior trade off tra statistica e realtà. Per tale motivo, si è ritenuto idoneo aver fatto le analisi sul modello scelto. Sul grafico del peso la dinamica, invece, trova riscontro su quanto documentato in precedenza.

Predizione (Esempio)

In conclusione, viene fornito un esempio di predizione del peso di una neonata corripondente ad una mamma che è alla 40esima settimana di gravidanza, basandoci sul modello scelto (mod10)

#Definizione della nuova osservazione (out-of-sample)
Gestazione_Val <- 40
Lunghezza_Val <- mean(dati$Lunghezza, na.rm=TRUE)
Cranio_Val <- mean(dati$Cranio, na.rm=TRUE)


nuova_obs_1 <- data.frame(
  Gestazione = Gestazione_Val,
  Lunghezza = Lunghezza_Val,
  Cranio = Cranio_Val,
  Sesso = factor("F", levels = levels(dati$Sesso)), 
  `I(Gestazione^2)` = Gestazione_Val^2,  
  `I(Lunghezza^2)` = Lunghezza_Val^2  
)

peso_previsto_1 <- predict(mod10, newdata = nuova_obs_1)
peso_previsto_1

##        1 
## 3263.544

La predizione del peso sulla base del modello scelto per una mamma alla 40esima settimana che aspetta una bambina è di 3263.54, appena al di sotto della media del dataset.

Modello Statistico per Prevedere il Peso Neonatale

Davide Mariano Carlucci

2025-02-23