---

Modello statistico per la previsione del peso dei neonati

Su richiesta dell’azienda Neonatal Health Solutions è stato sviluppato questo documento che riporta il processo di analisi e costruzione di un modello statistico in grado di prevedere il peso dei neonati.

Dopo un’analisi descrittiva dei dati forniti con particolare enfasi sulla durata della gestazione e sul fumo materno, si effettueranno test statistici sui dati. Infine, dopo aver costruito il modello e dopo averne valutata la qualità, si svolgeranno predizioni su dati nuovi e si analizzeranno i risultati ottenuti.

Analisi descrittiva del modello

Questa sezione sarà dedicata all’importazione e a un’analisi descrittiva dei dati del campione.

I dati sono stati forniti dall’azienda committente tramite un file CSV che è stato importato nel Dataframe R df.

df <- read.csv("neonati.csv", sep = ",")

Visualizziamone le dimensioni e memorizziamo la grandezza del campione.

dim(df)

## [1] 2500   10

n <- dim(df)[1]

Il campione riguarda le misurazioni di 10 variabili effettuate su 2500 bambini. La cella seguente mostra le prime cinque righe del campione.

head(df,5)

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

Le variabili saranno esplorate in tre sottosezioni. Nella prima ci focalizzeremo sulla variabile target del futuro modello, nella seconda valuteremo le variabili numeriche ed effettueremo una prima selezione per il modello, mentre nella terza ci occuperemo delle variabili di controllo che possono influire sulle variabili numeriche.

Il target Peso

In questa sottosezione ci occuperemo della variabile Peso, ossia il target del modello.

I valori di questa variabile, espressi in grammi, rappresentano una variabile quantitativa di tipo continuo. Di seguito si mostrano il minimo e il massimo valore di questa colonna.

range(df$Peso)

## [1]  830 4930

La cella seguente divide la colonna Peso in classi di 250 grammi e mostra la distribuzione di queste classi in un grafico a barre. Per facilitare la lettura del grafico, i valori sono stati convertiti in ettogrammi.

Peso100 <- df$Peso/100

min_limit <- floor(min(Peso100))
max_limit <- ceiling(max(Peso100))
Peso_CL <- cut(Peso100, breaks = seq(min_limit, max_limit, 2.50))

barplot(table(Peso_CL), las=2)

Nonostante la forma a campana, la distribuzione presenta uno sbilanciamento a favore dei pesi più elevati. La non normalità della distribuzione è verificabile tramite il test di Shapiro-Wilk nella cella che segue.

shapiro.test(df$Peso)

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

Il valore prossimo allo zero del p-value ci spinge a rifiutare l’ipotesi nulla e a confermare che i dati di peso non sono distribuiti normalmente.

Per visualizzare la relazione tra Peso e le altre variabili quantitative, è stata costruita la matrice di correlazione ed è stata estratta la colonna relativa al peso. La cella seguente ordina questi valori in modo decrescente.

mat.cor <- cor(df[c("Anni.madre",
         "N.gravidanze",
         "Gestazione",
         "Peso",
         "Lunghezza",
         "Cranio")])

as.data.frame(sort(mat.cor[,"Peso"], decreasing = T))

##              sort(mat.cor[, "Peso"], decreasing = T)
## Peso                                      1.00000000
## Lunghezza                                 0.79603676
## Cranio                                    0.70480151
## Gestazione                                0.59176872
## N.gravidanze                              0.00240730
## Anni.madre                               -0.02247017

Queste relazioni, insieme alle singole variabili quantitative, saranno esplorate nella prossima sottosezione.

Le variabili quantitative

Questa sottosezione è dedicata alle variabili quantitative del dataset e alla loro relazione con la variabile target. Nella cella seguente si mostrano le prime cinque righe di questi valori.

head(df[c("Peso","Lunghezza","Cranio","Gestazione","N.gravidanze","Anni.madre")],5)

##   Peso Lunghezza Cranio Gestazione N.gravidanze Anni.madre
## 1 3380       490    325         42            0         26
## 2 3150       490    345         39            2         21
## 3 3640       500    375         38            3         34
## 4 3690       515    365         41            1         28
## 5 3700       480    335         38            0         20

Di ogni variabile saranno visualizzati il minimo e il massimo valore, la loro distribuzione e la loro relazione con la variabile target.

Gestazione

Seguendo le indicazioni fornite dall’azienda committente, inizieremo la nostra analisi dalla colonna Gestazione, una variabile quantitativa discreta che rappresenta la durata in settimane della gravidanza. Di seguito se ne mostra il range di valori.

range(df$Gestazione)

## [1] 25 43

I valori mostrano un range tra le 25 settimane (intorno al secondo trimestre) e le 43 settimane (oltre il terzo trimestre), il che giustifica il range di peso visto nella sottosezione precedente.

La cella seguente mostra in un grafico a barre la distribuzione di questi dati.

barplot(table(df$Gestazione),las=2)

Il grafico a dispersione che segue mostra invece la relazione tra Gestazione e Peso.

plot(df$Gestazione, df$Peso, pch=20)

La maggioranza dei dati si concentra dopo la 35° settimana di gravidanza e presenta una relazione positiva con il peso. Inoltre, dalla forma del grafico, si potrebbe ipotizzare una dipendenza non lineare tra le due variabili probabilmente influenzata da altri fattori come il fumo materno.

Lunghezza

La colonna Lunghezza è una variabile quantitativa continua che misura la lunghezza dei neonati in millimetri. La cella seguente mostra il range di valori di questa colonna, compresi tra i 31 e i 57 centimetri.

range(df$Lunghezza)

## [1] 310 565

La cella seguente mostra la distribuzione di Lunghezza in classi da 1 centimetro.

min_limit <- round(min(df$Lunghezza), digits = -2)
max_limit <- round(max(df$Lunghezza), digits = -2)
Lunghezza_CL <- cut(df$Lunghezza, breaks = seq(min_limit, max_limit, 10))

barplot(table(Lunghezza_CL), las=2)

Il grafico a dispersione mostra invece la colonna Lunghezza in relazione a Peso.

plot(df$Lunghezza, df$Peso, pch=20)

Rispetto a Gestazione, con cui condivide la correlazione positiva con Peso, i dati variano in una fascia più stretta intorno all’ipotetica retta che mette in relazione le due colonne.

Cranio

La colonna Cranio esprime in millimetri la circonferenza del cranio del neonato. La cella che segue mostra il range di valori, compresi tra i 23 e i 39 cm.

range(df$Cranio)

## [1] 235 390

La cella che segue mostra la distribuzione di questa colonna nel dataset in classi di un centimetro.

min_limit <- round(min(df$Cranio), digits = -2)
max_limit <- round(max(df$Cranio), digits = -2)
Cranio_CL <- cut(df$Cranio, breaks = seq(min_limit, max_limit, 10))

barplot(table(Cranio_CL), las=2)

Infine, si mostra la relazione tra Cranio e Peso con un grafico a dispersione.

plot(df$Cranio, df$Peso, pch=20)

I risultati ottenuti con la colonna Cranio assomigliano molto a quelli di Lunghezza, il che fa sospettare che le due variabili possano essere correlate tra loro. Per verificarlo, abbiamo confrontato le correlazioni tra queste colonne con Peso prima singolarmente e poi come interazione.

paste("Peso~Lunghezza",cor(df$Lunghezza, df$Peso))

## [1] "Peso~Lunghezza 0.796036760877904"

paste("Peso~Cranio",cor(df$Cranio, df$Peso))

## [1] "Peso~Cranio 0.704801506528901"

paste("Peso~Lunghezza*Cranio",cor(df$Lunghezza*df$Cranio, df$Peso))

## [1] "Peso~Lunghezza*Cranio 0.840261911395334"

Unendo le colonne Lunghezza e Cranio otteniamo una maggiore correlazione con il peso rispetto alle stesse colonne singolarmente, il che ci porta a considerare nel modello la loro interazione.

