E.Ceci - Modello Previsione Peso Neonati

ANALISI PRELIMINARE

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

kable(summary(neonati))

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

Il summary mostra immediatamente un’anomalia nei dati, ovvero che il minimo della variabile “Anni.madre” è 0. Andando a vedere il dataset, emerge che sia il valore 1380 che il valore 1152 soffrono sicuramente di un errore di registrazione del dato, avendo come valori rispettivamente 0 e 1.

# Sviluppo la funzione per il calcolo degli indici
funzione_indici_forma <- function(x, nome_var) {
  data.frame(
    Variabile = nome_var,
    Skewness = round(skewness(x),2),
    Kurtosis = round(kurtosis(x)-3,2)
  )
}

# Seleziono le variabili d'interesse
vars <- neonati[c("Anni.madre", "Peso", "Lunghezza", "Cranio", "Gestazione")]

# Calcolo indici per ogni variabile ottenendo lista di vettori
results_list <- lapply(names(vars), function(n) funzione_indici_forma(vars[[n]], n))

# Converto la lista in tabella
kable(results_tab <- do.call(rbind, results_list))

Variabile	Skewness	Kurtosis
Anni.madre	0.04	0.38
Peso	-0.65	2.03
Lunghezza	-1.51	6.49
Cranio	-0.79	2.95
Gestazione	-2.07	8.26

La variabile risposta “Peso” è asimmetrica negativa anche se in maniera leggera, testimoniando una leggera prevalenza dei valori più elevati ed ha un coefficiente di kurtosi > 0 risultando quindi leptocurtica.

La variabile che più si avvicina ai valori di simmetria e curtosi della normale è “Anni.madre” risultando quasi perfettamente simmetrica e di poco leptocurtica

“Cranio”, mostra una asimmetria di -0,79, mentre “Lunghezza” ne mostra una più accentuata, quasi doppia rispetto a “Cranio”, con una forma leptocurtica comune ad entrambe (maggiore per “Lunghezza”).

#Grafici di densita

par(mfrow=c(2,3))

cont_vars = neonati[c("Anni.madre", "Peso", "Lunghezza", "Cranio")]

for (cont_var in names(cont_vars)) {
  x = neonati[[cont_var]]
  plot(density(x, na.rm = TRUE),
       main = paste("Densità di", cont_var),
       xlab = cont_var)
  
  
  abline(v = mean(x, na.rm = TRUE), col = 2)
}


#Distribuzione di frequenze
disc_vars <- neonati[c("N.gravidanze", "Gestazione")]

for (disc_var in names(disc_vars)) {
  x = neonati[[disc_var]]
  hist(x, na.rm = TRUE,
       main = paste("Distribuzione di frequenza", disc_var),
       xlab = disc_var,
       col = "lightblue",
       border = "blue")
  
  abline(v = mean(x, na.rm = TRUE), col = 2, lwd = 2)
}

L’asimmetria negativa di “Gestazione” (-2.07) è ovviamente all’origine delle asimmetrie negative delle variabili “Peso”, “Cranio”, “Lunghezza”: gestazione più lunghe, a parità di condizioni, produrranno individui più sviluppati, ed essendo maggiormente presenti nel dataset gestazioni più lunghe avremo anche valori più elevati per le misure antropometriche.

Verifica delle seguenti ipotesi con i test adatti.

In alcuni ospedali si fanno più parti cesarei

tabella.tipo_parto <- table(neonati$Ospedale, neonati$Tipo.parto)

df.tipo_parto <- as.data.frame(tabella.tipo_parto)
names(df.tipo_parto) <- c("Ospedale", "Tipo_Parto", "Frequenza")

kable(tabella.tipo_parto <- table(neonati$Ospedale, neonati$Tipo.parto))

	Ces	Nat
osp1	242	574
osp2	254	595
osp3	232	603

ggpubr::ggballoonplot(data=df.tipo_parto,
                      fill="Tipo_Parto")

Ciò che emerge è che l’Ospedale 2 è quello presso il quale si svolge il maggior numero di parti cesarei. Mentre l’ospedale 3 è quello in cui si svolge il maggior numero di parti naturali.

Vediamo se la differenza per i parti cesarei sia statisticamente significativa con un test del Chi-quadro.

H0 = il tipo di parto è indipendente dall’ospedale

H1 = il tipo di parto dipende dall’ospedale

test = chisq.test(tabella.tipo_parto)

kable(data.frame(
  X2 = round(test$statistic,2),
  df = test$parameter,
  p_value = round(test$p.value,2)
),
caption = "Test Chi-quadro per tipo di parto")

Test Chi-quadro per tipo di parto
	X2	df	p_value
X-squared	1.1	2	0.58

