1. Raccolta dei Dati e Struttura del Dataset

Per costruire il modello predittivo, abbiamo raccolto 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).

#Importo il CSV


dati = read.csv("neonati.csv", 
                sep= ",",
                stringsAsFactors = T, # all'interno del dataset ci sono variabili qualitative che vanno lette come fattori, e non come semplici stringhe
                fileEncoding = "ISO-8859-1")

attach(dati)

2. Analisi e Modellizzazione

Analisi Preliminari

- Comprensione della distribuzione e identificazione di outlier o anomalie:

# Visiono il summury e calcolo il numero di campioni

summary(dati)
n = nrow(dati)

head(dati,5)


# Esploro le osservazioni con valori minimi e massimi per le diverse variabili quantitative
which.min(dati$Peso)
which.max(dati$Peso)

which.min(dati$Lunghezza)
which.max(dati$Lunghezza)

which.min(dati$Cranio) 
which.max(dati$Cranio)

which.min(dati$Anni.madre) 
which.max(dati$Anni.madre)

which.min(dati$Gestazione) 
which.max(dati$Gestazione)

# Anni.madre presenta alcuni valori errati
dati[dati$Anni.madre == 0, ]
dati[dati$Anni.madre < 18, ]

dati[dati$Gestazione < 28, ]
dati[dati$Gestazione > 40, ]

# verifico la presenza di valori NA nel dataset
any(is.na(dati))

# correggo i valori errati relativi alla variabile Anni.madre che ho notato nel dataset andando a sostituirli con il valore della media di Anni.madre.


dati[1380, "Anni.madre"] <- round(mean(Anni.madre),0)
dati[1152, "Anni.madre"] <- round(mean(Anni.madre),0)

- Prima valutazione sulle variabili:

Variabile	Tipo	Descrizione
Anni.madre	Quantitativa continua	Misura dell’età della madre in anni.
N.gravidanze	Quantitativa discreta	Numero di gravidanze avute dalla madre.
Fumatrici	Qualitativa dicotomica	Variabile dummy (0 = madre non fumatrice, 1 = madre fumatrice).
Gestazione	Quantitativa continua	Numero di settimane di gestazione.
Peso	Quantitativa continua	Peso alla nascita del neonato in grammi.
Lunghezza	Quantitativa continua	Lunghezza alla nascita del neonato in millimetri.
Diametro	Quantitativa continua	Diametro craniale alla nascita del neonato in millimetri.
Tipo.parto	Qualitativa nominale	Tipo di parto: naturale o cesareo.
Ospedale	Qualitativa nominale	Ospedale dove è avvenuta la nascita (1, 2, o 3).
Sesso	Qualitativa dicotomica	Sesso del neonato: Maschio (M) o Femmina (F).

# Valuto la simmetria e la curtosi della distribuzione

moments::skewness(Peso)

## [1] -0.6470308

moments::kurtosis(Peso)-3

## [1] 2.031532

- Analisi della distribuzione della variabile risporta:

skewness = -0.65, assimetrica negativa
kurtosis = 2.03, coefficiente positivo, è una distribuzione leptocurtica

# Effettuo lo Shapiro test per saggiare la normalità della variabile risposta, questo perchè la sua anormalità ricade spesso anche sui residui.

shapiro.test(Peso)

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

plot(density(Peso))

qqnorm(Peso)
qqline(Peso, col = "red")

- Ristultato Shapiro test sulla variabile risposta:

p-value < 2.2e-16, valore estremamente piccolo, (inferiore alla soglia di significatività del 0.05). Pertanto si rifiuta l’ipotesi di normalità, quindi la variabile peso non segue una distribuzione normale. Verificheremo possibili impatti e problemi di anormalità nei residui nei test dedicati.

Dal grafico possiamo notare come la distribuzione abbia qualche problema sulle code, questi valori molto probabilmente sono la causa per cui il test di normalità fallisce. Trascurando le code, tutto sommato la distribuzione sembra seguire una distribuzione normale.

- Matrice di correlazione:

# Verifico la matrice di correlazione, che mostra la correlazione a due a due fra le variabili. Per semplicità applico un arrotondamento alla seconda cifra decimale.


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(dati, upper.panel = panel.smooth, lower.panel = panel.cor)

- Analisi della matrice di correlazione:

Per prima cosa andiamo a controllare la relazione tra il peso (la nostra variabile risposta) e tutte le altre variabili esplicative.

Il coefficiente di correlazione lineare fra peso e lunghezza è alto e positivo, possiamo quindi affermare che le due variabili sono molto correlate, anche lo scatter plot mostra una nuvola allungata ed in pendenza.
Si evidenzia una correlazione positiva anche tra il peso e il diametro del cranio, anche in questo caso lo scatter plot mostra una nuvola allungata e in pendenza, non sembra però una relazone perfettamente lineare, in quanto sembra cambiare leggermente pendenza.
Si evidenzia una correlazione positiva anche fra il peso e le settimane di gestazione, anche in questo caso dallo scatterplot non sembra una relazone perfettamente lineare, si evidenzia una leggera curva.
Si evidenzia anche una correlazione positiva tra il peso e il sesso.
Si rileve una correlazione tra alcuni regressori, in particolare la lunghezza, il diametro del cranio e le settimane di gestazione sono correlate fra loro. Approfondiremo le relazioni con test specifici per evitare di incorrere in problematiche di multicolinearità.
Si rileva anche una leggera correlazione positiva fra gli anni della madre, il numero di gravidanze e le settimane di gestazione.
Ad una prima analisi l’essere una fumatrice, il tipo di parto e l’ospedale di nascita non sembrano avere nessun tipo di correlazione (sono indipendenti) con il peso o con gli altri regressori.
Le relazioni relative alle variabili qualitative non rendono bene con gli scatter plot, inoltre anche il valore della correlazione non ha molto significato in questi casi, di seguito andiamo a valutare le variabili qualitative con i boxplot.