Numero di gravidanze pregresse

La colonna N.gravidanze è una variabile discreta che conta il numero di gravidanze precedenti a quella osservata nel dataset. Di seguito si mostra il range di valori.

range(df$N.gravidanze)

## [1]  0 12

Dal momento che non abbiamo ulteriori informazioni sulle madri, possiamo solo dedurre che il range di valori, benché insolito, sia comunque plausibile e spiegabile da fattori che esulano dallo scopo del progetto (ad esempio, gravidanze non portate a termine, famiglie numerose o anche parti gemellari). Vediamone comunque la distribuzione nella cella che segue.

barplot(table(df$N.gravidanze))

Osserviamo che il numero di neonati diminuisce all’aumentare delle gravidanze e che la maggioranza dei dati riguarda madri che hanno avuto al massimo due gravidanze precedenti.

La cella che segue esamina N.gravidanze in relazione con Peso.

plot(df$N.gravidanze, df$Peso, pch=20)

Dal grafico emerge che, all’aumentare delle gravidanze, i range di peso si restringono verso valori più elevati. Ciononostante, come anche evidenziato nella matrice di correlazione, le variabili N.gravidanze e Peso presentano una scarsissima correlazione tra di loro, il che ci porta a non dare priorità a questa variabile nella costruzione del modello.

Età materna

La colonna Anni.madre indica l’età materna al momento del parto ed è una variabile discreta. Di seguito mostriamo il range di valori.

range(df$Anni.madre)

## [1]  0 46

Dal momento che è poco plausibile una madre neonata, la conclusione più ovvia è che il minimo valore della colonna sia in realtà un errore di inserimento dei dati. Per escludere questo tipo di errore anche per altri record, mostriamo di seguito le prime cinque righe del dataset ordinate in modo crescente in base all’età materna.

head(sort_by(df,df$Anni.madre),5)

##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1380          0            0         0         39 3060       490    330
## 1152          1            1         0         41 3250       490    350
## 138          13            0         0         38 2760       470    325
## 1075         14            1         0         39 3510       490    365
## 1532         14            0         0         39 3550       500    355
##      Tipo.parto Ospedale Sesso
## 1380        Nat     osp3     M
## 1152        Nat     osp2     F
## 138         Nat     osp2     F
## 1075        Nat     osp2     M
## 1532        Ces     osp1     M

Da questi risultati confermiamo quindi che si è verificato un errore di inserimento per cui i due valori più piccoli di Anni.madre fossero 0 e 1, mentre è plausibile la presenza di madri minori di quindici anni.

Onde evitare che questo errore possa inficiare nel modello, si è deciso di sostituire i due valori errati con la media dei valori corretti arrotondata all’intero e di ricalcolare la matrice di correlazione.

df[df$Anni.madre %in% c(0,1),"Anni.madre"]=round(sum(df$Anni.madre-1)/(n-2))

mat.cor <- cor(df[c("Anni.madre",
         "N.gravidanze",
         "Gestazione",
         "Peso",
         "Lunghezza",
         "Cranio")])

La cella che segue mostra le prime cinque righe del dataset corretto ordinato per età materna.

head(sort_by(df,df$Anni.madre),5)

##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 138          13            0         0         38 2760       470    325
## 1075         14            1         0         39 3510       490    365
## 1532         14            0         0         39 3550       500    355
## 66           15            0         0         37 2810       490    325
## 415          15            0         0         39 3440       510    347
##      Tipo.parto Ospedale Sesso
## 138         Nat     osp2     F
## 1075        Nat     osp2     M
## 1532        Ces     osp1     M
## 66          Ces     osp3     F
## 415         Ces     osp1     F

La cella che segue mostra nuovamente le correlazioni con il peso.

as.data.frame(sort(mat.cor[,"Peso"], decreasing = T))

##              sort(mat.cor[, "Peso"], decreasing = T)
## Peso                                      1.00000000
## Lunghezza                                 0.79603676
## Cranio                                    0.70480151
## Gestazione                                0.59176872
## N.gravidanze                              0.00240730
## Anni.madre                               -0.02373525

Ora che il dataset è corretto, possiamo osservare come l’età materna sia più correlata con il peso, anche se in modo negativo. Vediamo dunque la distribuzione per fasce di età di cinque anni, dove emerge un’asimmetria dei dati verso età più giovani delle madri.

Anni_CL <- cut(df$Anni.madre, breaks = seq(10, 50, 5))

barplot(table(Anni_CL), las=2)

Infine, visualizziamo l’età materna in relazione al peso.

plot(df$Anni.madre, df$Peso, pch=20)

Il grafico a dispersione mostra una scarsissima correlazione tra l’età materna e il peso. La correlazione negativa potrebbe essere dovuta al fatto che molti dei neonati sotto i 2000 grammi si trovano nella fascia di età tra i 25 e i 40 anni.

Conclusioni sulle variabili quantitative

Per fare un primo punto dell’analisi descrittiva, dai risultati ottenuti possiamo dedurre quanto segue sul modello:

• il peso è più correlato alle misure antropometriche rispetto alla durata della gestazione;

• le misure di lunghezza e cranio non incidono sul peso solo singolarmente, ma anche in interazione tra di loro;

• la correlazione negativa tra età materna e peso è spiegabile dalla forte presenza di neonati sotto i 2000 grammi nella fascia di età tra i 25 e i 40 anni.

In virtù di queste considerazioni, quindi, possiamo già concludere che nel modello saranno prese in considerazione le colonne Gestazione, Lunghezza e Cranio, mentre le colonne N.gravidanze e Anni.madre avranno una minore priorità.

Le variabili di controllo

Questa sottosezione è dedicata alle variabili di controllo, ossia le variabili categoriche all’interno del dataset. La cella che segue mostra le prime cinque righe di queste variabili insieme alla colonna Peso.

head(df[c("Peso","Fumatrici","Tipo.parto","Ospedale","Sesso")],5)

##   Peso Fumatrici Tipo.parto Ospedale Sesso
## 1 3380         0        Nat     osp3     M
## 2 3150         0        Nat     osp1     F
## 3 3640         0        Nat     osp2     M
## 4 3690         0        Nat     osp2     M
## 5 3700         0        Nat     osp3     F

Per ogni variabile mostreremo la distribuzione e la sua influenza sulla relazione tra il peso e le colonne Gestazione, Lunghezza e Cranio sia dal punto di vista grafico che dal punto di vista delle correlazioni.

Fumo materno

Sempre su richiesta dell’azienda committente, inizieremo la nostra analisi dalla colonna Fumatrici, una variabile binaria in cui è indicato se la madre fuma (1) o meno (0). La cella seguente mostra la distribuzione dei valori di questa colonna nel dataset.

barplot(table(df$Fumatrici))

Si osserva una disparità tra fumatrici (in minoranza) e non fumatrici. Non avendo ulteriori informazioni sulle madri, questo sbilanciamento può essere spiegato da molti fattori (ad esempio, madri che hanno mentito durante l’anamnesi o che hanno smesso da poco di fumare) che tuttavia esulano dallo scopo di questo progetto.

Le tre celle seguenti filtrano i grafici a dispersione visti in precedenza sulla base del fumo materno.

library(ggplot2)

ggplot(df)+
  geom_point(aes(x=Gestazione, y=Peso))+
  facet_grid(~Fumatrici)

ggplot(df)+
  geom_point(aes(x=Lunghezza, y=Peso))+
  facet_grid(~Fumatrici)

ggplot(df)+
  geom_point(aes(x=Cranio, y=Peso))+
  facet_grid(~Fumatrici)

I tre grafici mostrano che, a parità di gestazione e misure, i neonati delle fumatrici presentano una fascia di peso meno estesa e più bassa rispetto alle non fumatrici.