Il p-value di 0,58 mostra che le differenze nel tipo di parto fra gli ospedali non sono statisticamente significative. Quindi queste differenze sono semplicemente dovute alla casualità del campionamento,e di conseguenza l’ospedale 2 non è caratterizzato da nessun tratto specifico che lo porti a svolgere un numero mediamente maggiore di parti cesarei rispetto agli altri due ospedali.

La media del peso e della lunghezza di questo campione di neonati è significativamente uguale a quelle della popolazione

Per saggiare questa ipotesi utilizzerò un t-test. Come valori di peso e lunghezza della popolazione ho utilizzato dei dati reperiti presso il sito di un noto ospedale pediatrico italiano (https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012/).

Il valore trovato per la lunghezza da riferire alla popolazione è: 500 mm

Il sistema d’ipotesi per la variabile “Lunghezza” è:

H0: mu = 500

H1: mu != 500

t.test_l = t.test(Lunghezza,
               mu = 500,
               conf.level = 0.95,
               alternative = "two.sided")

t_table_l = data.frame(
  Statistica_t = round(t.test_l$statistic, 3),
  Media_campionaria = round(t.test_l$estimate, 2),
  Media_attesa = 500,
  P_value = round(t.test_l$p.value, 4),
  IC_inferiore = round(t.test_l$conf.int[1], 2),
  IC_superiore = round(t.test_l$conf.int[2], 2)
)

kable(t_table_l, caption = "Risultati del test t per la variabile Lunghezza")

Risultati del test t per la variabile Lunghezza
	Statistica_t	Media_campionaria	Media_attesa	P_value	IC_inferiore	IC_superiore
t	-10.084	494.69	500	0	493.66	495.72

Il risultato mostra che la lunghezza media del campione è significativamente diversa da quella della popolazione. Infatti il p-value estremamente piccolo e un intervallo di confidenza che non contiene il valore della media attesa, portano a rifiutare l’ipotesi nulla.

Il sistema d’ipotesi per la variabile “Peso” è:

H0: mu = 3500

H1: mu != 3500

t.test_p = t.test(Peso,
       mu=3500,
       conf.level = 0.95,
       alternative = "two.sided")

t_table_p = data.frame(
  Statistica_t = round(t.test_p$statistic, 3),
  Media_campionaria = round(t.test_p$estimate, 2),
  Media_attesa = 3500,
  P_value = round(t.test_p$p.value, 4),
  IC_inferiore = round(t.test_p$conf.int[1], 2),
  IC_superiore = round(t.test_p$conf.int[2], 2)
)

kable(t_table_p, caption = "Risultati del test t per la variabile Peso")

Risultati del test t per la variabile Peso
	Statistica_t	Media_campionaria	Media_attesa	P_value	IC_inferiore	IC_superiore
t	-20.562	3284.08	3500	0	3263.49	3304.67

Anche in questo caso il valore ridotto del p-value e il posizionamento dell’intervallo di confidenza spingono a rifiutare l’ipotesi nulla portando ad affermare che il peso medio del campione è significativamente diverso da quella della popolazione.

Le misure antropometriche sono significativamente diverse tra i due sessi

Boxplot per visualizzare graficamente le differenze e prospetto delle statistiche di sintesi

Variabile “Peso” condizionata alla variabile “Sesso”

boxplot(Peso~Sesso,
        col=c("pink", "lightblue"))

# summary peso-maschi
summary_p_m = summary(neonati$Peso[neonati$Sesso == "M"])

df_summary_p_m = as.data.frame(t(as.numeric(summary_p_m)))
colnames(df_summary_p_m)  = names(summary_p_m)

kable(df_summary_p_m,
      caption = "Statistiche descrittive del Peso (Maschi)",
      digits = 0)

Statistiche descrittive del Peso (Maschi)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
980	3150	3430	3408	3720	4810

# summary peso-femmine
summary_p_f = summary(neonati$Peso[neonati$Sesso == "F"])

df_summary_p_f = as.data.frame(t(as.numeric(summary_p_f)))
colnames(df_summary_p_f)  = names(summary_p_f)

kable(df_summary_p_f,
      caption = "Statistiche descrittive del Peso (Femmine)",
      digits = 0)

Statistiche descrittive del Peso (Femmine)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
830	2900	3160	3161	3470	4930

Variabile “Lunghezza” condizionata alla variabile “Sesso”

boxplot(Lunghezza~Sesso,
        col=c("pink", "lightblue"))

# summary lunghezza-maschi
summary_l_m = summary(neonati$Lunghezza[neonati$Sesso == "M"])

df_summary_l_m = as.data.frame(t(as.numeric(summary_l_m)))
colnames(df_summary_l_m)  = names(summary_l_m)

kable(df_summary_l_m,
      caption = "Statistiche descrittive della Lunghezza (Maschi)",
      digits = 0)

Statistiche descrittive della Lunghezza (Maschi)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
320	490	500	500	515	560