- Saggiamo le seguenti ipotesi:

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

- Relazione tra tipo di parto e ospedali

# Indago la relazione fra il tipo di parto e l'ospedale di nascita per verificare se in alcuni ospedali si fanno più parti cesarei

# Creo una tabella di contingenza
tabella_contingenza = table(Tipo.parto, Ospedale)

# test chi-quadro
chisq.test(tabella_contingenza)

## 
##  Pearson's Chi-squared test
## 
## data:  tabella_contingenza
## X-squared = 1.0972, df = 2, p-value = 0.5778

- Risultato test chi-quadrato per la variabile tipo di parto e ospedale:

Il valore della statistica chi-quadrato è circa 1.1.
I gradi di libertà sono 2

Poiché il p-value (0.58) è maggiore di 0.05, non rifiutiamo l’ipotesi nulla. Non ci sono quindi prove statistiche di una relazione significativa tra il tipo di parto e l’ospedale di nascita, pertanto le due variabili risultano essere indipendenti fra loro.

- Relazione tra il campione e la popolazione

# Andiamo ora ad indagare se la media del peso e della lunghezza di questo campione di neonati sia significativamente uguale a quello della popolazione

peso_medio_popolazione = 3300
lunghezza_media_popolazione = 500

# T-test che confronta il peso del campione con quello della popolazione
t.test(Peso, mu = peso_medio_popolazione)

## 
##  One Sample t-test
## 
## data:  Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

# T-test che confronta la lunghezza del campione con quella della popolazione
t.test(Lunghezza, mu = lunghezza_media_popolazione)

## 
##  One Sample t-test
## 
## data:  Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

- Risultato test-t per valutare il peso medio del campione rispetto a quello della popolazione:

La statistica t è pari a -1.52. Il p-value è circa 0.13, maggiore della soglia di significatività del 5%, quindi non si rifiutare l’ipotesi nulla. La media osservata per il campione (3284.08g) risulta essere leggermente inferiore a quella della popolazione, ma non significativamente diversan (3300g). La differenza quindi non è statisticamente rilevante.

- Risultato test-t per valutare la lunghezza media del campione rispetto a quella della popolazione:

La statistica t è pari a -10.08, indicando una sensibile differenza tra la media del campione e quella della popolazione. Il p-value è estremamente basso (< 2.2e-16), molto inferiore alla soglia di significatività del 5%, pertanto rifiutiamo l’ipotesi nulla. La lunghezza media del campione (494.69mm) è significativamente minore rispetto a quella della popolazione (500 mm).

- Relazione tra le misure antropometriche e il sesso

par(mfrow = c(1,2)) 


boxplot(Peso~Sesso)

boxplot(Lunghezza~Sesso)

boxplot(Cranio~Sesso)


# Effettuo anche un t-test, con il quale miro a saggiare l'ugualizia tra medie per gruppi indipendenti
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(Lunghezza~Sesso)

## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.582, df = 2459.3, 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.929470  -7.876273
## sample estimates:
## mean in group F mean in group M 
##        489.7643        499.6672

t.test(Cranio~Sesso)

## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M 
##        337.6330        342.4486

- Risultati delle statistiche t-test:

Variabile	Statistica.t.test	p.value
Peso	t = -12.106	p-value < 2.2e-16
Lunghezza	t = -9.582	p-value < 2.2e-16
Cranio	t = -7.4102	p-value = 1.718e-13

- Considerazioni sulle misure antropometriche in relazione al sesso del neonato:

Tutti i t-test risultano negativi ed indicando che le misure antropometriche medie dei neonati femmina sono significativamente inferiori rispetto a quelle dei neonati maschi.

Tutti i p-value sono estremamente piccoli (praticamente 0), molto inferiore al livello di significatività alfa del 5%. Pertanto in tutti e tre i casi si rifiuta l’ipotesi nulla con altissima confidenza e possiamo affermare che le medie delle misure antropometriche sono significativamente diverse fra i sessi.

Possiamo quindi aspettarci che il sesso abbia un beta di regressione significativo.

Ha particolaremente senso mantenere questa varibile nel modello come variabile di controllo, così ci consentirà di effettuare predizzioni più accurata del peso del neonato in funzione del sesso.

par(mfrow = c(1,2)) 

boxplot(Peso, ylab = "Peso alla nascita (grammi)", main = "Distribuzione del Peso")

boxplot(Peso ~ Fumatrici,
        ylab = "Peso alla nascita (grammi)",        # Label asse Y
        xlab = "Fumatrici (0 = Non Fumatrici, 1 = Fumatrici)", # Label asse X
        main = "Peso rispetto Fumo della Madre")

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

- Considerazioni sul t-test relative al peso del neonato condizionato all’essere o meno una madre fumatrice:

Il test ha confermato quello che avevmo intuito dal risultato della matrice di correlazione e da una prima visualizzazione dello scatter plot.

p-value = 0.3033, è molto superiore al livello di significatività alfa del 5%, pertanto non si rifiuta l’ipotesi nulla, non c’è quindi un’evidenza statistica per affermare che il peso medio dei neonati differisca tra i figli di madri fumatrici e non fumatrici.

Creazione del Modello di Regressione