Per quantificare l’effetto del fumo materno, dividiamo il dataset in base alla colonna Fumatrici e calcoliamo le correlazioni tra queste variabili. I risultati sono mostrati nella cella che segue

df_splitted <- split(df, df$Fumatrici)

data.frame(
  Colonna = c("Gestazione","Lunghezza","Peso"),
  Correlazione0 = c(cor(df_splitted$`0`$Gestazione,df_splitted$`0`$Peso),
                     cor(df_splitted$`0`$Lunghezza,df_splitted$`0`$Peso),
                     cor(df_splitted$`0`$Cranio,df_splitted$`0`$Peso)),
  
  Correlazione1 = c(cor(df_splitted$`1`$Gestazione,df_splitted$`1`$Peso),
                     cor(df_splitted$`1`$Lunghezza,df_splitted$`1`$Peso),
                     cor(df_splitted$`1`$Cranio,df_splitted$`1`$Peso))
)

##      Colonna Correlazione0 Correlazione1
## 1 Gestazione     0.5991131     0.3884366
## 2  Lunghezza     0.7950173     0.8306399
## 3       Peso     0.7067947     0.6459086

Tranne che per la lunghezza, la correlazione con il peso diminuisce sulle variabili quantitative in esame e la variazione è più evidente con la colonna Gestazione, il che ci porta a considerare nel modello un’interazione tra le due colonne che spieghi anche la non linearità vista nella sezione dedicata a Gestazione.

Tipo di parto

La colonna Tipo.parto indica se il parto è avvenuto naturalmente o per mezzo di taglio cesareo. La cella seguente ne mostra la distribuzione nel dataset.

barplot(table(df$Tipo.parto))

Anche qui osserviamo una disparità tra i due valori, anche se meno netta rispetto alla colonna Fumatrici. Le tre celle che seguono esplorano come il tipo di parto influenzi le relazioni tra il peso, la gestazione e le misure del neonato.

ggplot(df)+
  geom_point(aes(x=Gestazione, y=Peso))+
  facet_grid(~Tipo.parto)

ggplot(df)+
  geom_point(aes(x=Lunghezza, y=Peso))+
  facet_grid(~Tipo.parto)

ggplot(df)+
  geom_point(aes(x=Cranio, y=Peso))+
  facet_grid(~Tipo.parto)

Mentre le misure antropometriche non sembrano influenzate dal tipo di parto, nei grafici relativi alla gestazione la relazione non lineare con il peso è più evidente tra i parti naturali. La cella che segue mostra come cambiano le correlazioni in base al tipo di parto.

df_splitted <- split(df, df$Tipo.parto)

data.frame(
  Colonna = c("Gestazione","Lunghezza","Peso"),
  CorrelazioneNat = c(cor(df_splitted$Nat$Gestazione,df_splitted$Nat$Peso),
                     cor(df_splitted$Nat$Lunghezza,df_splitted$Nat$Peso),
                     cor(df_splitted$Nat$Cranio,df_splitted$Nat$Peso)),
  CorrelazioneCes = c(cor(df_splitted$Ces$Gestazione,df_splitted$Ces$Peso),
                     cor(df_splitted$Ces$Lunghezza,df_splitted$Ces$Peso),
                     cor(df_splitted$Ces$Cranio,df_splitted$Ces$Peso))
)

##      Colonna CorrelazioneNat CorrelazioneCes
## 1 Gestazione       0.6096860       0.5387185
## 2  Lunghezza       0.8011037       0.7837514
## 3       Peso       0.7171237       0.6687096

Anche in questo caso, la correlazione tra Gestazione e Peso è quella che varia più sensibilmente rispetto alle altre variabili, anche se non come nel caso del fumo materno. Inoltre, in tutte e tre le colonne c’è una diminuzione della correlazione con il peso in caso di parto cesareo, il che porta a ipotizzare un’infuenza negativa del tipo di parto sul peso.

Ospedale

La colonna Ospedale indica in quale ospedale del circuito in cui opera l’azienda committente sono nati i bambini. Anche se il risultato può sembrare banale dal punto di vista intuitivo, verifichiamo comunque che l’ospedale non incida sul peso.

Prima di tutto, visualizziamo la distribuzione dei dati.

barplot(table(df$Ospedale))

I tre ospedali sono rappresentati in modo pressoché equo all’interno del dataset. I grafici seguenti mostrano le relazioni tra gestazione e misure con il peso in base all’ospedale in cui si è svolto il parto.

ggplot(df)+
  geom_point(aes(x=Gestazione, y=Peso))+
  facet_grid(~Ospedale)

ggplot(df)+
  geom_point(aes(x=Lunghezza, y=Peso))+
  facet_grid(~Ospedale)

ggplot(df)+
  geom_point(aes(x=Cranio, y=Peso))+
  facet_grid(~Ospedale)

I tre grafici mostrati per ogni coppia di variabili sono pressoché identici, fatto salvo per alcune differenze dovute soprattutto alla casualità. Verifichiamolo ulteriormente con le correlazioni.

df_splitted <- split(df, df$Ospedale)

data.frame(
  Colonna = c("Gestazione","Lunghezza","Peso"),
  CorrelazioneOsp1 = c(cor(df_splitted$osp1$Gestazione,df_splitted$osp1$Peso),
                     cor(df_splitted$osp1$Lunghezza,df_splitted$osp1$Peso),
                     cor(df_splitted$osp1$Cranio,df_splitted$osp1$Peso)),
  CorrelazioneOsp2 = c(cor(df_splitted$osp2$Gestazione,df_splitted$osp2$Peso),
                     cor(df_splitted$osp2$Lunghezza,df_splitted$osp2$Peso),
                     cor(df_splitted$osp2$Cranio,df_splitted$osp2$Peso)),
  CorrelazioneOsp3 = c(cor(df_splitted$osp3$Gestazione,df_splitted$osp3$Peso),
                     cor(df_splitted$osp3$Lunghezza,df_splitted$osp3$Peso),
                     cor(df_splitted$osp3$Cranio,df_splitted$osp3$Peso))
)

##      Colonna CorrelazioneOsp1 CorrelazioneOsp2 CorrelazioneOsp3
## 1 Gestazione        0.6445638        0.5501260        0.5785392
## 2  Lunghezza        0.7988384        0.8105646        0.7811182
## 3       Peso        0.7064814        0.7171440        0.6909397

Le correlazioni con il peso variano in modo minimo tra i tre ospedali, pertanto possiamo escludere definitivamente l’ospedale di nascita dal modello per prevedere il peso dei neonati.

Sesso

Concludiamo l’analisi descrittiva con la colonna Sesso, che indica con un carattere il sesso del neonato. Visualizziamo la distribuzione nella cella che segue.

barplot(table(df$Sesso))

Il dataset è distribuito in modo quasi uniforme tra i due sessi. Utilizziamo questa colonna per filtrare i grafici a dispersione tra il peso, la gestazione e le misure antropometriche.

ggplot(df)+
  geom_point(aes(x=Gestazione, y=Peso))+
  facet_grid(~Sesso)

ggplot(df)+
  geom_point(aes(x=Lunghezza, y=Peso))+
  facet_grid(~Sesso)

ggplot(df)+
  geom_point(aes(x=Cranio, y=Peso))+
  facet_grid(~Sesso)

A parità di regressori, osserviamo che i maschi tendono a concentrarsi in fasce di peso più elevate rispetto alle femmine. Inoltre, La maggioranza dei maschi si concentra nelle fasce di lunghezza e diametro del cranio più alte rispetto alle femmine. Da questo ipotizziamo che il sesso potrebbe non influire sulla gestazione, ma piuttosto sulle misure antropometriche.

Verifichiamo quanto ipotizzato con le variazioni di correlazione tra le colonne calcolate nella cella che segue.

df_splitted <- split(df, df$Sesso)