# summary lunghezza-femmine
summary_l_f = summary(neonati$Lunghezza[neonati$Sesso == "F"])

df_summary_l_f = as.data.frame(t(as.numeric(summary_l_f)))
colnames(df_summary_l_f)  = names(summary_l_f)

kable(df_summary_l_f,
      caption = "Statistiche descrittive della Lunghezza (Femmine)",
      digits = 0)

Statistiche descrittive della Lunghezza (Femmine)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
310	480	490	490	505	565

Variabile “Cranio” condizionata alla variabile “Sesso”

boxplot(Cranio~Sesso,
        col=c("pink", "lightblue"))

# summary cranio-maschi
summary_c_m = summary(neonati$Cranio[neonati$Sesso == "M"])

df_summary_c_m = as.data.frame(t(as.numeric(summary_c_m)))
colnames(df_summary_c_m)  = names(summary_c_m)

kable(df_summary_c_m,
      caption = "Statistiche descrittive del diametro del Cranio (Maschi)",
      digits = 0)

Statistiche descrittive del diametro del Cranio (Maschi)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
265	334	343	342	352	390

# summary cranio-femmine
summary_c_f = summary(neonati$Cranio[neonati$Sesso == "F"])

df_summary_c_f = as.data.frame(t(as.numeric(summary_c_f)))
colnames(df_summary_c_f)  = names(summary_c_f)

kable(df_summary_c_f,
      caption = "Statistiche descrittive del diametro del Cranio  (Femmine)",
      digits = 0)

Statistiche descrittive del diametro del Cranio (Femmine)
Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
235	330	340	338	348	390

In tutte le misure antropometriche i maschi mostrano valori maggiori; più evidenti rispetto al peso (3408 g di media contro 3161 g), e meno rispetto a lunghezza (500 mm di media contro 490 mm) e diametro del cranio (media di 342 mm contro 338 mm).

Di seguito verifichiamo se le differenze evidenti nei boxplot e rintracciabili nell medie mostrate siano significative.

L’H0 di questo test per tutte e tre le misure antropometriche è che la differenza in media tra i due sessi sia 0. Per verificare tale ipotesi utilizzeremo un t.test.

T.test per la variabile “Peso”

t.test_p_s = t.test(data=neonati,
                Peso~Sesso)

t.tab_p_s = data.frame(
  Statistica_t = round(t.test_p_s$statistic, 3),
  gradi_libertà = round(t.test_p_s$parameter, 1),
  p.value = round(t.test_p_s$p.value,3),
  IC_superiore = round(t.test_p_s$conf.int[1], 2),
  IC_inferiore = round(t.test_p_s$conf.int[2], 2),
  media_F = round(t.test_p_s$estimate[1], 3),
  media_M = round(t.test_p_s$estimate[2], 3)
)

kable(t.tab_p_s, caption = "Risultati del test t per la variabile Peso fra maschi e femmine")

Risultati del test t per la variabile Peso fra maschi e femmine
	Statistica_t	gradi_libertà	p.value	IC_superiore	IC_inferiore	media_F	media_M
t	-12.106	2490.7	0	-287.11	-207.06	3161.132	3408.215

Sia il p-value minore del livello di significatività che un intevallo di confidenza non contenente lo 0 previsto in H0, portano a rifutare l’ipotesi nulla

T.test per la variabile “Lunghezza”

t.test_l_s = t.test(data=neonati,
                Lunghezza~Sesso)

t.tab_l_s = data.frame(
  Statistica_t = round(t.test_l_s$statistic, 3),
  gradi_libertà = round(t.test_l_s$parameter, 1),
  p.value = round(t.test_l_s$p.value,3),
  IC_superiore = round(t.test_l_s$conf.int[1], 2),
  IC_inferiore = round(t.test_l_s$conf.int[2], 2),
  media_F = round(t.test_l_s$estimate[1], 3),
  media_M = round(t.test_l_s$estimate[2], 3)
)

kable(t.tab_l_s, caption = "Risultati del test t per la variabile Lunghezza fra maschi e femmine")

Risultati del test t per la variabile Lunghezza fra maschi e femmine
	Statistica_t	gradi_libertà	p.value	IC_superiore	IC_inferiore	media_F	media_M
t	-9.582	2459.3	0	-11.93	-7.88	489.764	499.667

Anche in questo caso sia il valore del p-value che la posizione dell’intervallo di confidenza portano a rifutare l’ipotesi nulla.

T.test per la variabile “Cranio”

t.test_c_s = t.test(data=neonati,
                Cranio~Sesso)

t.tab_c_s = data.frame(
  Statistica_t = round(t.test_c_s$statistic, 3),
  gradi_libertà = round(t.test_c_s$parameter, 1),
  p.value = round(t.test_c_s$p.value,3),
  IC_superiore = round(t.test_c_s$conf.int[1], 2),
  IC_inferiore = round(t.test_c_s$conf.int[2], 2),
  media_F = round(t.test_c_s$estimate[1], 3),
  media_M = round(t.test_c_s$estimate[2], 3)
)

