Importo il dataset csv e controllo la eventuale presenza di dati mancanti.

dati <- read.csv("neonati.csv", stringsAsFactors = T, sep = ",")
colSums(is.na(dati))

##   Anni.madre N.gravidanze    Fumatrici   Gestazione         Peso    Lunghezza 
##            0            0            0            0            0            0 
##       Cranio   Tipo.parto     Ospedale        Sesso 
##            0            0            0            0

nrow(dati)

## [1] 2500

1) Raccolta dat 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. (Variabile quantitativa discreta)
Numero di gravidanze: Quante gravidanze ha avuto la madre. (Variabile quantitativa discreta)
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice). (Variabile quantitativa discreta da trattare come indicatore qualitativo binario)
Durata della gravidanza: Numero di settimane di gestazione. (Variabile quantitativa discreta)
Peso del neonato: Peso alla nascita in grammi. (Variabile quantitativa continua)
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie. (Variabili quantitativa continue)
Tipo di parto: Naturale o cesareo. (Variabile qualitativa su scala nominale)
Ospedale di nascita: Ospedale 1, 2 o 3. (Variabile qualitativa su scala nominale)
Sesso del neonato: Maschio (M) o femmina (F). (Variabile qualitativa su scala nominale)

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.

2) Analisi e modellizzazione

Analisi preliminare

Tabella indicatori statistici delle variabili
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
mean	28.16	0.98	38.98	3284.08	494.69	340.03
Q.0%	0.00	0.00	25.00	830.00	310.00	235.00
Q.25%	25.00	0.00	38.00	2990.00	480.00	330.00
Q.50%	28.00	1.00	39.00	3300.00	500.00	340.00
Q.75%	32.00	1.00	40.00	3620.00	510.00	350.00
Q.100%	46.00	12.00	43.00	4930.00	565.00	390.00
sd	5.27	1.28	1.87	525.04	26.32	16.43
cv	18.72	130.51	4.79	15.99	5.32	4.83
range	46.00	12.00	18.00	4100.00	255.00	155.00
min	0.00	0.00	25.00	830.00	310.00	235.00
max	46.00	12.00	43.00	4930.00	565.00	390.00
skewness	0.04	2.51	-2.07	-0.65	-1.51	-0.79
curtosi	3.38	13.99	11.26	5.03	9.49	5.95

Nota: nonostante nel testo sia scritto esplicitamente che le misure del cranio riguardano il diametro e non sia specificata l’unità di misura, suppongo vi sia un errore visto che il raffronto con le tabulazioni in letteratura delle dimensioni in mm del diametro craniale sono molto minori, mentre quelle relative alla circonferenza risultano più compatibili con i dati presenti nel nostro campione. Lavoreremo su questi dati trattandoli come circonferenza del cranio

E’ necessaria una analisi preventiva del dataset per verificare la presenza di alcune misurazioni estreme anomale frutto di errori. Sicuramente si vede subito un dato anomalo legato all’età della madre che ha come minimo 0 anni. Eliminiamo i dati che hanno età troppo bassa e proseguiamo a controllare il resto.

dati <- subset(dati, Anni.madre > 10)

Verifico graficamente se gli altri dati mostrano evidenti misurazioni estreme.

attach(dati)
par(mfrow=c(1,2))
plot(Anni.madre, pch=20, col = "LightBlue")
plot(N.gravidanze, pch=20, col = "LightBlue")

par(mfrow=c(1,2))
plot(Gestazione, pch=20, col = "LightBlue")
plot(Peso, pch=20, col = "LightBlue")

par(mfrow=c(1,2))
plot(Lunghezza, pch=20,col = "LightBlue")
plot(Cranio, pch=20, col = "LightBlue")

I grafici mostrano alcuni dati lontani dal valore medio ma non escludibili a priori come errore di misurazione.

Possiamo verificare se ci sono dati anomali anche nella correlazione tra le diverse variabili.

par(mfrow=c(2,2))
plot(Peso, Lunghezza, pch=20, col = "LightBlue")
plot(Peso, Cranio, pch=20, col = "LightBlue")
plot(Peso, Gestazione, pch=20, col = "LightBlue")
plot(Anni.madre, N.gravidanze, pch=20, col = "LightBlue")

C’è un dato evidentemente anomalo nel rapporto peso-lunghezza, che mostra un peso alto nonostante una lunghezza bassa. Verifichiamo questo dato:

which(Peso>4000 & Lunghezza<350)

## [1] 1549

dati <- dati[-1549, ]

Per il resto i grafici incrociati mostrano dei possibili dati anomali soprattutto nel rapporto Peso del neondato-Dimensione del cranio, che possono però evidenziare delle anomalie nel neonato più che definire un evidente errore di misurazione.

Come ultima analisi descrittiva dei dati possiamo visualizzare la distribuzione delle variabili quantitative.

par(mfrow=c(1,2))
hist(Anni.madre, ylim = c(0,500), main = "", xlab = "Età della madre", ylab = "Frequenza", col = "LightBlue")
hist(N.gravidanze, ylim = c(0,2000), main = "", xlab = "Numero gravidanze precedenti", ylab = "Frequenza", col = "LightBlue")