data.frame(
  Colonna = c("Gestazione","Lunghezza","Peso"),
  CorrelazioneM = c(cor(df_splitted$M$Gestazione,df_splitted$M$Peso),
                     cor(df_splitted$M$Lunghezza,df_splitted$M$Peso),
                     cor(df_splitted$M$Cranio,df_splitted$M$Peso)),
  CorrelazioneF = c(cor(df_splitted$F$Gestazione,df_splitted$F$Peso),
                     cor(df_splitted$F$Lunghezza,df_splitted$F$Peso),
                     cor(df_splitted$F$Cranio,df_splitted$F$Peso))
)

##      Colonna CorrelazioneM CorrelazioneF
## 1 Gestazione     0.5575226     0.6022855
## 2  Lunghezza     0.7855003     0.7901206
## 3       Peso     0.6914895     0.7021958

Le variazioni di correlazione confermano a grandi linee quanto visto nei grafici a dispersione, ossia un’influenza del sesso del neonato sulle sue misure. Questa ipotesi sarà esplorata nel dettaglio nella prossima sezione.

Conclusioni sulle variabili di controllo

Concludiamo l’analisi descrittiva del dataset assumendo quanto segue sul modello e sulle variabili di controllo:

• il fumo materno interagisce con la gestazione abbassando il peso del neonato a parità di settimane di gravidanza;
• il tipo di parto influisce, anche se in minima parte, sulla relazione tra gestazione e peso;
• il sesso del neonato influisce sulle misure antropometriche, che sono più elevati tra i maschi con la stessa durata del parto;
• l’ospedale in cui si svolge il parto non ha alcuna influenza sul peso del neonato.

Nella costruzione del modello terremo dunque in considerazione la colonna Fumatrici per la sua relazione con Gestazione e la colonna Sesso per la sua influenza su Lunghezza e Cranio.

I test di ipotesi

Questa sezione è dedicata all’esecuzione di test statistici per la verifica delle seguenti ipotesi:

• l’associazione tra ospedale e tipo di parto;
• la bontà della stima della media della popolazione di lunghezza e peso dei neonati da parte del campione;
• l’associazione tra sesso e misure di lunghezza, cranio e peso.

A ciascuna di queste ipotesi sarà dedicata una sottosezione di questo documento.

Associazione tra ospedale e tipo di parto

Questa sottosezione verifica l’ipotesi per cui in certi ospedali si preferisca praticare un determinato tipo di parto. Questa ipotesi equivale a quella per cui esista una possibile associazione tra ospedale e tipo di parto.

Il primo passo è stato ottenere una tabella di contingenza dalle colonne Ospedale e Tipo.parto. I risultati sono riepilogati nel balloon plot della cella seguente.

osservazioni <- table(df$Ospedale,df$Tipo.parto)

ggpubr::ggballoonplot(data = as.data.frame(osservazioni))

Osserviamo che la prevalenza del parto naturale rispetto al cesareo si mantiene per tutti e tre gli ospedali e che, anzi, nell’ospedale 3 i parti cesarei sono addirittura meno frequenti rispetto agli altri due ospedali. Graficamente, quindi, possiamo supporre che non esiste associazione tra le due variabili in esame.

Verifichiamo numericamente questa supposizione costruendo la tabella delle osservazioni attese. Osserviamo che il balloon plot ottenuto è simile a quello tratto dalla tabella di contingenza.

attese <- osservazioni
for(i in 1:nrow(osservazioni)){
  for(j in 1:ncol(osservazioni)){
    attese[i,j] <- table(df$Ospedale)[i]*table(df$Tipo.parto)[j]/n
  }
}

ggpubr::ggballoonplot(data = as.data.frame(attese))

La cella che segue calcola il valore della statistica \(X^2\) di Pearson sulla tabella di contingenza e sulle frequenze attese.

(X2 <- sum((osservazioni-attese)^2/attese))

## [1] 1.097202

Il grafico seguente mostra dove si posiziona questo valore (in verde) nella distribuzione \(\chi^2\) con i gradi di libertà ottenuti dal numero di righe e di colonne della tabella di contingenza. La retta tratteggiata di rosso delimita la zona del grafico oltre la quale rifiutiamo l’ipotesi nulla.

plot(density(rchisq(n, df=(ncol(osservazioni)-1)*(nrow(osservazioni)-1))),main = "Relazione tra ospedale e tipo di parto")
abline(v=qchisq(0.95, df=(ncol(osservazioni)-1)*(nrow(osservazioni)-1)),col="red",lty=2)
abline(v=X2, col="green")

Osserviamo che il valore ottenuto ricade nella zona di accettazione, pertanto possiamo accettare l’ipotesi nulla per cui Ospedale e Tipo.parto sono due variabili indipendenti. I risultati ottenuti sono confermati dal test \(\chi^2\) nella cella seguente, dove il p-value supera la significatività al \(5\%\).

chisq.test(osservazioni)

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

Media di lunghezza e peso

In questa sottosezione verificheremo se la media di lunghezza e peso del campione siano una buona stima delle medie della popolazione.

Visualizziamo le prime cinque righe delle colonne Lunghezza e Peso.

head(df[c("Lunghezza","Peso")],5)

##   Lunghezza Peso
## 1       490 3380
## 2       490 3150
## 3       500 3640
## 4       515 3690
## 5       480 3700

Nella cella seguente fissiamo le stime delle due colonne a 50 centimetri e 3300 grammi rispettivamente.

LUNGHEZZA.POP <- 500
PESO.POP <- 3300

Dal momento che non conosciamo la varianza delle due variabili nella popolazione, opteremo per un test T a due code con una significatività del \(5\%\).

La cella seguente effettua il test su Lunghezza e ne memorizza i risultati.

t.lunghezza <- t.test(df$Lunghezza, mu = LUNGHEZZA.POP, conf.level = 0.95, alternative = "two.sided")

Il grafico che segue mostra la posizione dei risultati ottenuti nel test (in verde) all’interno della distribuzione T di Student.

plot(density(rt(n, df = n-1)), main="Differenza tra media della popolazione e media campionaria")
abline(v=c(qt(0.025,df=n-1), qt(0.975,df=n-1)),col="red",lty=2)
abline(v=qt(t.lunghezza$p.value, df=n-1), col="green")

L’assenza della retta verde è segno che il p-value ricade al di fuori della zona di accettazione, il che ci porta a rifiutare l’ipotesi nulla per cui la media di lunghezza del campione è significativamente uguale a 50 centimetri.

In modo analogo, eseguiamo lo stesso tipo di test sulla colonna Peso.

t.peso <- t.test(df$Peso, mu = PESO.POP, conf.level = 0.95, alternative = "two.sided")

Il grafico nella cella seguente mostra graficamente i risultati del test.

plot(density(rt(n, df = n-1)), main = "Differenza tra media della popolazione e media campionaria")
abline(v=c(qt(0.025,df=n-1), qt(0.975,df=n-1)),col="red",lty=2)
abline(v=qt(t.peso$p.value, df=n-1), col="green")

Il p-value del test ricade nella zona di accettazione, perciò non possiamo rifiutare l’ipotesi nulla per cui il campione ha una media significativamente uguale a quella della popolazione.

Sesso e misure antropometriche

In questa sottosezione, riprendiamo l’associazione tra sesso e misure dei neonati ipotizzata nell’analisi descrittiva.

Riportiamo nella cella che segue le prime cinque righe delle colonne Sesso, Peso, Lunghezza e Cranio.

head(df[c("Sesso","Peso","Lunghezza","Cranio")],5)

##   Sesso Peso Lunghezza Cranio
## 1     M 3380       490    325
## 2     F 3150       490    345
## 3     M 3640       500    375
## 4     M 3690       515    365
## 5     F 3700       480    335

Dividiamo il campione in due sottoinsiemi in base al sesso dei neonati.

df_splitted <- split(df[c("Sesso","Peso","Lunghezza","Cranio")], df$Sesso)