mod1 = lm(Peso ~. , data = dati)

summary(mod1)

## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1123.3  -181.2   -14.6   160.7  2612.6 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.1677   141.3977 -47.633  < 2e-16 ***
## Anni.madre        0.7983     1.1463   0.696   0.4862    
## N.gravidanze     11.4118     4.6665   2.445   0.0145 *  
## Fumatrici       -30.1567    27.5396  -1.095   0.2736    
## Gestazione       32.5265     3.8179   8.520  < 2e-16 ***
## Lunghezza        10.2951     0.3007  34.237  < 2e-16 ***
## Cranio           10.4725     0.4261  24.580  < 2e-16 ***
## Tipo.partoNat    29.5027    12.0848   2.441   0.0147 *  
## Ospedaleosp2    -11.2216    13.4388  -0.835   0.4038    
## Ospedaleosp3     28.0984    13.4972   2.082   0.0375 *  
## SessoM           77.5473    11.1779   6.938 5.07e-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.1 on 10 and 2489 DF,  p-value: < 2.2e-16

mod2 = update(mod1, ~.- Fumatrici)
summary(mod2) # togliendo la variabile Fumatrici l'R quadro aggiustato non cambia, era una variabile irrilevante

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1122.69  -181.85   -15.23   161.38  2615.71 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6734.8325   141.4031 -47.629  < 2e-16 ***
## Anni.madre        0.8083     1.1463   0.705   0.4808    
## N.gravidanze     11.1515     4.6606   2.393   0.0168 *  
## Gestazione       32.2721     3.8109   8.468  < 2e-16 ***
## Lunghezza        10.3092     0.3004  34.314  < 2e-16 ***
## Cranio           10.4768     0.4261  24.591  < 2e-16 ***
## Tipo.partoNat    29.2488    12.0830   2.421   0.0156 *  
## Ospedaleosp2    -11.1647    13.4392  -0.831   0.4062    
## Ospedaleosp3     28.3685    13.4955   2.102   0.0356 *  
## SessoM           77.3679    11.1772   6.922 5.65e-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.7287, Adjusted R-squared:  0.7278 
## F-statistic: 743.3 on 9 and 2490 DF,  p-value: < 2.2e-16

mod2b = update(mod2, ~.- Anni.madre)
summary(mod2b) # togliendo la variabile Anni.madre l'R quadro aggiustato non cambia, era una variabile irrilevante

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

mod3 = lm(Peso ~  Gestazione + Lunghezza + Cranio + Sesso, data = dati)

summary(mod3)

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

mod4 = lm(Peso ~  Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso, data = dati)
summary(mod4) # Inserendo la variabile N.gravidanze l'R quadro aggiustato aumenta

## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + 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 ***
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## 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

mod5 = lm(Peso ~  Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + Ospedale, data = dati)
summary(mod5) # Indago la variabile ospedale anche se ritengo possa essere ininfluente

## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + Ospedale, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1132.59  -183.63   -16.59   163.83  2620.70 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6682.0094   135.6710 -49.252  < 2e-16 ***
## Gestazione      32.0456     3.7933   8.448  < 2e-16 ***
## N.gravidanze    12.0824     4.3352   2.787  0.00536 ** 
## Lunghezza       10.2720     0.3003  34.204  < 2e-16 ***
## Cranio          10.5257     0.4256  24.729  < 2e-16 ***
## SessoM          77.4967    11.1865   6.928 5.42e-12 ***
## Ospedaleosp2   -11.0741    13.4494  -0.823  0.41036    
## Ospedaleosp3    29.1990    13.4998   2.163  0.03064 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.2 on 2492 degrees of freedom
## Multiple R-squared:  0.7281, Adjusted R-squared:  0.7273 
## F-statistic: 953.1 on 7 and 2492 DF,  p-value: < 2.2e-16

Selezione del Modello Ottimale

anova(mod4,mod3) # il test indica che  N.gravidanze ha una rilevanza statistica

## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F   Pr(>F)   
## 1   2494 188065546                                
## 2   2495 188688687 -1   -623141 8.2637 0.004079 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mod5,mod4) # il test indica che l'ospedale di nascita ha una rilevanza statistica

## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + 
##     Ospedale
## Model 2: Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F   Pr(>F)   
## 1   2492 187340680                                
## 2   2494 188065546 -2   -724866 4.8211 0.008133 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vif(mod3)

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

vif(mod4)

##   Gestazione N.gravidanze    Lunghezza       Cranio        Sesso 
##     1.669189     1.023475     2.074689     1.624465     1.040054

# Il VIF test serve a individuare e quantificare la multicollinearità. La presenza di multicollinearità tra i regressori può causare problemi nell'interpretazione del modello e nell'affidabilità delle stime.
# In entrambi i casi i valori sono inferiori a 5, pertanto non si presentano problemi multicollinearità.

BIC(mod5,mod4,mod3, mod1)

##      df      BIC
## mod5  9 35226.09
## mod4  7 35220.10
## mod3  6 35220.54
## mod1 12 35241.97

# Il Bayesian Information Criterion viene utilizzato per confrontare diversi modelli statistici che cercano di spiegare gli stessi dati. Nella sua stima tiene conto sia della bontà di adattamento ai dati osservati, sia della complessità del modello. 

# Il criterio indica che il miglior modello di regressione tra quelli proposti sia il mod4. Nonostante il BIC "premi" la semplicità del modello, il mod4 batte il mod3 e possiede il migliore compromesso in termini di adattamento e gradi di libertà.

- Considerazioni sul confronto tra modelli:

Nonostante il test anova indichi che l’ospedale di nascita porti un incremento ritenuto significativo per il modello, ho deciso di ignorare questa variabile in quanto logicamente possiamo considerarla ininfluente sulle misure antropometriche dei neonati.

La stessa considerazione vale per il tipo di parto effettuato, che possiamo anch’esso escludere logicamente.

Fatte queste considerazioni, il modello che prenderò in considerazione per i prossimi step di analisi sarà il mod4, che è stato anche individuato come migliore dal criterio BIC.

Analisi della Qualità del Modello

# I grafici diagnostici forniscono informazioni sulla distribuzione dei residui, l'eteroschedasticità e la normalità.
par(mfrow = c(2,2))
plot(mod4)

Considerazioni dei grafici diagnostici:

- Residuals vs Fitted:

Il Residuals vs Fitted mostra i residui del modello sull’asse y in funzione dei valori predetti (fitted) sull’asse x. La relazione potrebbe non essere perfettamente lineare. Si potrebbe valutare l’aggiunta di termini quadrati per le variabili predittive o applicare un GLM link log per provare a migliorare l’adattamento.

- Q-Q plot:

Il Q-Q plot valuta se un insieme di dati segue una distribuzione normale. Il grafico confronta i quantili dei residui standardizzati del modello con i quantili teorici di una distribuzione normale standard.

Complessivamente i punti sono distribuiti sulla bisettrice del grafico, però si possono notare delle deviazioni dalla normalità soprattutto nella coda destra. Questo indica che ci sono più valori estremi positivi di quanti ne preveda la normalità. Anche nella coda sinistra c’è un leggero discostamento, sebbene meno pronunciato.

- Scale-Location:

Ad un primo sguardo, la varianza nello grafico Scale-Location sembra relativamente costante. Tuttavia, ci sono alcuni aspetti da considerare che potrebbero spiegare perché il test Breusch-Pagan rileva eteroschedasticità.

La linea rossa nel grafico (che rappresenta la curva che meglio si adatta alla radice quadrata dei residui standardizzati rispetto ai valori predetti dal modello) non è perfettamente orizzontale, ma presenta una leggera curvatura all’estremità sinistra del grafico (valori predetti bassi).

Inoltre, anche se la dispersione verticale dei punti non mostra un andamento a imbuto evidente, si può comunque notare che la nuvola di punti tende ad “ispessirsi” leggermente man mano che i valori fitted aumentano. Questa curvatura, sebbene non molto pronunciata, suggerisce che possa esistere una relazione tra la varianza dei residui e i valori adattati.

- Residuals vs Leverage

La maggior parte delle osservazioni è concentrata verso il basso a sinistra (leverage basso e residui piccoli), suggerendo un modello relativamente stabile.

La linea rossa orizzontale mostra la media dei residui standardizzati, che dovrebbe essere vicina allo zero in un modello ben specificato. In questo caso specifico tende leggermente verso l’alto a destra (per leverage alti), indicando che i residui medi aumentano per osservazioni con leverage più elevato. Questo potrebbe succedere per la presenza di eteroschedasticità, oppure una leggera non linearità. Un altra possibile ragione potrebbe riguardare le ossevazioni influenti, ma tratto questo tema a breve andando a calcolare la distanza di Cook.

# Shapiro-Wilk, test statistico formale per saggiare la normalità:
shapiro.test(residuals(mod4))

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

plot(density(residuals(mod4)))

- Considerazioni dello Shapiro test sui residui del modello mod4

Rissultato: p-value < 2.2e-16: Un p-value molto piccolo indica un forte rifiuto dell’ipotesi nulla, pertanto i residui non seguono una distribuzione normale.

Nel grafico viene mostrata la distribuzione di densità dei residui, nonostante non superi lo shapiro test, i residui sembrano tutto sommato distribuirsi in maniera abbastanza normale. Probabilmente alcuni outlier presenti nella y non sono stati filtrati dal modello e contribuiscono a rendere non-normale anche la distribuzione dei residui.

- Misure statistiche:

# Test Breusch-Pagan per il rilevamento di eteroschedasticità:
lmtest::bptest(mod4)

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

- Considerazioni dello Breusch-Pagan test sui residui del modello mod4

Rissultato: p-value : 2.2e-16. Un p-value estremamente piccolo indica un forte rifiuto dell’ipotesi nulla. Ciò significa che la varianza degli errori non è costante, quindi vi è eteroschedasticità.

# Durbin-Watson, test statistico per verificare la presenza di autocorrelazione di primo ordine nei residui:

lmtest::dwtest(mod4)

## 
##  Durbin-Watson test
## 
## data:  mod4
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0

- Considerazioni dello Durbin-Watson test sui residui del modello mod4

Rissultato: p-value = 0.1224: non si rifiuta l’ipotesi nulla di incorrelazione, indicando quindi che i residui non presentano autocorrelazione significativa.

lev = hatvalues(mod4) 
plot(lev)

p = sum(lev) 
soglia = 2*  p/n
abline(h = soglia, col = 2)

#identifichiamo numericamente i valori di leva
lev[lev>soglia]

- Considerazioni sui valori di leva del modello mod4:

L’analisi ha identificato 138 valori di leva. Essi sono osservazioni che hanno un valore di x (variabile indipendente) insolitamente grande o piccolo rispetto alle altre osservazioni. Essi possono potenzialmente influenzare la retta di regressione, ma non necessariamente lo fanno.

plot(rstudent(mod4))
abline(h = c(-2,2),col = 2)

# identifico numericamente gli outlier
car::outlierTest(mod4)