par(mfrow=c(1,2))
hist(Gestazione,ylim = c(0,700), xlim = c(25, 45),  main = "", xlab = "Settimane di gestazione", ylab = "Frequenza", col = "LightBlue")
hist(Peso,  main = "", xlab = "Peso alla nascita", ylab = "Frequenza", col = "LightBlue")

par(mfrow=c(1,2))
hist(Lunghezza, ylim = c(0,1000), main = "", xlab = "Lunghezza del neonato", ylab = "Frequenza",col = "LightBlue")
hist(Cranio, ylim = c(0,700), xlim = c(220,400),  main = "", xlab = "Dimensione del cranio", ylab = "Frequenza",col = "LightBlue")

Nalle curve di distribuzione sono resi evidenti i dati di curtosi e skewness tabulati nella tabella precedente: - la distribuzione della età della madre è estremamente simmetrica. La distribuzione relativa al numero di gravidanze è invece asimmetrica positiva, mentre tutti le altre distribuzioni mostrano una asimmetria negativa. Gli inidici di curtosi positivi per tutte le variabili misurano delle distribuzioni leptocurtiche.

Studio anche le variabili qualitative.

attach(dati)
tab1 <- as.data.frame(table(Fumatrici))
tab1$Fumatrici <- c("Non fumatrici", "Fumatrici")
tab1$freq_rel <- round(table(Fumatrici)/length(Fumatrici)*100,2)

tab2 <- as.data.frame(table(Tipo.parto))
tab2$freq_rel <- round(table(Tipo.parto)/length(Tipo.parto)*100,2)

tab3 <- as.data.frame(table(Ospedale))
tab3$freq_rel <- round(table(Ospedale)/length(Ospedale)*100,2)

tab4 <- as.data.frame(table(Sesso))
tab4$freq_rel <- round(table(Sesso)/length(Sesso)*100,2)

Fumo durante la gestazione
Categoria	Frequenza	Frequenza relativa
Non fumatrici	2393	95.84
Fumatrici	104	4.16

Parto cesareo o naturale
Categoria	Frequenza	Frequenza relativa
Ces	728	29.15
Nat	1769	70.85

Ospedale di provenienza
Categoria	Frequenza	Frequenza relativa
osp1	816	32.68
osp2	848	33.96
osp3	833	33.36

Sesso del neonato
Categoria	Frequenza	Frequenza relativa
F	1254	50.22
M	1243	49.78

Per le variabili qualitative possiamo vedere una distribuzione di valori molto dissimile per le variabili Fumatrici e Tipo di parto, mentre un forte equilibrio nelle variabili Ospedale e Sesso.

In ultimo ricompilo la tabella degli indici di posizione del dataframe da cui sono stati tolti i valori estremi

Tabella indicatori statistici delle variabili
	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio
mean	28.18	0.98	0.04	38.98	3283.75	494.77	340.02
Q.0%	13.00	0.00	0.00	25.00	830.00	310.00	235.00
Q.25%	25.00	0.00	0.00	38.00	2990.00	480.00	330.00
Q.50%	28.00	1.00	0.00	39.00	3300.00	500.00	340.00
Q.75%	32.00	1.00	0.00	40.00	3620.00	510.00	350.00
Q.100%	46.00	12.00	1.00	43.00	4930.00	565.00	390.00
sd	5.22	1.28	0.20	1.87	524.88	26.09	16.42
cv	18.51	130.53	479.78	4.80	15.98	5.27	4.83
range	33.00	12.00	1.00	18.00	4100.00	255.00	155.00
min	13.00	0.00	0.00	25.00	830.00	310.00	235.00
max	46.00	12.00	1.00	43.00	4930.00	565.00	390.00
skewness	0.15	2.51	4.59	-2.07	-0.65	-1.43	-0.79
curtosi	2.90	13.98	22.05	11.26	5.03	8.95	5.95

Test sui parti cesarei per ospedale

verifica_parto <-round(prop.table(table(Ospedale, Tipo.parto), margin = 1)*100,2)

Distribuzione percentuale dei tipi di parto per ospedale
Ospedale	Cesareo	Naturale
osp1	29.66	70.34
osp2	29.95	70.05
osp3	27.85	72.15

Dalla tabella precedenti non si evidenziano particolari differenze tra i tre ospedali, se non una lievissima variazione nell’osp3. Possiamo utilizzare il test del Chi Quadro

chisq.test(verifica_parto)

## 
##  Pearson's Chi-squared test
## 
## data:  verifica_parto
## X-squared = 0.1254, df = 2, p-value = 0.9392

Il p value del Test di 0.93 è di molto maggiore della soglia di accettazione di 0.05 per cui non si può affermare che ci sia una correlazione tra il numero di parti cesarei e l’ospedale in cui viene e effettuato il parto.

Analisi peso e lunghezza del campione

Ci proponiamo di verificare se la media del peso e della lunghezza di questo campione di neonati siano significativamente uguali a quelle della popolazione.

https://www.who.int/tools/child-growth-standards/standards/weight-for-age

Il link precedente rimanda al sito della World Health Organization da cui abbiamo estratto le tabulazioni riguardanti peso e lunghezza dei neonati in ogni settimana del primo anno di vita.

Le seguenti due sono le tabelle relative al peso.