Iniziamo dalla colonna Peso ed eseguiamo un test T non appaiato sui due gruppi, dove l’ipotesi nulla è che non ci sia una differenza significativa della media del peso tra i due sessi.

t.test(df_splitted$M$Peso, df_splitted$F$Peso, paired = F, alternative = "two.sided")

## 
##  Welch Two Sample t-test
## 
## data:  df_splitted$M$Peso and df_splitted$F$Peso
## t = 12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  207.0615 287.1051
## sample estimates:
## mean of x mean of y 
##  3408.215  3161.132

Dal p-value prossimo allo 0 non possiamo che rifiutare l’ipotesi nulla a favore dell’alternativa, per cui il peso dei neonati è influenzato dal loro sesso.

In modo analogo, eseguiamo lo stesso tipo di test su Lunghezza e Cranio. I risultati sono riportati nelle due celle seguenti.

t.test(df_splitted$M$Lunghezza, df_splitted$F$Lunghezza, paired = F, alternative = "two.sided")

## 
##  Welch Two Sample t-test
## 
## data:  df_splitted$M$Lunghezza and df_splitted$F$Lunghezza
## t = 9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.876273 11.929470
## sample estimates:
## mean of x mean of y 
##  499.6672  489.7643

t.test(df_splitted$M$Cranio, df_splitted$F$Cranio, paired = F, alternative = "two.sided")

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

I due test, pur essendo eseguiti su colonne diverse, ci portano a concludere che le misure antropometriche dei neonati sono condizionate dal loro sesso e pertanto terremo in considerazione l’interazione tra queste variabili nel modello.

Costruzione e scelta del modello

In questa sezione costruiremo il modello statistico commissionato dall’azienda committente.

Il modello sarà costruito con un approccio a passi in avanti, ossia a partire da un modello con un’unica variabile ne aumenteremo la complessità a ogni passo. Anche alla luce dell’analisi descrittiva svolta in precedenza, le variabili che saranno coinvolte in ordine di priorità sono le seguenti:

• la durata della gravidanza, sia singolarmente che in interazione con il fumo materno;
• la lunghezza e le dimensioni del cranio del neonato, sia singolarmente che come interazione tra loro;
• il sesso, sia singolarmente che in interazione con le misure antropometriche.

Prima di costruire i modelli inizializziamo tre vettori che a ogni modello associno i seguenti valori:

• il coefficiente di determinazione aggiustato;
• il criterio di informazione bayesiano;
• la radice dell’errore quadratico medio.

Radjs <- c()
BICs <- c()
RMSEs <- c()

Iniziamo dal primo modello, dove includiamo la sola variabile Gestazione.

model1 <- lm(Peso~Gestazione, data = df)
summary(model1)

## 
## Call:
## lm(formula = Peso ~ Gestazione, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1603.61  -287.34   -11.07   280.46  1897.75 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3197.250    176.851  -18.08   <2e-16 ***
## Gestazione    166.272      4.532   36.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 423.3 on 2498 degrees of freedom
## Multiple R-squared:  0.3502, Adjusted R-squared:  0.3499 
## F-statistic:  1346 on 1 and 2498 DF,  p-value: < 2.2e-16

Radjs[1] <- summary(model1)$adj.r.squared
BICs[1] <- BIC(model1)
RMSEs[1] <- sqrt(mean(model1$residuals^2))

Nel secondo modello aggiorniamo il modello appena ottenuto con la colonna Fumatrici.

model2 <- update(model1, ~.*Fumatrici)
summary(model2)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Gestazione:Fumatrici, 
##     data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1609.03  -289.03   -11.54   280.97  1898.44 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -3240.771    178.861 -18.119   <2e-16 ***
## Gestazione             167.495      4.585  36.534   <2e-16 ***
## Fumatrici             1343.426   1164.746   1.153    0.249    
## Gestazione:Fumatrici   -36.764     29.646  -1.240    0.215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 422.9 on 2496 degrees of freedom
## Multiple R-squared:  0.352,  Adjusted R-squared:  0.3513 
## F-statistic:   452 on 3 and 2496 DF,  p-value: < 2.2e-16

Radjs[2] <- summary(model2)$adj.r.squared
BICs[2] <- BIC(model2)
RMSEs[2] <- sqrt(mean(model2$residuals^2))

Il terzo modello aggiorna il secondo con le colonne Lunghezza e Cranio singolarmente.

model3 <- update(model2, ~.+Lunghezza+Cranio)
summary(model3)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Gestazione:Fumatrici, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1105.74  -184.29   -13.72   167.14  2620.98 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6789.3469   136.7107 -49.662   <2e-16 ***
## Gestazione              32.2388     3.8632   8.345   <2e-16 ***
## Fumatrici              521.3624   765.0499   0.681    0.496    
## Lunghezza               10.4112     0.3024  34.433   <2e-16 ***
## Cranio                  10.7854     0.4283  25.180   <2e-16 ***
## Gestazione:Fumatrici   -13.8705    19.4735  -0.712    0.476    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.7 on 2494 degrees of freedom
## Multiple R-squared:  0.7208, Adjusted R-squared:  0.7202 
## F-statistic:  1288 on 5 and 2494 DF,  p-value: < 2.2e-16

Radjs[3] <- summary(model3)$adj.r.squared
BICs[3] <- BIC(model3)
RMSEs[3] <- sqrt(mean(model3$residuals^2))

Dal terzo modello otteniamo tre aggiornamenti:

• il quarto modello aggiunge l’interazione tra Lunghezza e Cranio;

• il quinto modello aggiunge Sesso singolarmente;

• il sesto modello aggiunge l’interazione tra Sesso e le misure antropometriche.

model4 <- update(model3, ~.+(Lunghezza*Cranio))
summary(model4)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Gestazione:Fumatrici + Lunghezza:Cranio, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1108.60  -182.60    -9.35   165.38  2873.42 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.345e+03  1.027e+03  -1.310   0.1905    
## Gestazione            3.873e+01  4.029e+00   9.613  < 2e-16 ***
## Fumatrici             7.603e+02  7.622e+02   0.998   0.3186    
## Lunghezza            -1.371e+00  2.223e+00  -0.617   0.5375    
## Cranio               -6.274e+00  3.217e+00  -1.950   0.0513 .  
## Gestazione:Fumatrici -1.990e+01  1.940e+01  -1.026   0.3052    
## Lunghezza:Cranio      3.521e-02  6.581e-03   5.350  9.6e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 276.2 on 2493 degrees of freedom
## Multiple R-squared:  0.7239, Adjusted R-squared:  0.7233 
## F-statistic:  1090 on 6 and 2493 DF,  p-value: < 2.2e-16

Radjs[4] <- summary(model4)$adj.r.squared
BICs[4] <- BIC(model4)
RMSEs[4] <- sqrt(mean(model4$residuals^2))

model5 <- update(model3, ~.+Sesso)
summary(model5)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Sesso + Gestazione:Fumatrici, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1138.31  -182.30   -18.44   164.51  2624.11 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6667.9088   136.4437 -48.869  < 2e-16 ***
## Gestazione              32.0177     3.8254   8.370  < 2e-16 ***
## Fumatrici              777.3046   758.4089   1.025    0.306    
## Lunghezza               10.1865     0.3011  33.836  < 2e-16 ***
## Cranio                  10.6641     0.4245  25.123  < 2e-16 ***
## SessoM                  79.8505    11.2261   7.113 1.48e-12 ***
## Gestazione:Fumatrici   -20.4706    19.3051  -1.060    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2493 degrees of freedom
## Multiple R-squared:  0.7263, Adjusted R-squared:  0.7257 
## F-statistic:  1103 on 6 and 2493 DF,  p-value: < 2.2e-16

Radjs[5] <- summary(model5)$adj.r.squared
BICs[5] <- BIC(model5)
RMSEs[5] <- sqrt(mean(model5$residuals^2))