kable(t.tab_c_s, caption = "Risultati del test t per la variabile Cranio fra maschi e femmine")

Risultati del test t per la variabile Cranio fra maschi e femmine
	Statistica_t	gradi_libertà	p.value	IC_superiore	IC_inferiore	media_F	media_M
t	-7.41	2491.4	0	-6.09	-3.54	337.633	342.449

Il risultato di questo test è il medesimo dei due precedenti, con un rifiuto dell’H0 dato un p-value ben inferiore a 0,05 ed un intervallo di confidenza che non contiene lo 0.

In generale, sia i p-value che gli intervalli di confidenza risultanti dai test delle misure antropometriche inducono a rifiutare l’H0. Le differenze fra sessi sono quindi statisticamente significative.

CREAZIONE DEL MODELLO DI REGRESSIONE

Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

Verifico la normalità della variabile risposta “Peso” con test di Shapiro.

shap.test.p = shapiro.test(Peso)

p_value = round(shap.test.p$p.value,3)

kable(data.frame("p-value" = p_value),
      caption = "Risultato dello Shapiro test")

Risultato dello Shapiro test
p.value
0

Il p-value porta a rifiutare l’ipotesi nulla di normalità della variabile “Peso”.

Valutiamo graficamente le relazioni fra le variabili

panel.cor  = function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr  = par("usr"); on.exit(par(usr))
  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(neonati[, !(names(neonati) %in% c("Sesso", "Ospedale", "Tipo.parto", "Fumatrici"))], upper.panel = panel.smooth, lower.panel = panel.cor)

Le correlazioni maggiori risultano fra peso e lunghezza, con un valore pari a 0.8, e fra peso e cranio (0.7). Risultano al di sopra di un valore di 0.5 le correlazioni fra: - gestazione e peso (0.59); - gestazione e lunghezza (0.62); - lunghezza e cranio (0.60).

Graficamente la relazione fra peso e lunghezza sembra essere perfettamente lineare, così come lo sembra quella fra peso e cranio, mentre gestazione-peso, gestazione-lunghezza, gestazione-cranio e lunghezza-cranio non sembrano essere delle semplici relazioni lineari ma sembra esserci una relazione quadratica.

ELABORAZIONE DEI MODELLI

Elaborazione modello di partenza

Per il modello di partenza escluderemo a priori le variabili “Ospedale” e “Tipo.parto”, essendo variabili che intuitivamente si prevede non avranno reale potere predittivo.

model1 = lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Sesso, data = neonati)