Effettuiamo il test T di student per verificare l’ipotesi di uguaglianza tra il campione e la popolazione.

t_test_F_peso = t.test(Peso[dati$Sesso == "F"],
       mu = 3200,
       conf.level = 0.95,
       alternative = "two.sided")


t_test_M_peso = t.test(Peso[dati$Sesso == "M"],
       mu = 3300,
       conf.level = 0.95,
       alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  Peso[dati$Sesso == "F"]
## t = -2.6883, df = 1253, p-value = 0.007276
## alternative hypothesis: true mean is not equal to 3200
## 95 percent confidence interval:
##  3130.978 3189.217
## sample estimates:
## mean of x 
##  3160.097

## 
##  One Sample t-test
## 
## data:  Peso[dati$Sesso == "M"]
## t = 7.7447, df = 1242, p-value = 1.978e-14
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3381.012 3435.979
## sample estimates:
## mean of x 
##  3408.496

Le seguenti due tabelle sono relative alla lunghezza.

Test t di student per i valori di lunghezza:

t_test_F_lunghezza = t.test(Lunghezza[dati$Sesso == "F"],
       mu = 491,
       conf.level = 0.95,
       alternative = "two.sided")


t_test_M_lunghezza = t.test(Lunghezza[dati$Sesso == "M"],
       mu = 499,
       conf.level = 0.95,
       alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  Lunghezza[dati$Sesso == "F"]
## t = -1.4323, df = 1253, p-value = 0.1523
## alternative hypothesis: true mean is not equal to 491
## 95 percent confidence interval:
##  488.4016 491.4054
## sample estimates:
## mean of x 
##  489.9035

## 
##  One Sample t-test
## 
## data:  Lunghezza[dati$Sesso == "M"]
## t = 0.98965, df = 1242, p-value = 0.3225
## alternative hypothesis: true mean is not equal to 499
## 95 percent confidence interval:
##  498.3369 501.0131
## sample estimates:
## mean of x 
##   499.675

I Test t di Student effettuati per peso e lunghezza sui campioni di neonati maschi e femmine mostrano tutti e quattro dei valori di p value molto bassi, che quindi portano in tutti e quattro i casi a rifiutare l’ipotesi nulla di somiglianza tra campione e popolazione. Viene invece accettata l’ipotesi alternativa di differenza significatica del campione dai valori della popolazione neonatale complessiva.

Analisi misure antropometriche

Ci proponiamo in ultimo di verificare se le misure antropometriche sono significativamente diverse tra i due sessi. Analizziamo perciò peso, lunghezza e dimensione del cranio in rapporto al sesso del neonato.

par(mfrow=c(1,3))
boxplot(Peso~Sesso, col = "LightBlue", main = "Peso")
boxplot(Lunghezza~Sesso, col = "LightBlue", main = "Lunghezza")
boxplot(Cranio~Sesso, col = "LightBlue", main = "Dimensione del cranio", ylab = "Cranio")

I boxplot mostrano delle distribuzioni di dati differenti per i maschi e per le femmine. Possiamo visualizzare questo anche in una tabulazione dei valori medi e verificare l’ipotesi con il test t di Student.

tab_antrop <- dati %>%
    group_by(Sesso) %>%
    summarize(
      peso_medio = round(mean(Peso),2),
      lunghezza_media= round(mean(Lunghezza),2),
      dimensione_cranio_medio = round(mean(Cranio),2)
    )

Tabella valori antropometrici medi
Sesso	peso_medio	lunghezza_media	dimensione_cranio_medio
F	3160.1	489.90	337.59
M	3408.5	499.67	342.46

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.17, df = 2487.9, 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:
##  -288.4204 -208.3762
## sample estimates:
## mean in group F mean in group M 
##        3160.097        3408.496

## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.5303, df = 2464.9, 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.782029  -7.760913
## sample estimates:
## mean in group F mean in group M 
##        489.9035        499.6750

## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4855, df = 2488.6, p-value = 9.835e-14
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.138779 -3.590159
## sample estimates:
## mean in group F mean in group M 
##        337.5941        342.4586

Tutti e tre i t test mostrano un p value molto basso e ci portano a selezionare per tutte e tre le variabili l’ipotesi alternativa di differenza significativa tra le misure antropometriche dei neonati dei due sessi.

Creazione modello di regressione

Ci proponiamo di trovare un modello di regressione lineare multipla che includa le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni.

Facciamo innanzitutto un test sulla distribuzione della variabile di output Peso per verificare cha abbia distribuzione normale e non produca effetti significativi nel calcolo dei residui.

Indicatori statistici della variabile peso
	Peso
mean	3283.75
Q.0%	830.00
Q.25%	2990.00
Q.50%	3300.00
Q.75%	3620.00
Q.100%	4930.00
sd	524.88
cv	15.98
range	4100.00
min	830.00
max	4930.00
skewness	-0.65
curtosi	5.03

Gli indicatori statistici del Peso ci mostrano una distribuzione leptocurtica con leggera asimmetria negativa.

shapiro.test(Peso)

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

I p value del test di Shapiro Wilk è molto basso, e quindi ci porta ad accettare la distribuzione del Peso come normale.

Parte deterministica del modello Iniziamo visualizzando la matrice di correlazione di tutte le variabili:

Possiamo sviluppare il modello di regressione valutando i coefficienti di regressione Beta, che rappresentano gli effetti marginali di ogni singola variabile sulla variabile Peso.

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
##     Fumatrici + Lunghezza + Cranio + Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1173.53  -179.25   -13.31   162.38  1407.71 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6704.9105   138.5985 -48.377  < 2e-16 ***
## Anni.madre       0.6050     1.1272   0.537  0.59153    
## N.gravidanze    12.4548     4.5873   2.715  0.00667 ** 
## Gestazione      30.0838     3.7655   7.989 2.05e-15 ***
## Fumatrici      -27.3365    27.0750  -1.010  0.31276    
## Lunghezza       10.8757     0.3022  35.988  < 2e-16 ***
## Cranio           9.9053     0.4233  23.402  < 2e-16 ***
## SessoM          78.2423    10.9960   7.116 1.45e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.4 on 2489 degrees of freedom
## Multiple R-squared:  0.7373, Adjusted R-squared:  0.7366 
## F-statistic: 998.1 on 7 and 2489 DF,  p-value: < 2.2e-16

Coefficienti del modello lineare
term	estimate	std.error	statistic	p.value
(Intercept)	-6704.91049	138.59849	-48.37650	0.00000
Anni.madre	0.60495	1.12719	0.53669	0.59153
N.gravidanze	12.45484	4.58734	2.71505	0.00667
Gestazione	30.08384	3.76550	7.98933	0.00000
Fumatrici	-27.33646	27.07504	-1.00966	0.31276
Lunghezza	10.87569	0.30220	35.98788	0.00000
Cranio	9.90532	0.42327	23.40217	0.00000
SessoM	78.24227	10.99599	7.11552	0.00000

I coefficienti ci indicano come variabili signficative il tempo di Gestazione, la Lunghezza del neonato, la dimensione del Cranio e il Sesso del neonato. Dalla matrice di correlazione possiamo vedere anche graficamente la correlazione tra Peso e Gestazione, Lunghezza e Cranio. Le curve di fit create automaticamente nella matrice sembrano indicare dei comportamenti non lineari che potrebbero essere però dovuti ad over fitting e che potremo verificare nella procedura stepwise.

L’R quadro aggiustato è di 0.7366 che ci indica un buon modello di descrizione dei dati.

Stepwise: vediamo di migliorare il modello attraverso la procedura stepwise.

Valutiamo il modello escludendo le variabili che sono meno significative:

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1165.68  -179.74   -12.42   162.92  1410.68 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6683.8326   133.1602 -50.194  < 2e-16 ***
## N.gravidanze    13.1448     4.2576   3.087  0.00204 ** 
## Gestazione      29.6341     3.7369   7.930 3.27e-15 ***
## Lunghezza       10.8899     0.3019  36.077  < 2e-16 ***
## Cranio           9.9192     0.4227  23.465  < 2e-16 ***
## SessoM          78.1376    10.9929   7.108 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7372, Adjusted R-squared:  0.7367 
## F-statistic:  1398 on 5 and 2491 DF,  p-value: < 2.2e-16

Togliendo le due variabili l’R quadro aggiustato è cresciuto a 0.7367

Valutiamo l’effetto della variabile N.gravidanze:

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1153.75  -182.01   -14.88   163.78  1391.18 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6652.115    132.990 -50.020  < 2e-16 ***
## Gestazione     28.533      3.726   7.657 2.69e-14 ***
## Lunghezza      10.841      0.302  35.903  < 2e-16 ***
## Cranio         10.059      0.421  23.893  < 2e-16 ***
## SessoM         79.321     11.005   7.208 7.51e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.8 on 2492 degrees of freedom
## Multiple R-squared:  0.7362, Adjusted R-squared:  0.7358 
## F-statistic:  1739 on 4 and 2492 DF,  p-value: < 2.2e-16

L’R quadro aggiustato è calato da 0.7367 a 0.7358 ma ci lascia un modello con quattro variabili (Gestazione, Lunghezza, Cranio e Sesso) tutte molto significative. mod3 potrebbe essere un buon modello.

Provo a verificare l’effetto di eventuali termini quadratici che avevamo osservato nella matrice di correlazione.

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Lunghezza^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.66  -182.74   -12.38   164.68  1307.78 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.666e+03  7.603e+02  -2.192 0.028481 *  
## Gestazione      3.644e+01  3.880e+00   9.391  < 2e-16 ***
## Lunghezza      -1.137e+01  3.348e+00  -3.395 0.000698 ***
## Cranio          1.030e+01  4.189e-01  24.577  < 2e-16 ***
## SessoM          7.358e+01  1.094e+01   6.723 2.19e-11 ***
## I(Lunghezza^2)  2.288e-02  3.436e-03   6.659 3.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7408, Adjusted R-squared:  0.7403 
## F-statistic:  1424 on 5 and 2491 DF,  p-value: < 2.2e-16

Con il termine quadratico relativo alla lunghezza l’R quadro aggiustato cresce a 0.7403 e non si vede perdita di sigificatività di altre variabili.

Per migliorare la descrizione del mod4 ed evitare di avere problemi di collinearità elevata tra Lunghezza e Lunghezza^2 utilizziamo il pacchetto poly per ridefinire mod4:

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.66  -182.74   -12.38   164.68  1307.78 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1673.8583   197.6607  -8.468  < 2e-16 ***
## Gestazione             36.4395     3.8803   9.391  < 2e-16 ***
## Cranio                 10.2953     0.4189  24.577  < 2e-16 ***
## SessoM                 73.5847    10.9444   6.723 2.19e-11 ***
## poly(Lunghezza, 2)1 13576.7873   398.9042  34.035  < 2e-16 ***
## poly(Lunghezza, 2)2  1888.3743   283.5889   6.659 3.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7408, Adjusted R-squared:  0.7403 
## F-statistic:  1424 on 5 and 2491 DF,  p-value: < 2.2e-16

Verifico anche il termine quadratico con Gestazione

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2) + I(Gestazione^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1174.11  -184.47   -12.53   163.22  1307.85 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -4654.2522  1265.6807  -3.677 0.000241 ***
## Gestazione            194.3051    66.3329   2.929 0.003429 ** 
## Cranio                 10.2454     0.4190  24.450  < 2e-16 ***
## SessoM                 74.7566    10.9452   6.830 1.06e-11 ***
## poly(Lunghezza, 2)1 13202.2085   428.3852  30.819  < 2e-16 ***
## poly(Lunghezza, 2)2  2498.7550   381.8715   6.543 7.27e-11 ***
## I(Gestazione^2)        -2.0729     0.8695  -2.384 0.017201 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.2 on 2490 degrees of freedom
## Multiple R-squared:  0.7414, Adjusted R-squared:  0.7408 
## F-statistic:  1190 on 6 and 2490 DF,  p-value: < 2.2e-16

L’aggiunta di un termine quadratico per gestazione mostra invece una sigificatività sotto al limite ma non incisiva e porta invece a far calare il coefficiente di sigifcatività di Gestazione. Scartiamo il mod5.

Verifichiamo eventuali interazioni tra variabili

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2) + Lunghezza + Gestazione:Lunghezza)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1175.45  -184.34   -11.54   166.47  1302.13 
## 
## Coefficients: (1 not defined because of singularities)
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1598.0934   206.6423  -7.734 1.51e-14 ***
## Gestazione             103.1573    53.2721   1.936 0.052929 .  
## Cranio                  10.2658     0.4195  24.471  < 2e-16 ***
## SessoM                  74.4364    10.9642   6.789 1.41e-11 ***
## poly(Lunghezza, 2)1  20192.6018  5283.5724   3.822 0.000136 ***
## poly(Lunghezza, 2)2   2554.4653   601.4744   4.247 2.25e-05 ***
## Lunghezza                    NA         NA      NA       NA    
## Gestazione:Lunghezza    -0.1381     0.1099  -1.256 0.209331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.5 on 2490 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:  0.7403 
## F-statistic:  1187 on 6 and 2490 DF,  p-value: < 2.2e-16