model6 <- update(model3, ~.+Sesso*(Lunghezza+Cranio))
summary(model6)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Sesso + Gestazione:Fumatrici + Lunghezza:Sesso + Cranio:Sesso, 
##     data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.35  -182.72   -19.27   163.64  2554.82 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6500.0268   175.2122 -37.098   <2e-16 ***
## Gestazione              32.2167     3.8279   8.416   <2e-16 ***
## Fumatrici              727.4740   758.5634   0.959   0.3376    
## Lunghezza                9.7985     0.3793  25.835   <2e-16 ***
## Cranio                  10.7071     0.5807  18.438   <2e-16 ***
## SessoM                -315.6436   250.2131  -1.261   0.2072    
## Gestazione:Fumatrici   -19.2286    19.3085  -0.996   0.3194    
## Lunghezza:SessoM         0.8983     0.5338   1.683   0.0925 .  
## Cranio:SessoM           -0.1457     0.8430  -0.173   0.8628    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7267, Adjusted R-squared:  0.7259 
## F-statistic: 828.1 on 8 and 2491 DF,  p-value: < 2.2e-16

Radjs[6] <- summary(model6)$adj.r.squared
BICs[6] <- BIC(model6)
RMSEs[6] <- sqrt(mean(model6$residuals^2))

Il grafico seguente mostra come varia il coefficiente di determinazione aggiustato in base al modello. La linea tratteggiata verde intercetta il modello che ha ottenuto il massimo valore.

plot(Radjs)
abline(h=max(Radjs), col="green",lty=2)

In modo analogo, il grafico che segue mostra i valori di criterio di informazione bayesiano intercettando il valore minimo.

plot(BICs)
abline(h=min(BICs), col="green", lty=2)

Il grafico che segue, invece, mostra i valori di radice di errore quadratico medio.

plot(RMSEs)
abline(h=min(RMSEs), col="green", lty=2)

I risultati ottenuti dai modelli ci portano a scegliere come candidato ideale il quinto modello. Questo modello infatti è quello che ha ottenuto il massimo in termini di bontà dell’adattamento, è il meno complesso e minimizza l’errore quadratico medio tra i modelli ottenuti.

Analisi dei residui del modello

In questa sezione ci occuperemo di analizzare la qualità del modello ottenuto nella precedente sezione dal punto di vista dei suoi residui.

Per semplicità, memorizziamo questo modello come final.model.

final.model <- model5

Estraiamo i residui del modello sui dati e visualizziamone la densità in un grafico a linee. La linea grigia tratteggiata rappresenta la distribuzione normale a parità di media e deviazione standard.

residui=final.model$residuals

plot(density(residui))
lines(density(rnorm(n, mean=mean(residui), sd=sd(residui))),col="lightgrey", lty=2)

Visivamente, il grafico della densità appare più leptocurtica e leggermente asimmetrica rispetto alla normale. Per stabilire se i residui hanno distribuzione normale, tuttavia, ci affideremo a un test di normalità di Shapiro-Wilk.

shapiro.test(residui)

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

Testiamo l’omoschedasticità tramite il test di Breusch-Pagan.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(final.model)

## 
##  studentized Breusch-Pagan test
## 
## data:  final.model
## BP = 89.347, df = 6, p-value < 2.2e-16

Infine, tramite il test di Durbin-Watson, verifichiamo l’assenza di correlazione tra i residui e i regressori del modello.

dwtest(final.model)

## 
##  Durbin-Watson test
## 
## data:  final.model
## DW = 1.9568, p-value = 0.14
## alternative hypothesis: true autocorrelation is greater than 0

I risultati di questi test ci dicono che, malgrado la bontà dell’adattamento e l’assenza di correlazione con i regressori, i residui non sono distribuiti normalmente e non hanno varianza costante. Questo può essere dovuto alla presenza di valori anomali tra i dati e, di conseguenza, tra i residui.

Verifichiamo la presenza di questi valori prima graficamente e poi esaminandoli nel dettaglio dal dataset. Il grafico che segue evidenzia i punti di leva dei residui (sopra la linea rossa).

plot(hatvalues(final.model))
abline(h=2*sum(hatvalues(final.model))/n, col="red")

Il grafico che segue, invece, delimita gli outliers al di fuori della fascia delimitata dalle due rette rosse.

plot(rstudent(final.model))
abline(h=c(-2,2), col="red")

Infine, esaminiamo i residui dal punto di vista della distanza di Cook.

plot(cooks.distance(final.model))
abline(h=0.5, col="red")

In tutti e tre i grafici emerge la presenza di un punto particolarmente estremo, insieme a una notevole quantità di valori anomali di ogni tipo.

Per analizzare questi punti, estraiamo dal dataset i record corrispondenti e le colonne del modello.

leverage_pts <- which(hatvalues(final.model) > 2*sum(hatvalues(final.model))/n)
outlier_pts <- which(rstudent(final.model)< -2 | rstudent(final.model)>2)
cook_pts <- which(cooks.distance(final.model)>0.5)

strange_vals <- unique(c(leverage_pts, outlier_pts, cook_pts))

df_strange <- df[strange_vals, c("Peso","Gestazione","Fumatrici","Lunghezza","Cranio","Sesso")]

dim(df_strange)

## [1] 259   6

I valori anomali riguardano 259 neonati, di cui di seguito riportiamo graficamente la distribuzione delle colonne coinvolte.

Peso100 <- df_strange$Peso/100

min_limit <- floor(min(Peso100))
max_limit <- ceiling(max(Peso100))
Peso_CL <- cut(Peso100, breaks = seq(min_limit, max_limit, 2.50))

barplot(table(Peso_CL), las=2)

barplot(table(df_strange$Gestazione))

barplot(table(df_strange$Fumatrici))

min_limit <- round(min(df_strange$Lunghezza), digits = -2)
max_limit <- round(max(df_strange$Lunghezza), digits = -2)
Lunghezza_CL <- cut(df_strange$Lunghezza, breaks = seq(min_limit, max_limit, 10))

barplot(table(Lunghezza_CL), las=2)

min_limit <- round(min(df_strange$Cranio), digits = -2)
max_limit <- round(max(df_strange$Cranio), digits = -2)
Cranio_CL <- cut(df_strange$Cranio, breaks = seq(min_limit, max_limit, 10))

barplot(table(Cranio_CL), las=2)

barplot(table(df_strange$Sesso))

Anche alla luce di quanto emerso dall’analisi descrittiva, osserviamo che il modello ha risentito dell’asimmetria dei valori di gestazione, fumo materno e misure nel dataset. Questo ci permette di identificare come anomali per il modello neonati nati da donne fumatrici prima della 35° settimana, con una lunghezza al di sotto del 44 centimetri e con una circonferenza cranica sotto i 30 centimetri.

Previsioni del modello

In questa sezione ci occuperemo di effettuare una previsione tramite il modello ottenuto.

Su richiesta dell’azienda committente, ci occuperemo di prevedere il peso di una neonata alla 39° settimana la cui madre ha già avuto due gravidanze. Useremo questi parametri per estrarre nascite simili dal dataset originale e considereremo in particolare le colonne prese in considerazione dal modello.

df_filtered <- df[df$Sesso=="F"&df$Gestazione==39&df$N.gravidanze==2,c("Peso","Gestazione","Fumatrici","Lunghezza","Cranio","Sesso")]

df_filtered