##       rstudent unadjusted p-value Bonferroni p
## 1551 10.051908         2.4906e-23   6.2265e-20
## 155   5.027798         5.3138e-07   1.3285e-03
## 1306  4.827238         1.4681e-06   3.6702e-03

# Calcolo la distanza di Cook
cook = cooks.distance(mod4)
plot(cook)

- Considerazioni sui valori outlier del mod4:

Gli outlier sono osservaziono che hanno un valore di y (variabile dipendente) che si discosta in maniera molto significativa rispetto alle altre osservazioni.

Possono influenzare in modo sproporzionato la pendenza e l’intercetta della retta di regressione, “tirandola” verso di loro. Di conseguenza possono aumentare l’errore standard dei coefficienti di regressione e ridurre l’R quadro.

L’outlier test ha individuato tre osservazioni (1551, 155 e 1306) potenzialmente problematiche.

La distanza di Cook identifica i punti influenti, ovvero le osservazioni che causano un cambiamento sostanziale nei coefficienti di regressione stimati. Un punto influente può essere un outlier nella variabile dipendente, un punto di leva nelle variabili indipendenti, o una combinazione di entrambi.

La distanza di Cook riporta come unica osservazione problematica la 1551, con una distanza di circa 0.8. L’osservazione 1551 ha superato la soglia di allerta di 0.5 e sia avvicina pericolosamente al valore 1.

# controllo le righe associate alle osservazioni identificate come outlier per verificare eventuali errori o indizzi che possano suggerire le motivazioni di tali valori 
dati[1551,]
dati[155,]
dati[1306,]

# esploro altre osservazioni che presentano le medesime lunghezze, dimensioni del cranio e peso riscontrata negli outlier, così da poter avere una maggiore comprensione dei dati.

dati[dati$Lunghezza == 410, ]
dati[dati$Lunghezza == 315, ]
dati[dati$Lunghezza == 510, ]

dati[dati$Cranio == 374, ]
dati[dati$Peso == 4370, ]

dati[dati$Peso < 900, ]
dati[dati$Peso > 4470, ]

- Esplorazione degli outlier:

Comparando l’osservazione 928 con la 1551 si possono effettuar alcune considerazioni, in particolare l’osservazione 928 presenta la lunghezza minima del neonato, pari a 310mm, che corrispondente ad una femmina di 7 mesi.

Mi risulta molto strano che l’osservazione 1551 presenti un peso così elevato con una lunghezza così ridotta (simile a quela dell’osservazione 928), sopratutto considerando anche che le settimane di gestazione sono 40.

In ogni caso l’osservazione 1551 (la più estrama ) anche se è oltre la prima soglia di allerta 0.5, non supera la distanza di Cook. Ritengo sia una osservazione che valga la pena approfondire per verificare se la misurazione sia stata presa correttamente oppure se il neonato presenti particolari condizioni cliniche che abbiano sensonso di essere riportate e considerate nel modello, così da valutare se mantenere o meno la sua osservazione.

Al momento ho optato per mantenerla nel dataset non avendo ulteriori evidenze. Prendiamo però atto del fatto che questo modello non supera lo Shapiro test e il Breusch-Pagan test.

- Esploro la possibilità di utilizzare un modello GLM

mod4_glm = glm(Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso, 
               family = gaussian(link = "log"), 
               data = dati)
summary(mod4_glm)

## 
## Call:
## glm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + 
##     Cranio + Sesso, family = gaussian(link = "log"), data = dati)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7459811  0.0508941  93.252  < 2e-16 ***
## Gestazione   0.0139455  0.0012520  11.138  < 2e-16 ***
## N.gravidanze 0.0047062  0.0013021   3.614 0.000307 ***
## Lunghezza    0.0033880  0.0000946  35.814  < 2e-16 ***
## Cranio       0.0032547  0.0001315  24.755  < 2e-16 ***
## SessoM       0.0200878  0.0033897   5.926 3.53e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 75249.61)
## 
##     Null deviance: 688888542  on 2499  degrees of freedom
## Residual deviance: 187671317  on 2494  degrees of freedom
## AIC: 35174
## 
## Number of Fisher Scoring iterations: 4

coef(mod4_glm)

##  (Intercept)   Gestazione N.gravidanze    Lunghezza       Cranio       SessoM 
##  4.745981113  0.013945478  0.004706151  0.003387979  0.003254715  0.020087845

anova(mod4_glm)

## Analysis of Deviance Table
## 
## Model: gaussian, link: log
## 
## Response: Peso
## 
## Terms added sequentially (first to last)
## 
## 
##              Df  Deviance Resid. Df Resid. Dev
## NULL                           2499  688888542
## Gestazione    1 228071063      2498  460817479
## N.gravidanze  1   2457279      2497  458360200
## Lunghezza     1 221295945      2496  237064255
## Cranio        1  46747985      2495  190316270
## Sesso         1   2644953      2494  187671317

anova(mod4_glm, mod4)

## Analysis of Deviance Table
## 
## Model: gaussian, link: log
## 
## Response: Peso
## 
## Terms added sequentially (first to last)
## 
## 
##              Df  Deviance Resid. Df Resid. Dev
## NULL                           2499  688888542
## Gestazione    1 228071063      2498  460817479
## N.gravidanze  1   2457279      2497  458360200
## Lunghezza     1 221295945      2496  237064255
## Cranio        1  46747985      2495  190316270
## Sesso         1   2644953      2494  187671317

# calcolo lo pseudo R quadrato per verificare il grado di adattamento del modello stimato