L’interazione Gestazione/Lunghezza non produce miglioramenti nel modello.

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2) + Lunghezza + Cranio:Lunghezza)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.4  -182.1   -13.3   165.6  1368.8 
## 
## Coefficients: (1 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.729e+03  1.996e+02  -8.661  < 2e-16 ***
## Gestazione           3.680e+01  3.883e+00   9.479  < 2e-16 ***
## Cranio              -2.118e+00  6.475e+00  -0.327   0.7436    
## SessoM               7.276e+01  1.095e+01   6.647 3.67e-11 ***
## poly(Lunghezza, 2)1  2.896e+03  5.574e+03   0.520   0.6034    
## poly(Lunghezza, 2)2  9.257e+02  5.757e+02   1.608   0.1079    
## Lunghezza                   NA         NA      NA       NA    
## Cranio:Lunghezza     2.530e-02  1.317e-02   1.921   0.0548 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.3 on 2490 degrees of freedom
## Multiple R-squared:  0.7412, Adjusted R-squared:  0.7406 
## F-statistic:  1188 on 6 and 2490 DF,  p-value: < 2.2e-16

L’interazione Cranio/Lunghezza non produce miglioramenti nel modello.

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2) + Lunghezza + Sesso:Lunghezza)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1160.65  -182.62   -12.86   166.26  1294.38 
## 
## Coefficients: (1 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1677.7177   197.8126  -8.481  < 2e-16 ***
## Gestazione             36.4897     3.8819   9.400  < 2e-16 ***
## Cranio                 10.3025     0.4192  24.579  < 2e-16 ***
## SessoM                190.9612   213.5338   0.894    0.371    
## poly(Lunghezza, 2)1 13706.0909   462.9878  29.604  < 2e-16 ***
## poly(Lunghezza, 2)2  1921.2943   289.8662   6.628 4.15e-11 ***
## Lunghezza                   NA         NA      NA       NA    
## SessoM:Lunghezza       -0.2371     0.4308  -0.550    0.582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.5 on 2490 degrees of freedom
## Multiple R-squared:  0.7408, Adjusted R-squared:  0.7402 
## F-statistic:  1186 on 6 and 2490 DF,  p-value: < 2.2e-16