##      Peso Gestazione Fumatrici Lunghezza Cranio Sesso
## 2    3150         39         0       490    345     F
## 25   2630         39         0       465    320     F
## 81   3090         39         0       475    340     F
## 103  3100         39         0       500    320     F
## 175  3750         39         0       520    353     F
## 194  3530         39         1       505    332     F
## 297  2750         39         0       470    320     F
## 380  3170         39         0       495    335     F
## 443  2910         39         0       480    328     F
## 452  3800         39         0       500    346     F
## 512  3120         39         0       500    335     F
## 540  3050         39         0       500    305     F
## 548  3250         39         0       480    340     F
## 580  2990         39         0       490    335     F
## 610  3220         39         0       515    346     F
## 710  3330         39         0       505    350     F
## 890  3660         39         0       490    318     F
## 903  4100         39         0       520    352     F
## 973  3340         39         0       510    345     F
## 1058 3220         39         0       475    320     F
## 1128 3580         39         0       500    341     F
## 1138 3220         39         0       480    335     F
## 1258 3440         39         0       500    332     F
## 1488 2260         39         0       460    330     F
## 1515 2850         39         0       495    344     F
## 1516 2850         39         0       480    330     F
## 1521 3090         39         0       500    325     F
## 1627 3330         39         0       480    350     F
## 1635 3430         39         0       445    322     F
## 1700 3210         39         0       500    340     F
## 1778 3180         39         0       470    325     F
## 1876 3180         39         0       510    340     F
## 1887 2850         39         0       500    340     F
## 1973 3080         39         0       480    332     F
## 1984 3170         39         0       510    332     F
## 2043 3240         39         0       490    342     F
## 2102 3550         39         0       495    362     F
## 2139 3440         39         0       505    335     F
## 2226 3540         39         0       500    354     F
## 2240 3250         39         0       490    336     F
## 2414 3740         39         0       490    345     F
## 2455 3250         39         0       500    330     F

Per le previsioni del modello sui nuovi dati, assumeremo quanto segue:

• I dati di sesso e gestazione del neonato rimarranno quelli forniti dall’azienda committente;

• Per il fumo materno saranno presi in esame sia il caso di madre fumatrice che quello di madre non fumatrice;

• Per le misure di peso, lunghezza e cranio calcoleremo la media arrotondata all’intero delle colonne Peso, Lunghezza e Cranio condizionate dalla colonna Fumatrici.

I dati così descritti saranno costruiti e visualizzati nella cella che segue.

Peso <- c(round(mean(df_filtered[df_filtered$Fumatrici==0,]$Peso)),
          round(mean(df_filtered[df_filtered$Fumatrici==1,]$Peso)))
Lunghezza <- c(round(mean(df_filtered[df_filtered$Fumatrici==0,]$Lunghezza)),
          round(mean(df_filtered[df_filtered$Fumatrici==1,]$Lunghezza)))
Cranio <- c(round(mean(df_filtered[df_filtered$Fumatrici==0,]$Cranio)),
          round(mean(df_filtered[df_filtered$Fumatrici==1,]$Cranio)))

df_to_predict <- data.frame(
  Sesso = rep("F",2),
  Gestazione = rep(39, 2),
  Peso = Peso,
  Lunghezza=Lunghezza,
  Cranio = Cranio,
  Fumatrici = c(0,1)
)

df_to_predict

##   Sesso Gestazione Peso Lunghezza Cranio Fumatrici
## 1     F         39 3228       492    336         0
## 2     F         39 3530       505    332         1

A questi dati aggiungiamo la colonna Peso.Pred per le previsioni del modello.

df_to_predict$Peso.Pred <- predict(final.model, df_to_predict)

df_to_predict

##   Sesso Gestazione Peso Lunghezza Cranio Fumatrici Peso.Pred
## 1     F         39 3228       492    336         0  3175.690
## 2     F         39 3530       505    332         1  3244.408

Effettuiamo ora un test statistico per stabilire se le previsioni del modello possano essere una buona stima del peso. Assumendo che la situazione in esame sia equiparabile a una stima della media del peso della popolazione di neonate alla 39° settimana, eseguiremo un test T a due code con media stimata pari al peso predetto e campione dato dal dataset filtrato.

Ciò premesso, la colonna PValue indicherà per i due record il p-value ottenuto dal test.

df_to_predict$PValue <- c(t.test(df_filtered$Peso, mu=df_to_predict$Peso.Pred[1], alternative = "two.sided")$p.value,
                          t.test(df_filtered$Peso, mu=df_to_predict$Peso.Pred[2], alternative = "two.sided")$p.value)

df_to_predict

##   Sesso Gestazione Peso Lunghezza Cranio Fumatrici Peso.Pred    PValue
## 1     F         39 3228       492    336         0  3175.690 0.2558641
## 2     F         39 3530       505    332         1  3244.408 0.8641681

Entrambi i p-value superano la significatività al \(5\%\), pertanto non rifiutiamo l’ipotesi nulla e assumeremo che le previsioni ottenute siano stime plausibili del peso dei neonati.

Analisi del modello

Concludiamo il progetto con un’analisi del modello ottenuto, di cui ne riepiloghiamo le caratteristiche nella cella che segue.

summary(final.model)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Fumatrici + Lunghezza + Cranio + 
##     Sesso + Gestazione:Fumatrici, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1138.31  -182.30   -18.44   164.51  2624.11 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6667.9088   136.4437 -48.869  < 2e-16 ***
## Gestazione              32.0177     3.8254   8.370  < 2e-16 ***
## Fumatrici              777.3046   758.4089   1.025    0.306    
## Lunghezza               10.1865     0.3011  33.836  < 2e-16 ***
## Cranio                  10.6641     0.4245  25.123  < 2e-16 ***
## SessoM                  79.8505    11.2261   7.113 1.48e-12 ***
## Gestazione:Fumatrici   -20.4706    19.3051  -1.060    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2493 degrees of freedom
## Multiple R-squared:  0.7263, Adjusted R-squared:  0.7257 
## F-statistic:  1103 on 6 and 2493 DF,  p-value: < 2.2e-16

In linea di massima, il modello ottenuto può essere interpretato come segue:

• I regressori relativi al fumo (Fumatrici e Gestazione:Fumatrici) sembrano i meno significativi per il modello, pur rappresentando il maggiore effetto positivo e negativo rispettivamente. A parità degli altri regressori, il figlio di una madre fumatrice guadagnerebbe circa 777 grammi, ma perderebbe 20 grammi per ogni settimana di gestazione. In altre parole, in termini di peso del neonato, la gestazione di una fumatrice equivale a sottrarre una settimana alla gestazione di una non fumatrice;

• Il maggior effetto positivo tra i regressori significativi è dato dalle settimane di gestazione, dove un neonato guadagna circa 32 grammi di peso per ogni settimana;

• La lunghezza e il cranio hanno pressapoco lo stesso effetto sul peso in termini di dimensione e di segno;

• A parità di fumo materno, gestazione e misure antropometriche, un maschio pesa circa 80 grammi in più rispetto a una femmina.

Utilizziamo il modello per effettuare delle previsioni sul peso dei neonati nel dataset fornito dall’azienda committente e visualizziamo le prime cinque righe.

df$Peso.Pred <- predict(final.model, df)

head(df,5)

##   Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1         26            0         0         42 3380       490    325        Nat
## 2         21            2         0         39 3150       490    345        Nat
## 3         34            3         0         38 3640       500    375        Nat
## 4         28            1         0         41 3690       515    365        Nat
## 5         20            0         0         38 3700       480    335        Nat
##   Ospedale Sesso Peso.Pred
## 1     osp3     M  3213.915
## 2     osp1     F  3251.294
## 3     osp2     M  3720.915
## 4     osp2     M  3863.125
## 5     osp3     F  3010.770

Il grafico che segue mostra dove si posizionano le previsioni del modello rispetto alla performance ideale in cui il peso predetto coincide con il peso reale.

ggplot(df)+
  geom_point(aes(x=Peso, y=Peso.Pred))+
  geom_line(aes(x=Peso, y=Peso), color="red")

In generale, il modello sovrastima i pesi dei neonati nel dataset, pur non mancando range di peso in cui questa tendenza si inverte come al di sotto dei 1500 grammi o sopra i 4000 grammi.