pseudoR2<-function(mod) {1-(deviance(mod)/mod$null.deviance)}
pseudoR2(mod4_glm)

## [1] 0.7275738

BIC(mod4_glm,mod5,mod4,mod3, mod1)

##          df      BIC
## mod4_glm  7 35214.85
## mod5      9 35226.09
## mod4      7 35220.10
## mod3      6 35220.54
## mod1     12 35241.97

vif(mod4)

##   Gestazione N.gravidanze    Lunghezza       Cranio        Sesso 
##     1.669189     1.023475     2.074689     1.624465     1.040054

vif(mod4_glm)

##   Gestazione N.gravidanze    Lunghezza       Cranio        Sesso 
##     1.341959     1.022240     1.669524     1.413039     1.043609

- Considerazioni sul modello GML link log:

In base ai risultati ottenuti possiamo affermare che:

Dal test BIC in cui vengono comparati i diversi modelli di regressione, mod4_glm risulta il miglior compromesso tra bontà di adattamento ai dati osservati e complessità del modello.
Il test Anova mostra che il miglioramento introdotto da mod4_glm rispetto a mod4 (il secondo miglior modello) è statisticamente significativo.
L’analisi VIF effettua mostra come in entrambi i modelli non vi siano problemi di multicollinearità (per tutte le variabili si misurano valori sensibilmente inferiori a 5).

- Esploro la possibilità di introdurre degli effetti di interazioni nel modello di regressione:

# Provo ad indagare un effetto quadratico per quelle variabili che da scatterplot e matrice di correlazione risultavano correlate con una lieve curvatura. 

# In questo modello considero l'interazione sia tra gestazione e lunghezza, sia tra gestazione e diametro del caranio. particolare  
mod4_quad <- lm(Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + 
                 I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + 
                 I(Gestazione * Lunghezza) + I(Gestazione * Cranio), 
               data = dati)

summary(mod4_quad)

## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + 
##     I(Gestazione * Lunghezza) + I(Gestazione * Cranio), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1211.55  -181.09    -9.75   163.42  1322.10 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -3.986e+02  1.196e+03  -0.333 0.738916    
## Gestazione                 2.779e+02  7.204e+01   3.857 0.000118 ***
## N.gravidanze               1.458e+01  4.236e+00   3.442 0.000587 ***
## Lunghezza                 -1.823e+01  5.315e+00  -3.429 0.000616 ***
## Cranio                    -1.471e+01  1.011e+01  -1.455 0.145712    
## SessoM                     7.307e+01  1.098e+01   6.656 3.44e-11 ***
## I(Gestazione^2)           -1.761e+00  1.525e+00  -1.154 0.248439    
## I(Lunghezza^2)             5.357e-02  6.517e-03   8.220 3.25e-16 ***
## I(Cranio^2)                4.411e-03  1.654e-02   0.267 0.789691    
## I(Gestazione * Lunghezza) -6.082e-01  1.869e-01  -3.255 0.001149 ** 
## I(Gestazione * Cranio)     5.720e-01  1.855e-01   3.083 0.002071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.7 on 2489 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:   0.74 
## F-statistic: 712.2 on 10 and 2489 DF,  p-value: < 2.2e-16

# In questo modello considero l'interazione solo tra gestazione e lunghezza.
mod4_quad2 = lm(Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + 
                  I(Gestazione^2) + I(Cranio^2) + I(Gestazione * Cranio), 
                data = dati)

summary(mod4_quad2)

## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Cranio^2) + I(Gestazione * Cranio), 
##     data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1140.34  -182.14   -14.38   165.80  2643.30 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             735.19209 1200.73603   0.612 0.540405    
## Gestazione               98.69507   64.11334   1.539 0.123837    
## N.gravidanze             13.27434    4.29903   3.088 0.002039 ** 
## Lunghezza                10.42260    0.30162  34.555  < 2e-16 ***
## Cranio                  -42.04189    9.05792  -4.641 3.64e-06 ***
## SessoM                   72.90075   11.14040   6.544 7.25e-11 ***
## I(Gestazione^2)          -3.69610    1.02053  -3.622 0.000298 ***
## I(Cranio^2)               0.04003    0.01627   2.460 0.013961 *  
## I(Gestazione * Cranio)    0.66182    0.17053   3.881 0.000107 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 271.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7327, Adjusted R-squared:  0.7318 
## F-statistic: 853.5 on 8 and 2491 DF,  p-value: < 2.2e-16

# In questo modello considero l'interazione solo tra gestazione e diametro del cranio.
mod4_quad3 = lm(Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + 
                  I(Gestazione^2) + I(Lunghezza^2) + I(Gestazione * Lunghezza), 
                data = dati)

summary(mod4_quad3)

## 
## Call:
## lm(formula = Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Gestazione * 
##     Lunghezza), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1208.06  -180.96   -12.49   165.40  1327.65 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -2.682e+03  9.194e+02  -2.917  0.00356 ** 
## Gestazione                 3.067e+02  6.446e+01   4.758 2.07e-06 ***
## N.gravidanze               1.442e+01  4.245e+00   3.398  0.00069 ***
## Lunghezza                 -2.845e+01  4.448e+00  -6.396 1.89e-10 ***
## Cranio                     1.037e+01  4.208e-01  24.641  < 2e-16 ***
## SessoM                     7.358e+01  1.100e+01   6.688 2.79e-11 ***
## I(Gestazione^2)           -1.348e+00  1.524e+00  -0.885  0.37647    
## I(Lunghezza^2)             5.332e-02  6.396e-03   8.336  < 2e-16 ***
## I(Gestazione * Lunghezza) -3.377e-01  1.713e-01  -1.971  0.04886 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7396, Adjusted R-squared:  0.7388 
## F-statistic: 884.6 on 8 and 2491 DF,  p-value: < 2.2e-16