Anche l’interazione Sesso/Lunghezza non produce miglioramenti nel modello.

Criteri di informazione di Akaike (AIC) e Bayes (BIC) per confronto dei modelli: utilizziamo i criteri AIC e BIC di confronto tra i modelli per vedere quale sia il più efficiente, e il modello VIF

Criterio di informazione AIC
	df	AIC
mod1	9	35043.32
mod2	7	35040.64
mod3	6	35048.18
mod4	7	35006.12
mod5	8	35002.43
mod6	8	35007.82

Criterio di informazione BIC
	df	BIC
mod1	9	35095.73
mod2	7	35081.40
mod3	6	35083.12
mod4	7	35046.88
mod5	8	35049.01
mod6	8	35054.40

Verifichiamo anche problemi di collinearità con la funzione vif

Verifica collinearità
	GVIF	Df	GVIF^(1/(2*Df))
Gestazione	1.835138	1	1.354672
Cranio	1.650092	1	1.284559
Sesso	1.044984	1	1.022245
poly(Lunghezza, 2)	2.391268	2	1.243532

Il mod4 sembra essere il più efficace in entrambi i criteri di verifica e i valori del test vif non mostrano problemi di collinearità avendo tutti valori inferiori a 5.

Facciamo una verifica anche con la funzione in r che calcola stepwise in automatico con il pacchetto MASS.