Per esaminare la situazione nel dettaglio, aggiungeremo una colonna Stima in cui distingueremo i record del dataset a seconda che si sia verificata una sovrastima o una sottostima. Nella stessa cella, visualizzeremo le prime cinque righe del dataset così modificato.

df$Stima = ifelse(df$Peso.Pred > df$Peso, "Over", "Under")

head(df,5)

##   Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1         26            0         0         42 3380       490    325        Nat
## 2         21            2         0         39 3150       490    345        Nat
## 3         34            3         0         38 3640       500    375        Nat
## 4         28            1         0         41 3690       515    365        Nat
## 5         20            0         0         38 3700       480    335        Nat
##   Ospedale Sesso Peso.Pred Stima
## 1     osp3     M  3213.915 Under
## 2     osp1     F  3251.294  Over
## 3     osp2     M  3720.915  Over
## 4     osp2     M  3863.125  Over
## 5     osp3     F  3010.770 Under

Il grafico a barre che segue mostra la distribuzione della colonna Stima.

barplot(table(df$Stima))

Il grafico a barre conferma la tendenza alla sovrastima del modello. Entreremo ulteriormente nel dettaglio nelle prossime sottosezioni, dove esamineremo il modello dal punto di vista dei singoli regressori.

Gestazione e fumo materno

Dedichiamo questa sottosezione alla gestazione e al fumo materno.

Il grafico che segue mostra i pesi reali del dataset (in nero) e li confronta con i pesi stimati dal modello (in rosso). La linea verde chiaro rappresenta la relazione tra peso e gestazione individuata dal modello.

ggplot(df)+
  geom_point(aes(x=Gestazione,y=Peso))+
  geom_point(aes(x=Gestazione,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Gestazione,y=Peso.Pred),color="lightgreen")

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Si nota che le fasce di peso per settimane di gestazione tendono a restringersi intorno alla linea individuata dal modello. Questa linea, poi, presenta una pendenza più elevata prima della 37° settimana. Visualizziamo lo stesso grafico, ma condizionandolo con il fumo materno.

ggplot(df)+
  geom_point(aes(x=Gestazione,y=Peso))+
  geom_point(aes(x=Gestazione,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Gestazione,y=Peso.Pred),color="lightgreen")+
  facet_grid(Fumatrici~.)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Failed to fit group -1.
## Caused by error in `smooth.construct.cr.smooth.spec()`:
## ! x has insufficient unique values to support 10 knots: reduce k.

La tendenza generale vista nel grafico precedente si mantiene soprattutto tra le madri non fumatrici, mentre non è possibile dire lo stesso delle madri fumatrici anche a causa dello sbilanciamento dei dati.

Riesaminiamo la situazione, ma in base al valore della colonna Stima.

ggplot(df)+
  geom_point(aes(x=Gestazione,y=Peso))+
  geom_point(aes(x=Gestazione,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Gestazione,y=Peso.Pred),color="lightgreen")+
  facet_grid(~Stima)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Osserviamo che mentre la relazione rimane pressoché identica nelle sottostime, in caso di sovrastima si ha un appiattimento della relazione tra peso e gestazione prima della 32° settimana.

Lunghezza

In questa sottosezione esamineremo il modello dal punto di vista della lunghezza del neonato, di cui riepiloghiamo la situazione generale.

ggplot(df)+
  geom_point(aes(x=Lunghezza,y=Peso))+
  geom_point(aes(x=Lunghezza,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Lunghezza,y=Peso.Pred),color="lightgreen")

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

La tendenza generale è un range più stretto tra i pesi previsti rispetto a quelli reali. Inoltre, la relazione individuata dal modello mostra una pendenza minore sotto i 35 cm.

Il grafico che segue mostra gli stessi dati in base alla colonna Stima.

ggplot(df)+
  geom_point(aes(x=Lunghezza,y=Peso))+
  geom_point(aes(x=Lunghezza,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Lunghezza,y=Peso.Pred),color="lightgreen")+
  facet_grid(~Stima)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Se la relazione del modello è più accentuata nelle sottostime, nelle sovrastime la relazione tra previsione e lunghezza è pù lineare.

Cranio

In questa sottosezione riepiloghiamo la performance del modello dal punto di vista della circonferenza cranica.

ggplot(df)+
  geom_point(aes(x=Cranio,y=Peso))+
  geom_point(aes(x=Cranio,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Cranio,y=Peso.Pred),color="lightgreen")

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

A differenza di Lunghezza, la pendenza della linea che mette in relazione Cranio e previsioni di peso diminuisce dopo i 30 cm per poi appiattirsi dopo i 37 cm.

Di seguito vediamo come si comporta il modello in base alla colonna Stima.

ggplot(df)+
  geom_point(aes(x=Cranio,y=Peso))+
  geom_point(aes(x=Cranio,y=Peso.Pred),color="red")+
  geom_smooth(aes(x=Cranio,y=Peso.Pred),color="lightgreen")+
  facet_grid(~Stima)

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Mentre la relazione del modello si accentua in caso di sottostima, nelle sovrastime si evidenzia una relazione di tipo non lineare tra il peso previsto e la circonferenza del cranio.

Sesso

In questa ultima sottosezione ci occuperemo di confrontare le stime di peso in funzione del sesso. Riprendiamo quindi il grafico che confronta le previsioni con i valori reali.

ggplot(df)+
  geom_point(aes(x=Peso, y=Peso.Pred))+
  geom_line(aes(x=Peso, y=Peso), color="red")

Condizioniamo questo grafico dal punto di vista del sesso.

ggplot(df)+
  geom_point(aes(x=Peso, y=Peso.Pred))+
  geom_line(aes(x=Peso, y=Peso), color="red")+
  facet_grid(~Sesso)

Dal grafico osserviamo che il modello tende a sovrastimare il peso, soprattutto tra i maschi, mentre è riuscito a collegare il sesso maschile alle fasce di peso più elevate. Un quadro più preciso della situazione è dato dal grafico a barre affiancate della cella seguente, che conta i valori della colonna Stima in base a Sesso.

ggplot(data = df)+
  geom_bar(aes(x = Stima, fill=Sesso), position = "dodge" , stat = "count" )

Dal grafico è evidente che il modello sovrastimi il peso più di frequente tra le femmine, mentre le sottostime sono più frequenti tra i maschi.

Conclusioni sul modello

Per riassumere, possiamo riepilogare il modello ottenuto come segue:

• la tendenza generale del modello è la sovrastima del peso rispetto a quello reale;

• la relazione generale tra peso e regressori è più accentuato nelle sottostime rispetto alle sovrastime;

• le sovrastime sono più frequenti tra le femmine rispetto ai maschi.

Conclusioni

Concludiamo il progetto riepilogando i risultati ottenuti nel corso di questo documento.

Il progetto è stato richiesto dal Neonatal Health Solutions e si poneva l’obiettivo di costruire un modello statistico in grado di prevedere il peso dei neonati. Nella prima sezione è stata svolta un’analisi descrittiva dei dati forniti dall’azienda riguardo 2500 neonati nati in 3 ospedali del circuito. Attraverso questa analisi, sono state individuate le possibili variabili più rilevanti per la costruzione del modello.

Dopo una serie di test statistici svolti sui dati, sono stati costruiti diversi modelli con un approccio per step in avanti. Dal confronto di diverse misure, è stato selezionato il modello che permette di prevedere il peso dalla durata della gravidanza, dal fumo materno, dal sesso del neonato e dalle sue misure antropometriche. La qualità del modello è stata ulteriormente valutata dall’analisi dei suoi residui, che tuttavia hanno risentito dei valori anomali presenti nel dataset.

Nonostante la diagnostica dei residui, le previsioni effettuate dal modello sui parametri forniti dall’azienda si sono rivelate stime plausibili del peso dei neonati. Il progetto si è concluso con una panoramica del modello ottenuto, con particolare focus sul tipo di stima ottenuto dal modello sui dati.