BIC(mod4_quad, mod4_quad2, mod4_quad3, mod4_glm)

##            df      BIC
## mod4_quad  12 35127.45
## mod4_quad2 10 35190.83
## mod4_quad3 10 35125.07
## mod4_glm    7 35214.85

anova(mod4_quad3, mod4_glm)

## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2) + I(Lunghezza^2) + I(Gestazione * Lunghezza)
## Model 2: Peso ~ Gestazione + N.gravidanze + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq     F    Pr(>F)    
## 1   2491 179359271                                 
## 2   2494 187671317 -3  -8312046 38.48 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

- Considerazioni sul modello con interazioni fra variabili:

Nonostante mod4_quad3 sia il modello valutato meglio dal criterio BIC , ritengo che il mod4_glm sia preferibile. Questo perchè presenta solamente 7 gradi di liberta rispetto ai 10 del modello con l’effetto quadratico, ma sopratutto perchè i GLM sono più semplici e migliorano l’interpretabilità del modello.

Grafici diagnostici e misure statisiche per mod4_GLM:

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

- Considerazioni generali sui grafici diagnostici del mod4_GLM:

Ad occhio il grafico Residuals vs Fitted mi sembra molto simile a quello relativo al modello lineare mod4, la differenza è che sembra “specchiato” sull’asse x, ma la linea rossa mi sembra molto simile.
Ad occhio anche il Q-Q plot mi sembra molto simile a quello del modello lineare mod4, in particolare sembra andare meglio sulla coda sinistra, ma probabilmente peggiora un po’ sulla coda destra.
Il grafico Scale-Location mi sembra leggermente più lineare, ma davvero in maniera impercettibile.
Il grafico Residuals vs Leverage migliora in quanto l’ossevazione 1551 scende sotto la soglia critica dello 0.5.

Guardando i grafici non percepisco un miglioramento sostanziale rispetto al modello lineare mod4.

shapiro.test(residuals(mod4_glm))

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

- Considerazioni dello Shapiro test sui residui del modello mod4_GLM:

Si rifiuta l’ipotesi nulla, pertanto i residui non seguono una distribuzione normale.

lmtest::bptest(mod4_glm)

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

- Considerazioni del Breusch Pagan test sui residui del modello mod4_GLM:

Si rifiuta l’ipotesi nulla, la varianza degli errori non è costante, vi è quindi eteroschedasticità.

lmtest::dwtest(mod4_glm)

## 
##  Durbin-Watson test
## 
## data:  mod4_glm
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0

- Considerazioni del Durbin Watson test sui residui del modello mod4_GLM:

Non si rifiuta l’ipotesi nulla di incorrelazione, quindi i residui non sono correlati.

lev = hatvalues(mod4_glm) 
plot(lev)

p = sum(lev)
soglia = 2*  p/n
abline(h = soglia, col = 2)

lev[lev>soglia]

- Risultato valori di leva del mod4_GLM:

Si possono identificare 156 valori leverage

plot(rstudent(mod4_glm))
abline(h = c(-2,2),col = 2)

car::outlierTest(mod4_glm)

##      rstudent unadjusted p-value Bonferroni p
## 1551 9.125365         7.1495e-20   1.7874e-16
## 155  4.887339         1.0221e-06   2.5552e-03
## 1306 4.730108         2.2440e-06   5.6100e-03

cook = cooks.distance(mod4_glm)
plot(cook)

- Considerazioni sui valori outlier del mod4_glm:

Il grafico mostra diversi valori che superano i valori soglia di -2 e 2.

Andando a misurare con il test degli outlier si evidenziano 3 osservazioni (1551, 155, 1306) potenzialmente problematiche.

Il valore massimo raggiunto per la desitanza di cook è di circa 0.25, il chè significa che gli outlier non dovrebbe essere un problema in quanto non supera la soglia di allerta dello 0.5. Non si evidenziano quindi valori particolarmente infuelnti sulle stime di regressione.

- Considerazioni generali sul mod4_glm

Prendiamo atto del fatto che anche questo modello non supera lo Shapiro test e il Breusch-Pagan test. Il test Durbin-Watson ci conferma che i residui non sono correlati. Dalla distanza di Cook non si rilevano valori problematici.

mod4 possiede un R quadrato di 0.727
mod4_glm possiede uno pseudo R quadrato di 0.7275

In fase di esplorazione non ho percepito un netto miglioramento delle performace portato dai modelli più complessi, per questa ragione ho optato di proseguire con la fase di predizione utilizzando il modello di regressione lineare semplice mod4.

Analisi della Qualità del Modello

# Creo le predizioni del modello
predizioni = predict(mod4, newdata = dati)

# Calcolo RMSE manualmente
RMSE = sqrt(mean((Peso - predizioni)^2)) 
print(paste("RMSE:", round(RMSE) ))

## [1] "RMSE: 274"

mean(Peso)

## [1] 3284.081

Considerazioni RMSE per il mod4:

Considerando la media del peso di 3284 g e un RMSE di 274 g, il modello sbaglia dell’8.34%, che può essere accettabile in contesti esplorativi, ma non per scopi clinici.

Il modello mira al prevedere con precisione il peso dei neonati alla nascita per migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale. Pertanto ritengo che l’errore sia significativamente non trascurabile e il modello non sia sufficientemente preciso per tale scopo.