## Start:  AIC=28001.73
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Fumatrici + Lunghezza + 
##     Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     20903 180645287 27994
## - Fumatrici     1     73977 180698361 27995
## - N.gravidanze  1    534941 181159325 28001
## <none>                      180624384 28002
## - Sesso         1   3674222 184298606 28044
## - Gestazione    1   4632038 185256422 28057
## - Cranio        1  39743276 220367660 28491
## - Lunghezza     1  93986157 274610541 29040
## 
## Step:  AIC=27994.19
## Peso ~ N.gravidanze + Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     74619 180719906 27987
## <none>                      180645287 27994
## - N.gravidanze  1    712872 181358159 27996
## + Anni.madre    1     20903 180624384 28002
## - Sesso         1   3679797 184325084 28037
## - Gestazione    1   4617826 185263113 28049
## - Cranio        1  39918188 220563475 28485
## - Lunghezza     1  93992253 274637540 29032
## 
## Step:  AIC=27987.4
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      180719906 27987
## - N.gravidanze  1    691520 181411426 27989
## + Fumatrici     1     74619 180645287 27994
## + Anni.madre    1     21545 180698361 27995
## - Sesso         1   3665488 184385394 28030
## - Gestazione    1   4562528 185282434 28042
## - Cranio        1  39945921 220665827 28478
## - Lunghezza     1  94424525 275144431 29029

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1165.68  -179.74   -12.42   162.92  1410.68 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6683.8326   133.1602 -50.194  < 2e-16 ***
## N.gravidanze    13.1448     4.2576   3.087  0.00204 ** 
## Gestazione      29.6341     3.7369   7.930 3.27e-15 ***
## Lunghezza       10.8899     0.3019  36.077  < 2e-16 ***
## Cranio           9.9192     0.4227  23.465  < 2e-16 ***
## SessoM          78.1376    10.9929   7.108 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7372, Adjusted R-squared:  0.7367 
## F-statistic:  1398 on 5 and 2491 DF,  p-value: < 2.2e-16

La funzione automatica di r ci restituisce un modello differente da mod4. Rieffettuiamo i test AIC e BIC sui modelli mod4 e stepwise.mod per decidere quale selezionare:

Criterio di informazione AIC
	df	AIC
mod4	7	35006.12
stepwise.MOD	7	35040.64

Criterio di informazione BIC
	df	BIC
mod4	7	35046.88
stepwise.MOD	7	35081.40