summary(model1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1161.56  -181.19   -15.75   163.70  2630.75 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6714.4109   141.1515 -47.569  < 2e-16 ***
## Anni.madre       0.9585     1.1347   0.845   0.3984    
## N.gravidanze    11.2756     4.6690   2.415   0.0158 *  
## Fumatrici      -30.2959    27.5971  -1.098   0.2724    
## Gestazione      32.9331     3.8267   8.606  < 2e-16 ***
## Lunghezza       10.2342     0.3009  34.009  < 2e-16 ***
## Cranio          10.5177     0.4268  24.642  < 2e-16 ***
## SessoM          78.0845    11.2039   6.969 4.06e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic:   949 on 7 and 2492 DF,  p-value: < 2.2e-16

Da questo modello di partenza le variabili più significative risultano: “Gestazione”, “Lunghezza”, “Cranio”, “Sesso”.

“Anni.madre” e “Fumatrici” risultano variabili non statisticamente significative, mentre la variabile “N.gravidanze” statisticamente significativa al 5%.

In termini di impatto, al netto della signficatività statistica, tutte le variabili hanno un impatto positivo sul peso finale del neonato, con “Gestazione” che ha l’impatto maggiore, eccetto “Fumatrici” che ha un impatto in termini assoluti quasi pari a quello di “Gestazione”.

Per arrivare a selezionare il modello finale si procederà in primis con l’esclusione delle variabili non significative e successivamente con una valutazione di eventuali andamenti quadratici di alcune variabili e di effetti d’interazione fra le stesse

L’R quadro corretto di questo modello è un buon valore, circa il 73% (0,7264), quindi oltre il 70% della variabilità di Y è spiegata dalle variabili prese in considerazione.

Eliminazione variabili non significative: Anni.madre, Fumatrici, N.gravidanze

model2 = update(model1,~.-Anni.madre - Fumatrici)

summary(model2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = neonati)
## 
## 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 potere esplicativo del modello è rimasto immutato come testimonia l’R quadro corretto prevedibilmente invariato (da 0.7264 a 0.7265) vista la non significatività delle variabili rimosse. Inoltre, i coefficienti hanno subito delle variazioni minime rispetto al modello precedente.

Si procederà adesso all’indagine di effetti quadratici

Indagine effetti quadratici per la variabile Gestazione

Poiché durante l’analisi grafica delle variabili è emersa una relazione non puramente lineare fra la variabile “Gestazione” e la variabile risposta, si è deciso di inserire nel modello tale variabile esplicativa al quadrato, mantenendo al contempo la sua versione di 1° grado.

Aggiunta variabile “Gestazione^2”

model3.g_g2 = update(model2,.~.+I(Gestazione^2))

summary(model3.g_g2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2), data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1144.11  -181.17   -13.16   165.73  2662.56 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4650.2268   898.1706  -5.177 2.43e-07 ***
## N.gravidanze       12.5660     4.3361   2.898  0.00379 ** 
## Gestazione        -81.0486    49.7127  -1.630  0.10316    
## Lunghezza          10.3534     0.3039  34.074  < 2e-16 ***
## Cranio             10.6363     0.4280  24.854  < 2e-16 ***
## SessoM             75.7900    11.2340   6.747 1.88e-11 ***
## I(Gestazione^2)     1.5136     0.6617   2.287  0.02226 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7269 
## F-statistic:  1110 on 6 and 2493 DF,  p-value: < 2.2e-16

Il nuovo modello mostra un incremento irrilevante del suo potere esplicativo (da 0.7265 a 0.7269) L’inserimento del termine quadratico di Gestazione, ha provocato la perdita di significatività e il cambio di segno del coefficiente lineare, mentre il termine quadratico risulta significativo al 5%. Questo suggerisce collinearità tra le due componenti (lineare e quadratica), e quindi entrambe catturano porzioni simili della variazione del peso. Alla luce di ciò risulta quindi più utile semplificare il modello, mantenendo solo il termine più informativo (“Gestazione^2).

Rimozione della variabile “Gestazione”

model4.g2 = update(model2,.~.-Gestazione+I(Gestazione^2))

summary(model4.g2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2), data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1146.87  -181.09   -15.12   165.40  2643.89 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -6.101e+03  1.235e+02 -49.380  < 2e-16 ***
## N.gravidanze     1.255e+01  4.338e+00   2.894  0.00383 ** 
## Lunghezza        1.026e+01  2.987e-01  34.358  < 2e-16 ***
## Cranio           1.056e+01  4.255e-01  24.816  < 2e-16 ***
## SessoM           7.733e+01  1.120e+01   6.906 6.32e-12 ***
## I(Gestazione^2)  4.379e-01  5.053e-02   8.667  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2494 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7267 
## F-statistic:  1330 on 5 and 2494 DF,  p-value: < 2.2e-16

L’eliminazione del termine lineare ha riportato la massima significatività al termine quadratico, con un R quadro corretto che si attesta a metà fra i due precedenti modelli (0.7259)

Elaborazione dei modelli con interazione

Interazione Gestazione-Lunghezza

model5.int.g_l = lm(Peso ~ Gestazione*Lunghezza+N.gravidanze+Cranio+Sesso)
summary(model5.int.g_l)

## 
## Call:
## lm(formula = Peso ~ Gestazione * Lunghezza + N.gravidanze + Cranio + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1133.50  -179.45   -11.74   168.77  2653.34 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.991e+03  9.202e+02  -2.163 0.030619 *  
## Gestazione           -9.396e+01  2.480e+01  -3.789 0.000155 ***
## Lunghezza            -8.111e-02  2.027e+00  -0.040 0.968082    
## N.gravidanze          1.304e+01  4.319e+00   3.020 0.002551 ** 
## Cranio                1.076e+01  4.262e-01  25.245  < 2e-16 ***
## SessoM                7.228e+01  1.120e+01   6.454 1.31e-10 ***
## Gestazione:Lunghezza  2.729e-01  5.295e-02   5.153 2.76e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2493 degrees of freedom
## Multiple R-squared:  0.7299, Adjusted R-squared:  0.7292 
## F-statistic:  1123 on 6 and 2493 DF,  p-value: < 2.2e-16

L’inserimento della variabile d’interazione ha aumentato, anche se di poco, l’R quadro corretto da 0.7265 a 0.7292, con la perdita di significatività per la variabile “Lunghezza”.

Interazione Gestazione-Cranio

model6.int.g_c = lm(Peso ~ Gestazione*Cranio+N.gravidanze+Lunghezza+Sesso)
summary(model6.int.g_c)

## 
## Call:
## lm(formula = Peso ~ Gestazione * Cranio + N.gravidanze + Lunghezza + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1137.16  -180.86   -12.44   167.43  2696.03 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -193.11866 1106.42819  -0.175  0.86145    
## Gestazione        -140.67784   29.52621  -4.765 2.00e-06 ***
## Cranio              -9.83681    3.47497  -2.831  0.00468 ** 
## N.gravidanze        13.14181    4.31191   3.048  0.00233 ** 
## Lunghezza           10.47042    0.30096  34.790  < 2e-16 ***
## SessoM              72.05680   11.17200   6.450 1.34e-10 ***
## Gestazione:Cranio    0.53339    0.09028   5.908 3.94e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.8 on 2493 degrees of freedom
## Multiple R-squared:  0.7308, Adjusted R-squared:  0.7301 
## F-statistic:  1128 on 6 and 2493 DF,  p-value: < 2.2e-16

L’interazione Gestazione-Cranio ha prodotto il modello con l’R quadro corretto più alto finora (0.7301) e a differenza del modello precedente, la variabile “Cranio” non ha perso significatività come avvenuto per “Lunghezza”.

Interazione Lunghezza-Cranio

model7.int.l_c = lm(Peso ~ Lunghezza*Cranio+N.gravidanze+Gestazione+Sesso)
summary(model7.int.l_c)

## 
## Call:
## lm(formula = Peso ~ Lunghezza * Cranio + N.gravidanze + Gestazione + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.74  -180.48   -13.62   165.82  2866.17 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.801e+03  1.018e+03  -1.770  0.07685 .  
## Lunghezza        -3.047e-01  2.202e+00  -0.138  0.88996    
## Cranio           -4.759e+00  3.191e+00  -1.491  0.13601    
## N.gravidanze      1.295e+01  4.321e+00   2.996  0.00276 ** 
## Gestazione        3.810e+01  3.964e+00   9.610  < 2e-16 ***
## SessoM            7.327e+01  1.119e+01   6.545 7.20e-11 ***
## Lunghezza:Cranio  3.158e-02  6.528e-03   4.838 1.39e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7295, Adjusted R-squared:  0.7289 
## F-statistic:  1121 on 6 and 2493 DF,  p-value: < 2.2e-16

Quest’interazione ha prodotto un abbassamento dell’R quadro corretto rispetto al modello precedente, ora a 0.7289. Inoltre, entrambe le variabili coinvolte nell’interazione hanno perso significatività.

SELEZIONE DEL MODELLO OTTIMALE

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

Calcolo del BIC

	df	BIC
model1	9	35234
model2	7	35220
model3.g_g2	8	35223
model4.g2	7	35218
model5.int.g_l	8	35201
model6.int.g_c	8	35193
model7.int.l_c	8	35205

In assoluto, il modello con il BIC minore è il modello 6 che prevede l’interazione fra le varibili “Gestazione” e “Cranio”. Fra i modelli puramente lineari, il modello con il BIC minore è quello con le sole variabili maggiormente significative, mentre per i modelli che includono una componente quadratica quello che comprende la sola variabile “Gestazione” al quadrato risulta avere il BIC minore.

Calcolo VIF

VIF del modello di partenza “model1”

vif.mod1 = vif(model1)

kable(round(vif.mod1),2)

	x
Anni.madre	1
N.gravidanze	1
Fumatrici	1
Gestazione	2
Lunghezza	2
Cranio	2
Sesso	1

VIF del secondo medello “model2” con componente quadratica (presente solo “Gravidanza” elevata al quadrato)

vif.mod4 = vif(model4.g2)

kable(round(vif.mod4),2)

	x
N.gravidanze	1
Lunghezza	2
Cranio	2
Sesso	1
I(Gestazione^2)	2

VIF quinto modello “model5.int.g_l” (interazione gestazione-lunghezza)

vif.mod5 = vif(model5.int.g_l, type = "predictor")

kable(vif.mod5)

	GVIF	Df	GVIF^(1/(2*Df))	Interacts With	Other Predictors
Gestazione	1.699079	3	1.092368	Lunghezza	N.gravidanze, Cranio, Sesso
Lunghezza	1.699079	3	1.092368	Gestazione	N.gravidanze, Cranio, Sesso
N.gravidanze	1.024145	1	1.012001	–	Gestazione, Lunghezza, Cranio, Sesso
Cranio	1.640863	1	1.280962	–	Gestazione, Lunghezza, N.gravidanze, Sesso
Sesso	1.050345	1	1.024863	–	Gestazione, Lunghezza, N.gravidanze, Cranio

VIF sesto modello “model6.int.g_c” (interazione gestazione-cranio)

vif.mod6 = vif(model6.int.g_c, type = "predictor")

kable(vif.mod6)

	GVIF	Df	GVIF^(1/(2*Df))	Interacts With	Other Predictors
Gestazione	2.135379	3	1.134782	Cranio	N.gravidanze, Lunghezza, Sesso
Cranio	2.135379	3	1.134782	Gestazione	N.gravidanze, Lunghezza, Sesso
N.gravidanze	1.024177	1	1.012016	–	Gestazione, Cranio, Lunghezza, Sesso
Lunghezza	2.107498	1	1.451722	–	Gestazione, Cranio, N.gravidanze, Sesso
Sesso	1.048535	1	1.023980	–	Gestazione, Cranio, N.gravidanze, Lunghezza

VIF settimo modello “model7.int.l_c” (interazione lunghezza-cranio)

vif.mod7 = vif(model7.int.l_c, type = "predictor")

kable(vif.mod6)

	GVIF	Df	GVIF^(1/(2*Df))	Interacts With	Other Predictors
Gestazione	2.135379	3	1.134782	Cranio	N.gravidanze, Lunghezza, Sesso
Cranio	2.135379	3	1.134782	Gestazione	N.gravidanze, Lunghezza, Sesso
N.gravidanze	1.024177	1	1.012016	–	Gestazione, Cranio, Lunghezza, Sesso
Lunghezza	2.107498	1	1.451722	–	Gestazione, Cranio, N.gravidanze, Sesso
Sesso	1.048535	1	1.023980	–	Gestazione, Cranio, N.gravidanze, Lunghezza

I VIF sono minori di 5 per tutte le variabili sia nei modelli puramente lineari, che in quelli con una componente quadratica o che prevedono interazione.

Scelta del modello

Il modello scelto è il modello 2 , “model2”, con le variabili: “N.gravidanze”, “Gestazione”, “Lunghezza”, “Cranio” e “Sesso”.

La scelta ricade su questa versione poiché pur essendo la più semplice, non contenendo componenti quadratiche o interazioni, ha un potere esplicativo non così inferiore rispetto a modelli più complessi come vedremo in seguito e un BIC poco più alto rispetto al modello 6 per il quale si riscontra il BIC minore (35220 contro 35193).

Si è ritenuto che il modello 2 rappresentasse il miglior compromesso fra potere predittivo e complessità (anche interpretativa).

ANALISI DELLA QUALITA’ DEL MODELLO

Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.

Analisi dell’R quadro corretto

Per i modelli di regressione lineare multipla è meglio considerare l’R quadro corretto poiché esso non aumenta con la sola aggiunta di regressori, ma può addiritura diminuire qualora il regressore aggiunto abbia scarso potere esplicativo o si unisca ad un numero già consistente di regressori.

tabella = matrix(data = c(0.7264,
                          0.7265,
                          0.7269,
                          0.7267,
                          0.7292,
                          0.7301,
                          0.7289),
                 nrow = 7,
                 ncol = 1
)

colnames(tabella) = "R^2 CORRETTO"

row.names(tabella) = c(
                      "model1",
                      "model2",
                      "model3.g_g2",
                      "model4.g2",
                      "model5.int.g_l",
                      "model6.int.g_c",
                      "model7.int.l_c")
                      
kable(tabella)

	R^2 CORRETTO
model1	0.7264
model2	0.7265
model3.g_g2	0.7269
model4.g2	0.7267
model5.int.g_l	0.7292
model6.int.g_c	0.7301
model7.int.l_c	0.7289

Come anticipato nella sezione precedente, il model2, pur essendo il meno complesso, ha un potere esplicativo di 0.7265, che è minore in maniera trascurabile rispetto a quello del modello più performante, ovvero il modello 6 con interazione Gestazione-Cranio e un R quadro corretto di 0.7301.

Il modello con la componente quadratica, “model4.g2” riporta un R quadro corretto leggermente più alto del model2 (0.7267) con un BIC leggermente inferiore (35218).

Analisi del RMSE

rmse = function(actual, predicted) {
  sqrt(mean((actual - predicted)^2))
}

kable(round(rmse(
  actual = neonati$Peso,
  predicted = predict(model3.g_g2, newdata = neonati)
),2))

x
273.99

Considerando che il peso medio risultante la dataset è di 3284,9 g, un errore medio di 274 g tra valori predetti e osservati, rappresenta un margine d’errore inferiore al 10% (8,3%), quindi più che accettabile.

Analisi dei residui

Analisi grafica

par(mfrow=c(2,2))

plot(model2)

resid_vals = residuals(model2)
dens = density(resid_vals)
plot(dens)

abline(v = mean(resid_vals), col = 2)

Residuals vs Fitted
La maggior parte dei punti è distribuita attorno alla media con un pattern leggermente ricurvo soprattutto in corrispondenza della coda sinistra curvo. Questo significa che il modello non riesce a catturare pienamente la struttura dei dari, perché la relazione puramente lineare non è sufficiente a riprodurre la reale relazione che intercorre fra i dati. Ci sono tre valori che che sono più lontanti dalla media, in particolare il valore 1551;
Q-Q residuals
I punti ricadono per la maggior parte sulla bisettrice eccetto per i valori 1551, 155 e 1306, che sono anche i punti più distanti nel Residuals vs Fitted:
Scale-Location
Il grafico per la valutazione della varianza mostra un pattern leggermente ricurvo, più accentuato sulla coda sinistra, con i 3 valori già segnalati che risultano più distanti. Questo potrebbe far dubitare sull’omoschedasticità del modello;
Residuals vs Leverage
Il valore 1551 supera la soglia di 0.5 della cook’s distance, segnalandosi come valore a cui prestare attezione, essendo anche molto vicino alla soglia di 1.
Distribuzione dei residui
Risultano non distribuiti normalmente soprattutto per via della coda di destra.

Analisi quantitativa

Shapiro Test

shapiro.test(residuals(model2))

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

Il risultato porta a rifiutare l’ipotesi nulla di normalità della distribuzione dei residui, come già anticipato dalla rappresentazione grafica.

Breusch-Pagan Test

bptest(model2)

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

Il risultato porta a rifiutare l’ipotesi nulla di omoschedasticità come già intuibile dal grafico Scale-Location

Darbin-Watsont Test

dwtest(model2)

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

I risultati portano a non rifiutare l’ipotesi nulla di indipendenza degli errori.

Analisi dei valori influenti

Leverage

lev = hatvalues(model2)
plot(lev)
p = sum(lev)
threshold.lev = 2 *p/n
abline(h=threshold.lev, col = 2)

I valori leverage risultano in gran numero. I più distanti probabilmente sono gli stessi valori risultati anomali anche negli altri test: 1551,155 e 1306. Il valore più distante sarà ragionevolmente il 1551.

Outliers

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

out.t = outlierTest(model2)
kable(out.t$rstudent)

	x
1551	10.051908
155	5.027798
1306	4.827238

Sono confermati i tre valori anomali già emersi precedentemente, con il valore #1551 come il più anomalo fra i 3.

Cook’s distance

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

max_index = which.max(cook)
max_value = cook[max_index]

cook_table = data.frame(
  "Osservazione" = max_index,
  "Cook's Distance" = round(max_value, 2)
)

kable(cook_table, caption = "Osservazione con massima Cook's Distance")

Osservazione con massima Cook’s Distance
	Osservazione	Cook.s.Distance
1551	1551	0.83

Il valore più preoccupante è il 1551. Questo valore ha superato la soglia di 0.5, trovandosi a metà fra quest’ultima e la soglia di 1, con una distanza di Cook pari a 0.77.

L’osservazione in questione è una neonata, con peso di 4370 g, un’altezza di 315 mm e un diametro craniale di 374 mm. Probabilmente il l’anomalia risiede in un peso che è superiore di quasi il 40% alla media femminile (3161 g) accompagnato da una lunghezza di 315 mm che è minore di oltre il 35% rispetto alla media femminile (di 490 mm), divario che aumenta se si considerano le sole femmine che si trovano in quel range di peso e che hanno una lunghezza di almeno 500 mm.

L’ipotesi più plausibile è un’errore di registrazione di uno dei due dati, perché se si guarda ai neonati con un peso maggiore o uguale 4000 g, non si scende al di sotto di una lunghezza di 450 mm (fra maschi e femmine); mentre per i neonati di lunghezza fra i 300 mm e i 400 mm, in un solo caso si superano i 2000g, con oscillazioni che variano fra gli 830 g e i 1780g.

PREVISIONI E RISULTATI

Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = neonati)
## 
## Coefficients:
##  (Intercept)  N.gravidanze    Gestazione     Lunghezza        Cranio  
##     -6681.14         12.47         32.33         10.25         10.54  
##       SessoM  
##        77.99

Tutti i coefficienti inseriti hanno un impatto marginale positivo, quindi, a parità di altre variabili, in media ogni ulteriore gravidanza incide per 12 grammi in più sul peso finale, ogni settimana di gestazione aggiuntiva incide per circa 32 grammi (coefficiente più alto tra le variabili quantitative), mentre ogni millimetro in più di lunghezza del neonato e di diametro craniale incidono di circa 10 grammi ognuno. Infine, rispetto alla variabile qualitativa “Sesso”, in media i maschi pesano circa 78 g in più rispetto alle femmine (modalità baseline)

Effettuo la previsione con il modello 2: model2

Utilizzo come valori per le altre variabili le medie ottenute dal dataset:

Lunghezza = 495 Cranio = 340

kable(predict(model2,newdata = data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = mean(neonati$Lunghezza),
  Cranio = mean(neonati$Cranio),
  Sesso = "F")))

x
3271.09

Il valore predetto per una femmina che verrà partorita da una donna alla terza gravidanza dopo 39 settimane di gestazione è pari a 3271.09g

E.Ceci - Modello Previsione Peso Neonati

Enrico Ceci

2025-10-28