3. Previsioni e Risultati

# Crea un dataframe con i dati per la predizione
prediction_data = data.frame(Gestazione = 39, 
                       N.gravidanze = 3, 
                       Lunghezza = 495, # lunghezza media
                       Cranio = 340, # diametro cranio medio
                       Sesso = "F") 

predict(mod4, newdata = prediction_data)

##        1 
## 3273.939

# La predizione restituisce un peso di circa 3274 g, che mi sembra in linea con le osservazioni del dataset.

# Calcolo una veloce misura di confronto approssimativa che restituisce il peso medio per nenonati femmine con 39 settimane di gestazione recuperato dal dataset, la quale mi restiuisce un valore medio di 3233.5 g.
mean(Peso[Sesso == "F" & Gestazione == 39])

## [1] 3233.531

Considerazioni sui valori predetti:

La predizione restituisce un peso di circa 3274 g, che mi sembra in linea con le osservazioni del dataset.

Calcolando una veloce misura di confronto che restituisce il peso medio per nenonati femmine con 39 settimane di gestazione recuperato dal dataset, questa restiuisce un valore medio di 3233.5 g, abbastanza in linea al valore predetto.

4. Visualizzazioni

# Grafico che mostra la relazione tra le variabili Peso e Gestazione
ggplot(data = dati)+
  geom_point(aes(x = Gestazione,
                 y = Peso,
                 col = Sesso),position = "jitter")+
  geom_smooth(aes(x = Gestazione,
                  y = Peso,
                  col = Sesso), se = F, method = "lm")+
 
  geom_smooth(aes(x = Gestazione,
                  y = Peso),col = "black", se = F, method = "lm")

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

# Grafico che mostra la relazione tra le variabili Peso e Fumatrici 
ggplot(data = dati)+
  geom_point(aes(x = Fumatrici,
                 y = Peso,
                 col = Sesso),position = "jitter")+
  geom_smooth(aes(x = Fumatrici,
                  y = Peso,
                  col = Sesso), se = F, method = "lm")+

  geom_smooth(aes(x = Fumatrici,
                  y = Peso),col = "black", se = F, method = "lm")

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Considerazioni sui grafici

Nel primo grafico possiamo notare come vi sia una forte relazione positiva tra il peso e le settimane di gestazione, condizionate al sesso, e come le rette di regressione descrivano coerentemente questa relazione con un forte pendenza. Queste considerazioni possono essere generalizzate anche per le altre due variabili antropometriche (lunghezza e diametro del cranio).

Nel secondo grafico possiamo notare come la variabile fumatrici non abbia alcuna relazione rispetto al peso del neoato. Le rette di regressione pertanto sono coerentemente orrizontali, mostrando l’indipendenza fre le varaibili.

Previsione del Peso Neonatale

Michele Bianchi

2024-12-11

1. Raccolta dei Dati e Struttura del Dataset

2. Analisi e Modellizzazione

Analisi Preliminari

- Comprensione della distribuzione e identificazione di outlier o anomalie:

- Prima valutazione sulle variabili:

- Analisi della distribuzione della variabile risporta:

- Ristultato Shapiro test sulla variabile risposta:

- Matrice di correlazione:

- Analisi della matrice di correlazione:

- Saggiamo le seguenti ipotesi:

- Relazione tra tipo di parto e ospedali

- Risultato test chi-quadrato per la variabile tipo di parto e ospedale:

- Relazione tra il campione e la popolazione

- Risultato test-t per valutare il peso medio del campione rispetto a quello della popolazione:

- Risultato test-t per valutare la lunghezza media del campione rispetto a quella della popolazione:

- Relazione tra le misure antropometriche e il sesso

- Risultati delle statistiche t-test:

- Considerazioni sulle misure antropometriche in relazione al sesso del neonato:

- Considerazioni sul t-test relative al peso del neonato condizionato all’essere o meno una madre fumatrice:

Creazione del Modello di Regressione

Selezione del Modello Ottimale

- Considerazioni sul confronto tra modelli:

Analisi della Qualità del Modello

Considerazioni dei grafici diagnostici:

- Residuals vs Fitted:

- Q-Q plot:

- Scale-Location:

- Residuals vs Leverage

- Considerazioni dello Shapiro test sui residui del modello mod4

- Misure statistiche:

- Considerazioni dello Breusch-Pagan test sui residui del modello mod4

- Considerazioni dello Durbin-Watson test sui residui del modello mod4

- Considerazioni sui valori di leva del modello mod4:

- Considerazioni sui valori outlier del mod4:

- Esplorazione degli outlier:

- Esploro la possibilità di utilizzare un modello GLM

- Considerazioni sul modello GML link log:

- Esploro la possibilità di introdurre degli effetti di interazioni nel modello di regressione:

- Considerazioni sul modello con interazioni fra variabili:

Grafici diagnostici e misure statisiche per mod4_GLM:

- Considerazioni generali sui grafici diagnostici del mod4_GLM:

- Considerazioni dello Shapiro test sui residui del modello mod4_GLM:

- Considerazioni del Breusch Pagan test sui residui del modello mod4_GLM:

- Considerazioni del Durbin Watson test sui residui del modello mod4_GLM:

- Risultato valori di leva del mod4_GLM:

- Considerazioni sui valori outlier del mod4_glm:

- Considerazioni generali sul mod4_glm

Analisi della Qualità del Modello

Considerazioni RMSE per il mod4:

3. Previsioni e Risultati

Considerazioni sui valori predetti:

4. Visualizzazioni

Considerazioni sui grafici