Selezione del modello Visto che in entrambi i test il valore più basso è legato al modello costruito manualmente con la procedura stepwise, confermiamo la selezione del modello 4 come migliore per descrivere il comportamento della variabile peso in relazione alle altre del nostro campione. mod4 ha come variabili correlate Gestazione, Cranio, Sesso e un polinomio di secondo grado della variabile Lunghezza.

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Sesso + poly(Lunghezza, 
##     2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.66  -182.74   -12.38   164.68  1307.78 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1673.8583   197.6607  -8.468  < 2e-16 ***
## Gestazione             36.4395     3.8803   9.391  < 2e-16 ***
## Cranio                 10.2953     0.4189  24.577  < 2e-16 ***
## SessoM                 73.5847    10.9444   6.723 2.19e-11 ***
## poly(Lunghezza, 2)1 13576.7873   398.9042  34.035  < 2e-16 ***
## poly(Lunghezza, 2)2  1888.3743   283.5889   6.659 3.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7408, Adjusted R-squared:  0.7403 
## F-statistic:  1424 on 5 and 2491 DF,  p-value: < 2.2e-16

Analisi della qualità del modello

Per valutare la bontà del modello possiamo controllare i valori di R^2, R^2 aggiustato e del Root Mean Square Error. Calcoliamo RMSE

data.frame <- dati
prediction <- predict(mod4, newdata = dati)
RMSE = sqrt(mean((dati$Peso-prediction)^2))
round(RMSE,4)

## [1] 267.1725

R^2 = 0.7408
R^2 aggiustato = 0.7403
RMSE = 267.1725 g

I valori di R^2 indicano che il modello descrive con una bontà del 74.03% l’andamento dei dati, mentre l’RMSE ci dice che l’errore medio tra predizione del modello e dati misurati sulla variabile peso è di circa 267 g (su un valore medio di Peso di 3283 g).

Parte erratica del modello Ricordiamo che dai dati iniziali abbiamo già eliminato 3 misurazioni: due per età errata della madre (0 e 2 anni), e una per una valutazione anomala del rapporto peso lunghezza

Proseguo valutando i residui del modello:

Il grafico in alto a sinistra mostra la dispersione attorno alla media, che è correttamente 0. La nuvola di punti non presenta pattern particolari, per cui non sembra esser stata persa informazione dal modello. La linea rossa, come ci aspettavamo, non mostra particolari flessioni.
Il grafico in alto a destra mette in relazione i residui con una distribuzione normale. I punti si adattano bene (giacciono sulla retta bisettrice) a parte per i valori più esterni della curva che mostrano delle flessioni rispetto alla bisettrice.
Nel grafico in basso a sinistra i punti si distribuiscono in maniera abbastanza casuale attorno ad un valore medio di y che sembra essere costante. Si osserva forse una leggera tendenza positiva della linea ma non sostanziale.
Nel grafico in basso a destra si valutano eventuali misurazioni outliers o leverage, ma nel nostro grafico non risultano esserci problemi di valori influenti.

Verifichiamo queste osservazioni anche in maniera più dettagliata

Valori di leva

lev <-hatvalues(mod4)
plot(lev, pch = 20, col = "LightBlue")
p = sum(lev)
soglia = 2*p/nrow(dati)

abline(h = soglia, col = "Red")

length(lev[lev>soglia])

## [1] 137

Sono state osservate 137 osservazioni sopra il valore di soglia.

Valori outliers

plot(rstudent(mod4), pch = 20, col = "LightBlue")
abline(h = c(-2,2), col = "Red")

outlierTest(mod4)

##       rstudent unadjusted p-value Bonferroni p
## 155   4.934357         8.5724e-07    0.0021405
## 1305  4.864857         1.2169e-06    0.0030385
## 1397 -4.370966         1.2877e-05    0.0321550
## 1691  4.331301         1.5410e-05    0.0384790

Nonostante il numero di outliers sia superiore, con il test t e la correzione di Bonferroni sono stati evidenziati solamente 4 osservazioni molto distanti dai valori di soglia.

Possiamo considerare Outliers e Leverage in contemporanea valutando la distanza di Cook:

cook <- cooks.distance(mod4)
plot(cook, pch = 20, col = "LightBlue")

max(cook)

## [1] 0.06576053

Il massimo della distanza di Cook è di 0.06, molto al di sotto della soglia di allarme di 0.5, per cui valutiamo che le misure estreme non siano realmente così influenti nella valutazione dei residui.

In ultimo facciamo i test di Breusch-Pagan e di Durbin-Watson per la omoschedasticità e la autocorrelazione.

bptest(mod4)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod4
## BP = 13.578, df = 5, p-value = 0.01853

dwtest(mod4)

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

shapiro.test(residuals(mod4))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod4)
## W = 0.98968, p-value = 1.982e-12

Risultati dei test:

il p value di BP molto basso ci mostra che c’è eteroschedasticità
la stastica DW vicino a 2 e il p value ci dicono che i residui non sono autocorrelati
Il test SW per la normalità ha un p value molto basso che ci porta a dire che i dati non sono assimilabili ad una distribuzione normale.

Conclusioni

Il modello 4 risulta in definitiva un modello abbastanza affidabile anche se con alcune imprecisioni le cui cause si possono studiare più approfonditamente. Va osservato che nel modello le variabili più rilevanti nella previsione del peso neonatale sono quelle antropometiche (Lunghezza, dimensione del cranio) più il sesso del bambino;solamente le settimane di gestazione sono misurabili senza una osservazione ecografica del feto. In più evidenziamo come le variabili legate alla madre (Fumatrice e gravidanze precedenti) risultino poco influenti, e quindi escluse dal modello, nella previsione del peso. Come ci si poteva aspettare in ultimo le variabili legate a questioni esterne alla crescita del feto come la scelta dell’ospedale o le metodologie mediche per il parto non hanno influenza sul peso.

3. Previsioni e Risultati

Utilizziamo il modello selezionato per effettuare delle previsioni di peso.

Previsione 1

Valutiamo la previsione del peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

Estraiamo i valori dai dati complessivi che verificano queste caratteristiche e ne tabuliamo alcuni valori di posizione come confronto con quello che il modello predirrà:

filtro <- dati[N.gravidanze == 3 & Gestazione == 39 & Sesso == "F", ]
filtro_tab <- filtro %>%
    summarize(
      peso_medio = round(mean(Peso),2),
      peso_mediano = round(median(Peso),2),
      peso_min = round(min(Peso),2),
      peso_max = round(max(Peso),2)
    )

Valori di partenza
peso_medio	peso_mediano	peso_min	peso_max
3248.82	3250	2680	3880

Previsione modello

test1 <- data.frame(N.gravidanze = 3, Gestazione = 39, Lunghezza = mean(Lunghezza), Cranio = mean(Cranio), Sesso = "F")

predict(mod4, newdata = test1)

##        1 
## 3232.285

Il modello ci da un previsione di 3232 g circa per il peso di una neonata Femmina alla 39esima settimana, terza gravidanza della madre.

Previsione 2

Facciamo una previsione partendo da una misurazione che contenga tutte le variabili del nostro modello:

gestazione = 27
cranio = 348
sesso = maschio
lunghezza = 510

Non ci sono dati nel dataset iniziale che verifichino queste quattro condizioni.

Previsione modello

test1 <- data.frame(Gestazione = 27, Lunghezza = 510, Cranio = 348, Sesso = "M")

predict(mod4, newdata = test1)

##        1 
## 3127.825

La previsione del modello è di 3127g per un neonato maschio alla 27esima settimana, con dimensioni antropometriche 510mm lunghezza e 348mm dimensione del cranio.

4. Visualizzazioni

Utilizziamo rappresentazioni grafiche per mettere in evidenza il comportamento del modello e alcune dipendenze sostanziali tra variabili.

I due scatterplot 3D sono utili per vedere graficamente la dipendenza del peso dalle variabili che abbiamo evidenziato come più influenti nel modello sviluppato.

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

Il grafico mette in evidenza l’aumento di peso neonatale al crescere delle settimane di gestazione. L’aspetto da evidenziare è la differente inclinazione legata alla variabile Fumatrice, che mostra un aumento di peso meno marcato per il campione di madri Fumatrici. Il dato può essere interessante, anche se nel nostro modello abbiamo evidenzziato come sia marginale. Inoltre il campione di madri fumatrici rappresenta il 4.16% del totale delle misurazioni, quindi un numero piccolo che può essere insufficiente per evidenziare in maniera esaustiva un differente andamento tra Fumatrici e non Fumatrici. Infatti nella tabulazione del Modello1 (modello iniziale che conteneva tutte le variabili) il p-value associato alla variabili Fumatrici era di 0.31, confermando il fatto che l’effetto che si può vedere non può essere considerato statisticamente significativo.

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

Il grafico moette in evidenza la differenza di peso tra maschi e femmine al crescere del tempo di gestazione. E’ subito evidente come le due linee di plot abbiano la stessa pendenza, ma intercette e valori più alti per i maschi. La differenza di andamento è quella che ci ha portato ad inserire la variabile Sesso come siginificativa nel modello.

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

Nel modello abbiamo incluso la varibile Lunghezza come polinomio quadratico. Il grafico mostra proprio l’andamento del Peso in funzione dalla Lunghezza al quadrato e la retta di plot in rosso conferma graficamente il legame tra queste due grandezze.

I due grafici precedenti mostrano l’andamento del Peso in funzione dell’Età della madre e del Numero di gravidanze. Sono due variabili che non abbiamo incluso nel modello come variabili significative nella valutazione del Peso; i due scatter plot non mostrano dipendenze evidenti e quindi confermano la valutazione del modello.

L’ultimo grafico mette a confronto le previsioni di peso estratte dal modello con quelle misurate nel database. Si può vedere il buon accordo del modello, soprattutto per le misurazioni di peso attorno al valore medio. Si può vedere forse una leggera sottostima del modello rispetto ai dati misurati per valori di peso molto alti.

Conclusioni

Il modello costruito si presta ad essere un buon modello per descrivere il comportamento del peso neonatale soprattutto per i valori centrali della crescita del neonato. Durante l’analisi preliminare abbiamo verificato come le tecniche ospedaliere dei tre siti in analisi siano già abbastanza omologate, aspetto fondamentale per poter studiare bene i dati e poter implementare metodologie comuni di lavoro. Sempre all’interno dell’analisi preliminare abbiamo valutato la differenza delle misure antropometriche del neonato tra i due sessi. Proprio le variabili antropometriche sono quelle che abbiamo determinato essere significative per lo sviluppo del modello. Nel modello infatti abbiamo selezionato Lunghezza, dimensione del Cranio, Sesso e settimane di Gestazione. Attreverso queste variabili il modello ha una variabilità spiegata del 74%. Lo sviluppo del modello aveva come scopo il miglioramento delle delle previsioni cliniche ospedaliere in modo da poter evidenziare anomalie e fattori di rischio, e proprio grazie ad una capacità di predizione buona per neonati in salute il modello può aiutarci (abbiamo eliminato nella verifica iniziale dei dati tre valori che mostravano misurazioni anomale, probabilemente date da errori di inserimento o anomalie nei neonati). Per valutare nel dettaglio l’effetto e la provenienza di eventuali anomalie rispetto al modello atteso servirebbe fare una anlisi specifica delle variazioni introdotte nel modello dalle variabili escluse, che pur non avendo effetti significativi nella stima del peso possono essere utili per determinare eventuali variazioni non sostanziali al peso e quindi alla salute del neonato.

Modello Statistico per la Previsione del Peso Neonatale

Giusti Nicola

2025-11-23

1) Raccolta dat e struttura del dataset

2) Analisi e modellizzazione

Analisi preliminare

Test sui parti cesarei per ospedale

Analisi peso e lunghezza del campione

Analisi misure antropometriche

Creazione modello di regressione

Analisi della qualità del modello

Conclusioni

3. Previsioni e Risultati

4. Visualizzazioni

Conclusioni