Raccolta dei Dati e Struttura del Dataset:

Viene caricato il dataset neonanti.csv e visualizzate le prime dieci righe, per avere un’idea generale della sua struttura.

dati<-read.csv("neonati.csv")
kable(head(dati, 10), format = "html", caption = "Prime 10 righe del dataset") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Prime 10 righe del dataset
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 42 3380 490 325 Nat osp3 M
21 2 0 39 3150 490 345 Nat osp1 F
34 3 0 38 3640 500 375 Nat osp2 M
28 1 0 41 3690 515 365 Nat osp2 M
20 0 0 38 3700 480 335 Nat osp3 F
32 0 0 40 3200 495 340 Nat osp2 F
26 1 0 39 3100 480 345 Nat osp3 F
25 0 0 40 3580 510 349 Nat osp1 M
22 1 0 40 3670 500 335 Ces osp2 F
23 0 0 41 3700 510 362 Ces osp2 F

Allo stesso tempo viene usata la funzione str per conoscere la lunghezza del dataset e il tipo di variabili presenti.

N=dim(dati)[1]
str(dati)
## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : chr  "Nat" "Nat" "Nat" "Nat" ...
##  $ Ospedale    : chr  "osp3" "osp1" "osp2" "osp2" ...
##  $ Sesso       : chr  "M" "F" "M" "M" ...

In questo caso si può osservare la presenza di 10 variabili e di 2500 casi. In particolare si ha:

In prima analisi vengono stampate le statistiche descrittive relative alle variabili quantitative:

vars <- dati %>%
  select(Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio)
tabella<-data.frame()
for (var in names(vars)){
  riga<-data.frame(
    Variabile=var,
    Media = round(mean(vars[[var]]),2),
    Q1 = round(as.numeric(quantile(vars[[var]], 0.25, na.rm = TRUE)),2),
    Mediana = round(median(vars[[var]]),2),
    Q3 = round(as.numeric(quantile(vars[[var]], 0.75, na.rm = TRUE)),2),
    Minimo = round(min(vars[[var]]),2),
    Massimo = round(max(vars[[var]]),2),
    Dev.St = round(sd(vars[[var]]),2),
    Varianza = round(var(vars[[var]]),2),
    Asimmetria = round(skewness(vars[[var]]),2),
    Curtosi = round(kurtosis(vars[[var]]) - 3),2)
tabella<-rbind(tabella,riga)
}

tabella <- tabella[ , -ncol(tabella)]

kable(tabella, format = "html", caption = "Statistiche descrittive variabili quantitative") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche descrittive variabili quantitative
Variabile Media Q1 Mediana Q3 Minimo Massimo Dev.St Varianza Asimmetria Curtosi
Anni.madre 28.16 25 28 32 0 46 5.27 27.81 0.04 0
N.gravidanze 0.98 0 1 1 0 12 1.28 1.64 2.51 11
Gestazione 38.98 38 39 40 25 43 1.87 3.49 -2.07 8
Peso 3284.08 2990 3300 3620 830 4930 525.04 275665.68 -0.65 2
Lunghezza 494.69 480 500 510 310 565 26.32 692.67 -1.51 6
Cranio 340.03 330 340 350 235 390 16.43 269.79 -0.79 3

Osservazione: Tutte le variabili presentano valori di curtosi positivi, si può quindi affermare che sono distribuzioni leptocurtiche, che risultano più allungate rispetto alla normale.

Stupisce il fatto che la variabile Anni.madre abbia come valore minimo 0, questo fa presagire che ci potrebbero essere stati errori durante l’inserimento del dato, in quanto una madre non può avere età pari a 0, per questo si può guardare la frequenza di ciascuna modalità, presente in questa variabile e sostituire agli anni che non sembrano adeguati il valore mediano di 28. Successivamente verrano visualizzate nuovamente le statistiche.

kable(table(dati$Anni.madre), format = "html", caption = "Frequenze assolute di Anni.madre") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Frequenze assolute di Anni.madre
Var1 Freq
0 1
1 1
13 1
14 2
15 6
16 13
17 18
18 24
19 45
20 66
21 74
22 100
23 115
24 131
25 180
26 184
27 197
28 172
29 174
30 200
31 147
32 159
33 110
34 96
35 66
36 64
37 41
38 38
39 27
40 19
41 13
42 8
43 2
44 4
45 1
46 1

Sia il valore 0, che il valore 1 sono impossibili, per questo vengono sostituiti. E ristampate la statistiche descrittive.

dati$Anni.madre[dati$Anni.madre %in% c(0, 1)] <- 28
vars2 <- dati %>%
  select(Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio)
tabella2<-data.frame()
for (var in names(vars2)){
  riga<-data.frame(
    Variabile=var,
    Media = round(mean(vars2[[var]]),2),
    Q1 = round(as.numeric(quantile(vars2[[var]], 0.25, na.rm = TRUE)),2),
    Mediana = round(median(vars2[[var]]),2),
    Q3 = round(as.numeric(quantile(vars2[[var]], 0.75, na.rm = TRUE)),2),
    Minimo = round(min(vars2[[var]]),2),
    Massimo = round(max(vars2[[var]]),2),
    Dev.St = round(sd(vars2[[var]]),2),
    Varianza = round(var(vars2[[var]]),2),
    Asimmetria = round(skewness(vars2[[var]]),2),
    Curtosi = round(kurtosis(vars2[[var]]) - 3),2)
tabella2<-rbind(tabella2,riga)
}

tabella2 <- tabella2[ , -ncol(tabella2)]

kable(tabella2, format = "html", caption = "Statistiche descrittive variabili quantitative") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche descrittive variabili quantitative
Variabile Media Q1 Mediana Q3 Minimo Massimo Dev.St Varianza Asimmetria Curtosi
Anni.madre 28.19 25 28 32 13 46 5.22 27.20 0.15 0
N.gravidanze 0.98 0 1 1 0 12 1.28 1.64 2.51 11
Gestazione 38.98 38 39 40 25 43 1.87 3.49 -2.07 8
Peso 3284.08 2990 3300 3620 830 4930 525.04 275665.68 -0.65 2
Lunghezza 494.69 480 500 510 310 565 26.32 692.67 -1.51 6
Cranio 340.03 330 340 350 235 390 16.43 269.79 -0.79 3

Si può dunque osservare per la variabile Anni.madre che il 50 % dei valori è situato tra i 25 e i 32 anni, con deviazione standard di circa 5 anni. In generale si ha un’asimettria leggermente positiva, quindi avremo la distribuzione leggermente spostata sulla sinistra e quasi simmetrica con valore di curtosi pari a 0, quindi sarà distribuita in modo molto simile alla normale.

La variabile N.gravidanze presenta sicuramente almeno un outlier, in quanto come si può osservare il 75% delle donne ha avuto al massimo un figlio, mentre il valore massimo arriva a 12. La deviazione standard supera la media, questo indica che il numero di gravidanze può variare ampiamente rispetto la media. L’indice di asimmetria è positivo, significa che la maggior parte delle madri ha avuto un numero basso di figli o si trova alla prima gravidanza. SI ha poi una curtosi molto elevata, quindi la distribuzione di tale variabile sarà molto alta, come se ci fosse una concentrazione rispetto ad un particolare valore.

Nel 50% dei casi le settimane di Gestazione si concentrano tra le 38 e le 40. Ammettendo una deviazione standard bassa. Il valore minimo (neonato prematuro) in questo caso è riportato a 25 settimane. Tale variabile presenta una asimettria negativa e alta, significa che in generale le donne hanno periodi di gestazione più lunghi e ci sono però alcuni casi significativi di prematuri. Anche in questo caso la distribuzione risulta leptocurtica e quindi più allungata rispetto alla normale.

Per qunato riguarda il Peso, si può osservare una distribuzione leggermente asimettrica negativa, quindi spostata verso destra, con un peso mediano di 3300 grammi e lievemente più allungata rispetto la normale. La deviazione standard è di 525 grammi, non è una grande variazione, si può controllare se sia dovuto al sesso del neonato.

La Lunghezza presenta un valore medio di 495 mm, e l’IQR, ossia il 50% dei casi si trova tra 480 e 510 mm, in generale la lunghezza varia di 26 mm rispetto alla media, quindi una variazione minima. L’asimmetria è negativa, quindi la distribuzione è spostata verso valori più alti, mentre anche in questo caso la curtosi positiva, mostra un valore alto positivo, che indica una distribuzione più appiattita rispetto alla normale.

La variabile Cranio contenente la lunghezza del diametro del cranio dei neonati, presenta asimettria leggermente negativa, quindi tenderà a concentrarsi verso destra con una lunga verso sinistra La deviazione standard rimane bassa ed il 50 % dei valori si concentra tra 330 e 350 millimetri.

L’utilizzo di bar plot per le variabili discrete, e delle distribuzioni per quelle continue, permette di visualizzare graficamente quanto affermato.

p1<-ggplot(data=dati)+
  geom_bar(aes(x=N.gravidanze),
           stat='count',
           col='black',
           fill='lightgreen')+
  labs(x="Numero di gravidanze",
       y="Frequenze assolute delle madri")+
  theme_classic()+
  theme(legend.position="bottom")

p2<-ggplot(data=dati)+
  geom_bar(aes(x=Gestazione),
           stat='count',
           col='black',
           fill='lightgreen')+
  labs(x="Settimane di gestazione",
       y="Frequenze assolute delle madri")+
  theme_classic()+
  theme(legend.position="bottom")

p3<-ggplot(data=dati)+
  geom_bar(aes(x=Anni.madre),
           stat='count',
           col='black',
           fill='lightgreen')+
  labs(x="Età madri",
       y="Totale")+
  theme_classic()+
  theme(legend.position="bottom")
p1+p2+p3

par(mfrow=c(1,2))
 ggplot(data = dati) +
  geom_density(aes(x=Peso), 
             fill = "lightblue", 
             alpha = 0.5, 
             color = "black") +
  geom_vline(xintercept=quantile(dati$Peso,0.25),col="red")+
  geom_vline(xintercept=quantile(dati$Peso,0.75),col="red")+
  labs(x = "Peso (g)",
       y = "Densità") +
  theme_classic()

 ggplot(data = dati) +
  geom_density(aes(x=Lunghezza), 
             fill = "lightblue", 
             alpha = 0.5, 
             color = "black") +
  geom_vline(xintercept=quantile(dati$Lunghezza,0.25),col="red")+
  geom_vline(xintercept =quantile(dati$Lunghezza,0.75),col="red")+
  labs(x = "Lunghezza (mm)",
       y = "Densità") +
  theme_classic()

Si può ora passare all’analisi delle variabili qualitative, calcolando la moda e l’indice di Gini per ciascuna.

moda<-function(x){
  return(names(which.max(table(x))))
}
gini.index<-function(x){
  ni=table(x)
  fi=ni/length(x)
  fi2=fi^2
  j=length(table(x))
  gini=1-sum(fi2)
  gini.normalizzato=round(gini/((j-1)/j),2)
  return(gini.normalizzato)
}
varsq <- dati %>%
  select(Fumatrici, Tipo.parto, Ospedale, Sesso)
tabella3<-data.frame()
for (var in names(varsq)){
  riga<-data.frame(
  Variabile=var,
  Moda = moda(varsq[[var]]),
  Indice_di_Gini = round(gini.index(varsq[[var]]),2)
  )
 tabella3<-rbind(tabella3,riga)
}

kable(tabella3, format = "html", caption = "Statistiche descrittive variabili qualitative") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche descrittive variabili qualitative
Variabile Moda Indice_di_Gini
Fumatrici 0 0.16
Tipo.parto Nat 0.83
Ospedale osp2 1.00
Sesso F 1.00

Osservazioni:

Osservando l’indice di Gini per la variabile Fumatrici, esso tende a 0, quindi essendo la moda 0, significa che è quasi omogenea per la modalità non fumatrici. Le variabili Ospedale e Sesso hanno la massima eterogeneità, in quanto l’indice di Gini è 1, saranno leggermente maggiori i neonati nati nel secondo ospedale e le nascite di femmine. In generale il tipo di parto naturale è stato quello più praticato.

Si possono visualizzare graficamente le distribuzioni relative attraverso dei grafici a torta.

par(mfrow=c(2,2))
for (var in c("Fumatrici","Tipo.parto","Ospedale","Sesso")){
  #cat("Frequenze assolute di", var, ":\n")
  fi<-table(dati[[var]])/2500
  etichette<-paste(rownames(fi),round(fi*100,1),"%")
  pie(fi, main=paste("Distribuzione di",var), radius=1,
      labels=etichette)
}

Si decide a questo punto di suddividere la variabile Anni.madre in classi di tre anni (a meno della prima e dell’ultima) per poi vedere le relazioni di queste con le altre variabili.

dati$Anni_cl<-cut(dati$Anni.madre, breaks=c(12,16,19,21,24,27,30,33,36,39,42,46))
ni<-table(dati$Anni_cl) #frequenza assoluta
fi<-round(table(dati$Anni_cl)/2500,2) #frequenza relativa
Ni<-cumsum(ni) #frequenza cumulata
Fi<-round(Ni/N,2) #frequenza cumulata relativa

distr_freq_Anni_cl<-as.data.frame(cbind(ni,fi,Ni,Fi))
kable(distr_freq_Anni_cl, format = "html", caption = "Tabella delle frequenze") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Tabella delle frequenze
ni fi Ni Fi
(12,16] 22 0.01 22 0.01
(16,19] 87 0.03 109 0.04
(19,21] 140 0.06 249 0.10
(21,24] 346 0.14 595 0.24
(24,27] 561 0.22 1156 0.46
(27,30] 548 0.22 1704 0.68
(30,33] 416 0.17 2120 0.85
(33,36] 226 0.09 2346 0.94
(36,39] 106 0.04 2452 0.98
(39,42] 40 0.02 2492 1.00
(42,46] 8 0.00 2500 1.00

Si sono create 11 classi. L’obiettivo del progetto si pone di vedere se a seconda degli anni e delle variabili Peso e Gestazione, ci siano delle influenze sul peso del neonato. Si può intanto vedere graficamente come si distribuiscono le classi degli anni, attraverso un bar plot.

ggplot(data=dati)+
  geom_bar(aes(x=Anni_cl),
           stat='count',
           col='black',
           fill='red')+
  labs(title="Distribuzione degli anni delle madri",
       x="Classi di Anni",
       y="Frequenze assolute")+
  theme_classic()+
  theme(legend.position="bottom")

Da questo grafico è possibile osservare che in generale la maggior parte delle gravidanze avviene tra i 24 ed i 30 anni. Quindi si nota un picco nelle classi centrali, ed una discesa per età infantili e più anziane. Il peso invece, è influenzato dall’età della madre? Si può provare a vedere come di comporta graficamente, attraverso un line plot, calcolandone la media per classe d’età.

summary_weight_age <- dati %>%
  group_by(Anni_cl) %>%
  summarise(
    mean_weight =round(mean(Peso, na.rm = TRUE),2),
    sd_weight = round(sd(Peso, na.rm = TRUE),2))

ggplot(data=summary_weight_age, aes(x=Anni_cl,y=mean_weight, group=1))+
  geom_line()+
  geom_point()+
  labs(title='Distribuzione delle medie di Peso secondo le varie classi',
       x='Classi Età',
       y='Media Peso')+
  geom_hline(yintercept = 3284.08 , col="red")

  theme(legend.position="bottom")
## List of 1
##  $ legend.position: chr "bottom"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

Osservazione Da questo grafico si può notare come nelle fascia d’età dai 19 ai 33 il peso sia uguale o superiore alla media di 3284.08 dell’intero dataset, mentre diminuisce negli anni precedenti e successivi.

Si può ora indagare come cambia la situazione a seconda della variabile Gestazione. Considerando le settimane di gestazione che vanno da 35 a 42 e ammettendo di definire un parto prematuro avvenuto dalla 35 alla 37 esima settimana, si può guardare come varia il peso medio in questo range ed includere il fatto che le donne siano o meno fumatrici.

summary_weight_gest <- dati %>%
  group_by(Gestazione,Fumatrici) %>%
  summarise(
    mean_weight =round(mean(Peso, na.rm = TRUE),2),
    sd_weight = round(sd(Peso, na.rm = TRUE),2),
    .groups = "drop")

ggplot(data=summary_weight_gest)+
  geom_col(aes(x=Gestazione,    y=mean_weight,fill=Fumatrici),col="black",position=position_dodge(width=0.5) )+
  geom_hline(yintercept = 3284, col="red")+
  labs(title='Distribuzione delle medie dei dei pesi per settimana di gestazione',
       x='Settimane di gestazione',
       y='Media peso')+
  scale_x_continuous(limits = c(34, 43),breaks=34:43)+

  theme(legend.position="bottom")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_col()`).

Osservazioni: Da questo grafico si può vedere, come in generale per tutte le settimane di gestazione la media del peso nel caso di donne non fumatrici sia superiore rispetto alle altre. In questo range le donne che hanno terminato la gravidanza dopo 35 settimane erano non fumatrici, quindi la nascita prematura potrebbe non essere definita dal fumo.

Analisi e Modellizzazione

Analisi Preliminare

Test Statistici per verificare se in alcuni ospedali si fanno più parti cesarei

Per verificarlo avendo due variabili qualitative nominali può essere utilizzato il test Chi-quadrato sulla seguente tabella di contingenza, che calcola le frequenze assolute dei tipi di parto, suddivisi per ospedale.

Imponendo come ipotesi nulla (Ho), che le due variabili siano indipendenti e come alternativa (H1), nel caso venga ammessa un’associazione.

tabella=table(dati$Tipo.parto, dati$Ospedale)
kable(tabella, format = "html", caption = "Tabella di contingenza") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Tabella di contingenza
osp1 osp2 osp3
Ces 242 254 232
Nat 574 595 603 
test.indipendenza<-chisq.test(tabella)
tab <- data.frame(
  Chi_quadrato = round(test.indipendenza$statistic,2),
  P_value = round(test.indipendenza$p.value,2),
  Gradi_libertà = test.indipendenza$parameter
)
kable(tab, format = "html", caption = "Risultati del test Chi-quadro") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati del test Chi-quadro
Chi_quadrato P_value Gradi_libertà
X-squared 1.1 0.58 2

Il valore di chi quadrato approssimato ad una cifra è di 1.1, mentre i gradi di libertà sono 2, ossia 3 numero ospedali - 1.

Considerando un valore di significatività alpha=0.05 non è possibile rifiutare l’ipotesi di indipendeza in quanto il valore di p value è maggiore. Per cui è possibile affermare che non vi sia una differenza significativa nel numero di parti cesarei nei tre ospedali. Graficamente è possibile vederlo anche nel ballon plot sottostante, dove si può notare che i pallini relativi al parto cesareo e naturale rispettivamente, non variano notevolmente le loro dimensioni a seconda dell’ospedale considerato.

ggpubr :: ggballoonplot(data=as.data.frame(tabella), fill="blue")

Si può anche tracciare il grafico e vedere graficamente dove viene collocato il valore di Chi quadrato rispetto alla soglia, calcolata con R per un alpha di 0.05

plot(density(rchisq(1000000,2)), xlim=c(0,15),
     main="Distribuzione Chi-quadro con 2 gradi di libertà",
     xlab= "Valore chi-quadro",
     ylab="Densità")
abline(v=qchisq(0.99,2), col=2)
points(test.indipendenza$statistic,0,cex=3,col=4,pch=20)

Il valore di Chi quadrato è inferiore a quello della soglia, quindi come si era già affermato nel test effettuato con il pvalue, non è possibile rifiutare l’ipotesi nulla.

Test statistico per verificare che la media del peso e della lunghezza del seguente campione di neonati sono significativamente uguali a quelle del resto della popolazione

Sono stati considerati il peso e la lunghezza medi relativi alla popolazione riportati in https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012/#:~:text=In%20media%20il%20peso%20nascita,pari%20mediamente%20a%2050%20centimetri.

In particolare si ha un peso medio di 3300 grammi ed un’altezza di 500 millimetri. Per confrontare le due medie si potrebbe effettuare un t.test, a tal fine è necessario controllare la normalità di ciascuna variabile considerata. Per farlo viene utilizzato sia un metodo analitico, il test di Shapiro-Wilk, che ha come ipotesi nulla (Ho), la normalità, sia una visualizzazione grafica, attraverso il Q-Q plot.

p<-shapiro.test(dati$Peso)
firstrow <- data.frame(
  Variabile="Peso",
  Statistica_W = round(p$statistic,3),
  P_value = round(p$p.value,3)
)
l<-shapiro.test(dati$Lunghezza)
secondrow <- data.frame(
  Variabile="Lunghezza",
  Statistica_W = round(l$statistic,3),
  P_value = round(l$p.value,3)
)
tab<-rbind(firstrow,secondrow)
kable(tab, format = "html", caption = "Risultati Shapiro_Test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati Shapiro_Test
Variabile Statistica_W P_value
W Peso 0.971 0
W1 Lunghezza 0.909 0

La statistica W può assumere valori tra 0 ed 1, scostandosi da 1, ci si discosta dalla distribuzione normale, e queso accade per entrambe le variabili. In particolare il valore del p-value è molto piccolo per entrambi i casi quindi si può con una certa significatività rifiutare l’ipotesi di normalità.

Lo stesso andamento è possibile dimostrarlo attraverso i QQ-plot, che rappresentano graficamente i quantili di ciascuna cistribuzione.

p1<-ggplot(data=dati, aes(sample = Peso)) +
  stat_qq() +
  stat_qq_line(col="red")+
  labs(title="Peso",
    x = "Quantili Teorici",
       y = "Quantili del Campione") +
  theme_classic()
p2<-ggplot(data=dati, aes(sample = Lunghezza)) +
  stat_qq() +
  stat_qq_line(col="red")+
  labs(title="Lunghezza",
    x = "Quantili Teorici",
       y = "Quantili del Campione") +
  theme_classic()
p1+p2

Se la variabile avesse distribuzione normale, i punti individuati nel seguente grafico, dovrebbero addensarsi attorno alla diagonale rossa, come si può vedere in entrambi i casi vi è una deviazione sia all’inizio che alla fine.

Nonostante non si abbia la completa normalità, ma avendo un campione sufficientemente grande, si può provare comunque ad effetuare il test t a due code, imponendo come ipotesi nulla (Ho), l’uguaglianza delle medie, e come ipotesi alternativa, la loro differenza.

weight<-t.test(dati$Peso,mu=3300, conf.level=0.95, alternative="two.sided")
primariga <- data.frame(
  Variabile="Peso",
  Statistica_t = round(weight$statistic,2),
  P_value = round(weight$p.value,2),
  Media_Campione = round(mean(dati$Peso),2),
  Media_Popolazione = 3300,
  Intervallo_Confidenza = paste(round(weight$conf.int[1], 2), "-", round(weight$conf.int[2], 2)),
  Gradi_libertà = weight$parameter
)
length<-t.test(dati$Lunghezza,mu=500, conf.level=0.95, alternative="two.sided")
secondariga <- data.frame(
  Variabile="Lunghezza",
  Statistica_t = round(length$statistic,2),
  P_value = round(length$p.value,2),
  Media_Campione = round(mean(dati$Lunghezza),2),
  Media_Popolazione = 500,
  Intervallo_Confidenza = paste(round(length$conf.int[1], 2), "-", round(length$conf.int[2], 2)),
  Gradi_libertà = length$parameter
)
tabella<-rbind(primariga,secondariga)
kable(tabella, format = "html", caption = "Risultati t-test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati t-test
Variabile Statistica_t P_value Media_Campione Media_Popolazione Intervallo_Confidenza Gradi_libertà
t Peso -1.52 0.13 3284.08 3300 3263.49 - 3304.67 2499
t1 Lunghezza -10.08 0.00 494.69 500 493.66 - 495.72 2499 
pone<-plot(density(rt(100000,2499),xlim=c(-4,4)),
     main="Distribuzione T del Peso",
     xlab= "Valore t",
     ylab="Densità")
## Warning: In density.default(rt(1e+05, 2499), xlim = c(-4, 4)) :
##  l'argomento extra 'xlim' sarà ignorato
abline(v=qt(0.025,2499), col=2)
abline(v=qt(0.975,2499), col=2)
points(weight$statistic,0,cex=3,col=4,pch=20)

ptwo<-plot(density(rt(10000000,2499),xlim=c(-11,11)),
     main="Distribuzione T della lunghezza",
     xlab= "Valore t",
     ylab="Densità")
## Warning: In density.default(rt(1e+07, 2499), xlim = c(-11, 11)) :
##  l'argomento extra 'xlim' sarà ignorato
abline(v=qt(0.025,2499), col=2)
abline(v=qt(0.975,2499), col=2)
points(length$statistic,0,cex=3,col=4,pch=20)

pone+ptwo
## integer(0)

Osservazione:

Dai risultati dei test t ed i relativi grafici effettuati per entrambe le variabili è possibile osservare che nel caso della lunghezza è possibile rifiutare con una certa significatività l’ipotesi di uguaglianza tra le medie ed ammettere l’esistenza di una differenza, in particolare la statistica è inferiore alla soglia fissata per un alpha di 0.05. Anche l’intervallo di confidenza esclude la media della popolazione. E come si può vedere dal grafico la statistica uscirebbe essendo -10 dalla soglia.

Nel caso invece del Peso non è possibile fare le stesse osservazioni, in quanto in questo caso non è possibile ammettere una differenza sufficientemente grande da rifiutare l’ipotesi nulla. Anche l’intervallo di confidenza in questo caso contiene la media della popolazione.

Vista la non normalità di entrambe le variabili, può esere effettuato anche un secondo test non parametrico di controllo, il Wilcoxon-Mann-Whitney, che si basa sui ranghi, a differenza del test t utilizza come indice di tendenza centrale la mediana. Può dar meno peso a differenze piccole tra i due gruppi.

w_peso<-wilcox.test(dati$Peso,mu=3300, conf.level=0.95, alternative="two.sided")
uno <- data.frame(
  Variabile="Peso",
  Statistica_W = round(w_peso$statistic,2),
  P_value = round(w_peso$p.value,2)
)
w_lun<-wilcox.test(dati$Lunghezza,mu=500, conf.level=0.95, alternative="two.sided")
due <- data.frame(
  Variabile="Lunghezza",
  Statistica_W = round(w_lun$statistic,2),
  P_value = round(w_lun$p.value,2)
)
tabella<-rbind(uno,due)
kable(tabella, format = "html", caption = "Risultati Wilcox-test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati Wilcox-test
Variabile Statistica_W P_value
V Peso 1495594 0.96
V1 Lunghezza 877236 0.00

Tale test conferma quanto anticipato con il T-test, ossia fissando una soglia alpha dello 0.05, si può concludere che nel caso della Variabile lunghezza esista una differenza significativa tra la media del campione e della popolazione, mentre questo effetto non è riscontrabile nel caso del Peso.

Test per vedere se le misure antropometriche siano significativamente diverse tra i due sessi (Peso, Lunghezza, Cranio)

Sarebbe utile osservare i box plot per ciascuna variabile, suddivisa secondo il sesso, per avere una prima idea visiva di eventuali differenze:

p<-ggplot(data=dati)+
 geom_boxplot(aes(x=Sesso,y=Peso, fill=Sesso))+
  labs(x="Sesso",
       y="Peso alla nascita (g)")+
  theme_classic()+
  theme(legend.position="bottom")
l<-ggplot(data=dati)+
 geom_boxplot(aes(x=Sesso,y=Lunghezza, fill=Sesso))+
  labs(x="Sesso",
       y="Lunghezza alla nascita (mm)")+
  theme_classic()+
  theme(legend.position="bottom")
c<-ggplot(data=dati)+
 geom_boxplot(aes(x=Sesso,y=Cranio, fill=Sesso))+
  labs(x="Sesso",
       y="Lunghezza del cranio (mm)")+
  theme_classic()+
  theme(legend.position="bottom")
p+l+c

In questi tre grafici sembrerebbe che la differenza più significativa sia presente nel caso del peso. Si può affermare che in generale in tutte e tre le variabili l’IQR è circa lo stesso per entrambi i sessi. E in generale la mediana di ciascuna variabile è superiore nel caso del sesso maschile.

In questo caso è possibile utilizzare un test t ed un test di Wilcox per confronti multipli.

wp<-pairwise.wilcox.test(dati$Peso,dati$Sesso, paired=F,pool.sd=T, p.adjust="bonferroni")
wl<-pairwise.wilcox.test(dati$Lunghezza,dati$Sesso, paired=F,pool.sd=T, p.adjust="bonferroni")
wc<-pairwise.wilcox.test(dati$Cranio,dati$Sesso, paired=F,pool.sd=T, p.adjust="bonferroni")

tab <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  P_value = round(c(wp$p.value, wl$p.value, wc$p.value), 4)
)

kable(tab, format = "html", caption = "Risultati Wilcox-test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati Wilcox-test
Variabile P_value
Peso 0
Lunghezza 0
Cranio 0

Avendo i p_value per tutte e tre le variabili inferiori a zero si può concludere che vi sia una differenza significativa in tutte le variabili a seconda del sesso.

Da ultimo è possibile effetuare anche un t-test per tutte e tre le variabili.

attach(dati)
weight<-t.test(Peso~Sesso, conf.level=0.95, alternative="two.sided")
primariga <- data.frame(
  Variabile="Peso",
  Statistica_t = round(weight$statistic,2),
  P_value = round(weight$p.value,2),
  
  Intervallo_Confidenza = paste(round(weight$conf.int[1], 2), "-", round(weight$conf.int[2], 2)),
  Gradi_libertà = weight$parameter
)
length<-t.test(Lunghezza~Sesso, conf.level=0.95, alternative="two.sided")
secondariga <- data.frame(
  Variabile="Lunghezza",
  Statistica_t = round(length$statistic,2),
  P_value = round(length$p.value,2),
  
  Intervallo_Confidenza = paste(round(length$conf.int[1], 2), "-", round(length$conf.int[2], 2)),
  Gradi_libertà = length$parameter
)
cran<-t.test(Cranio~Sesso, conf.level=0.95, alternative="two.sided")
terzariga <- data.frame(
  Variabile="Cranio",
  Statistica_t = round(cran$statistic,2),
  P_value = round(cran$p.value,2),
 
  Intervallo_Confidenza = paste(round(length$conf.int[1], 2), "-", round(length$conf.int[2], 2)),
  Gradi_libertà = cran$parameter
)
tabella<-rbind(primariga,secondariga,terzariga)
kable(tabella, format = "html", caption = "Risultati t-test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati t-test
Variabile Statistica_t P_value Intervallo_Confidenza Gradi_libertà
t Peso -12.11 0 -287.11 - -207.06 2490.716
t1 Lunghezza -9.58 0 -11.93 - -7.88 2459.302
t2 Cranio -7.41 0 -11.93 - -7.88 2491.389

Anche in questo caso, osservando la colonna relativa al p_value, per tutte e tre le variabili è possibile rifiutare l’ipotesi nulla di uguaglianza, ammettendo una differenza tra le medie a seconda del sesso per ciascuna variabile.

Creazione del Modello di Regressione

Si vuole sviluppare un modello di regressione lineare multipla che includa tutte le variabili, in modo da quantificare l’impatto di ciascuna di queste sul peso del neonato ed eventuali interazioni.

Che distribuzione ha la variabile Y (il peso del neonato), che si vuole prevedere? Per capirlo vengono effettuati il test di Fisher, per vedere eventuali simmetrie e calcolata la curtosi, per avere un’idea della forma della distribuzione rispetto alla normale.

moments :: skewness(dati$Peso)
## [1] -0.6470308
moments :: kurtosis(dati$Peso)-3
## [1] 2.031532

L’indice di fisher risulta negativo, quindi il peso presenterà una distribuzione asimmetrica negativa, per cui la media è minore della mediana. Invece la curtosi è positiva, quindi sarà più allungata verso l’alto rispetto alla normale. Si può ossservare tale comportamento anche plottando la distribuzione:

ggplot()+
  geom_density(aes(x=dati$Peso),col="black",fill="lightblue")+
  geom_vline(xintercept = mean(dati$Peso), color = "red", linetype = "dashed", linewidth = 1)+
  geom_vline(xintercept = median(dati$Peso), color = "black", linetype = "dashed", linewidth = 1)+
  labs(x="Peso",
       y="Densità")

Al fine di costruire il modello di regressione multipla è importante capire come interagiscono le variabili tra di loro e con la variabile Peso. Per questo viene costruito il grafico sottostante, che nella matrice triangolare superiore presnta gli scatter plot tra le variabili e nella parte inferiore restituisce i coefficienti di correlazione. Vengono considerate le variabili quantitative.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(vars,lower.panel = panel.cor,upper.panel = panel.smooth)

Osservazioni:

Vengono riportate le correlazioni seguenti:

  • Peso e anni madre= -0.022 debole negativa, infatti lo scatter plot si presenta orizzontale e sparso;
  • Peso e N.gravidanze= 0.0024 debole positiva;
  • Peso e Gestazione= 0.59 positiva, sembra abbastanza logaritmica;
  • Peso e Lunghezza= 0.8 correlazione fortemente positiva, non propriamente lineare;
  • Peso e Cranio = 0.7 positiva, anche in questo caso abbastanza lineare;

In generale le relazioni con correlazione maggiore con il Peso, sono quelle relative alle variabili Gestazione, Lunghezza e Cranio, che tenderanno ad influenzare di più il modello.

Si può invece osservare come si comportano le variabili qualitative rispetto al peso attraverso dei boxplot, riportiamo ancora quello relativo al sesso, anche se già visto in precedenza.

par(mfrow=c(2,2))
ggplot(data=dati)+
 geom_boxplot(aes(x=Sesso,y=Peso, fill=Sesso))+
  labs(x="Sesso",
       y="Peso alla nascita (g)")+
  theme_classic()+
  theme(legend.position="bottom")

ggplot(data=dati)+
 geom_boxplot(aes(x=factor(Fumatrici),y=Peso, fill=Fumatrici))+
  labs(x="Fumatrici",
       y="Peso alla nascita (g)")+
  theme_classic()+
  theme(legend.position="bottom")

ggplot(data=dati)+
 geom_boxplot(aes(x=Tipo.parto,y=Peso, fill=Tipo.parto))+
  labs(x="Tipo di parto",
       y="Peso alla nascita (g)")+
  theme_classic()+
  theme(legend.position="bottom")

ggplot(data=dati)+
geom_boxplot(aes(x=Ospedale,y=Peso, fill=Ospedale))+
  labs(x="Ospedale",
       y="Peso alla nascita (g)")+
  theme_classic()+
  theme(legend.position="bottom")

Osservazioni:

  • Sesso: come si era già accennato in precedenza, il peso, varia a seconda del sesso e lo si può vedere dalla mediana che nei maschi è posta più alta rispetto alle femmine. L’IQR invece sembra essere abbastanza simile, si può notare che per entrambi i sessi sono presenti alcuni outlier sia inferiormente che superiormente;
  • Fumatrici: il peso è leggermente più elevato nel caso di bambini nati da donne non fumatrici, in questo caso si ha più variabilità, bisogna anche ricordarsi però la presenza di più campioni.
  • Tipo.parto: in questo caso i pesi sono molto simili, la differenza più evidente sta nel numero di outlier presenti, che nel caso di parti naturali aumenta.
  • Ospedale: per questa variabile non si hanno troppe differenze, a parte un leggero calo della mediana per il terzo ospedale. In generale in tutti e tre casi gli outlier sono più presenti nella parte inferiore.

Si può anche usare un t.test, in particolare viene esclusa la variabile Ospedale in quanto si osservavano le differenze minori,per vedere se le differenze tra le medie siano significative.

weights<-t.test(Peso~Sesso)
primariga <- data.frame(
  Variabile="Sesso",
  Statistica_t = round(weights$statistic,2),
  P_value = round(weights$p.value,2),
  Media1 = round(weights$estimate[1],2),
  Media2 = round(weights$estimate[2],2)
)
weightf<-t.test(Peso~Fumatrici)

secondariga <- data.frame(
  Variabile="Fumatrici",
  Statistica_t = round(weightf$statistic,2),
  P_value = round(weightf$p.value,2),
  Media1=round(weightf$estimate[1],2),
  Media2 = round(weightf$estimate[2],2)
  
)
weightt<-t.test(Peso~Tipo.parto)
terzariga <- data.frame(
  Variabile="Tipo.parto",
  Statistica_t = round(weightt$statistic,2),
  P_value = round(weightt$p.value,2),
  Media1= round(weightt$estimate[1],2),
  Media2 = round(weightt$estimate[2],2)
  
)
tabella<-rbind(primariga,secondariga,terzariga)
kable(tabella, format = "html", caption = "Risultati t-test") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati t-test
Variabile Statistica_t P_value Media1 Media2
t Sesso -12.11 0.0 3161.13 3408.22
t1 Fumatrici 1.03 0.3 3286.15 3236.35
t2 Tipo.parto -0.13 0.9 3282.05 3284.92

Dal t-test è possibile rifiutare l’ipotesi nulla di uguaglianza solo per la varibaile sesso, mentre per le altre il valore di p_value supera la soglia di 0.05 quindi non è possibile affermare che esista una differenza significativa.

Selezione del modello ottimale

A questo punto, si può creare il primo modello, considerando tutte le variabili, tranne quella relativa alla classe degli anni.

mod1<-lm(Peso~.-Anni_cl,data=dati)
kable(summary(mod1)$coefficients, format = "html", caption = "Statistiche primo modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche primo modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6735.1677173 141.3977471 -47.6327795 0.0000000
Anni.madre 0.7983057 1.1463259 0.6964038 0.4862410
N.gravidanze 11.4117638 4.6664758 2.4454780 0.0145349
Fumatrici -30.1566601 27.5395750 -1.0950300 0.2736095
Gestazione 32.5265096 3.8178659 8.5195527 0.0000000
Lunghezza 10.2950617 0.3006984 34.2371680 0.0000000
Cranio 10.4724852 0.4260520 24.5802992 0.0000000
Tipo.partoNat 29.5026702 12.0847811 2.4413078 0.0147034
Ospedaleosp2 -11.2216217 13.4387547 -0.8350195 0.4037869
Ospedaleosp3 28.0984131 13.4971965 2.0817963 0.0374630
SessoM 77.5472538 11.1779356 6.9375292 0.0000000 
tabmod1 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod1)$r.squared, summary(mod1)$adj.r.squared, summary(mod1)$fstatistic["value"]), 4)
)

kable(tabmod1, format = "html", caption = "Statistiche Mod1") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod1
Variabile Valore
Rquadro 0.7289
Rquadro_adj 0.7278
F-statistic 669.1376

Il valore di Rquadro è di 0.7289, mentre quello aggiustato di 0.7278, si potrebbe ancora migliorare, in quanto al momento il modello spiega il 72.78 % della variabilità della Risposta Peso. La statistica è 669.13, con un p-value vicino allo 0, quindi il modello risulta essere significativo. In generale si può osservare che solo Anni.madre, Fumatrici, Ospedaleosp2 risultano avere dei pvalue elevati e superiori alla soglia di 0.05

Per migliorare ed ottimizzare il modello, si può provare ad eliminare la variabile Anni.madre, visto il basso coefficiente e l’alto p-value, ci si aspetta che influisca poco nel modello.

mod2<-update(mod1,~.-Anni.madre)
kable(summary(mod2)$coefficients, format = "html", caption = "Statistiche secondo modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche secondo modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6708.10650 135.939383 -49.3463068 0.0000000
N.gravidanze 12.60846 4.338111 2.9064414 0.0036879
Fumatrici -30.30924 27.535855 -1.1007191 0.2711253
Gestazione 32.25015 3.796793 8.4940502 0.0000000
Lunghezza 10.29445 0.300666 34.2388134 0.0000000
Cranio 10.48764 0.425452 24.6505749 0.0000000
Tipo.partoNat 29.53506 12.083442 2.4442587 0.0145839
Ospedaleosp2 -11.08165 13.435862 -0.8247812 0.4095748
Ospedaleosp3 28.36596 13.490332 2.1026882 0.0355932
SessoM 77.62052 11.176284 6.9451098 0.0000000 
tabmod2 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod2)$r.squared, summary(mod2)$adj.r.squared, summary(mod2)$fstatistic["value"]), 4)
)

kable(tabmod2, format = "html", caption = "Statistiche Mod2") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod2
Variabile Valore
Rquadro 0.7288
Rquadro_adj 0.7278
F-statistic 743.5861

Si può notare che non sono stati prodotti cambiamenti significativi, il valore di Rquadro si è mantenuto, e si è alzata la F-statistic.

Si può utilizzare il test di ANOVA per guardare se il cambiamento della varianza spiegata sia stato significativo o meno.

c<-anova(mod2,mod1)
round(c[2, "Pr(>F)"],2)
## [1] 0.49

Pvalue 0.49, quindi l’aumento della varianza spiegata non è significativo, e si può decidere di utilizzare il secondo modello, escludendo la variabile Anni.madre. Si può poi andare ad utilizzare un altro criterio di valutazione: BIC. In generale penalizza modelli con più varibaili, incoraggiando la scelta di modelli più semplici, che spieeghino adeguatamente i dati senza inutili complessità.

kable(BIC(mod2,mod1), format = "html", caption = "BIC") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
BIC
df BIC
mod2 11 35234.64
mod1 12 35241.97

Come si può vedere il secondo modello che utilizza 11 variabili presenta un valore di BIC minore, per questo è preferibile.

Possono essere calcolati anche gli indicatori di multicollinearità (VIF). La multicollinearità si verifica quando due o più variabili indipendenti sono altamente correlate, attraverso questi coefficienti è possibile valutare quanto aumenta lla varianza di un coefficiente di rigressione stimato quando i predittori sono correlati.

vif_values <- vif(mod2)
vif_table <- data.frame(
  Variabile = rownames(vif_values),
  VIF = round(vif_values[,1], 2)
)
kable(vif_table[-1], format = "html", caption = "VIF") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
VIF
VIF
N.gravidanze 1.03
Fumatrici 1.01
Gestazione 1.68
Lunghezza 2.09
Cranio 1.63
Tipo.parto 1.00
Ospedale 1.00
Sesso 1.04

Sarebbe preferibile trovarli vicino a 0, in ogni caso è preferibile che siano tutti sotto al valore 5. In particolare in questo caso è indicata una correlazione moderata perché hanno valori compresi tra 1 e 5.

Si può ora provare a costruire un terzo modello rimuovendo la variabile Ospedale.

mod3<-update(mod2,~.-Ospedale)
kable(round(summary(mod3)$coefficients,3), format = "html", caption = "Coefficienti terzo modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti terzo modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6708.074 135.984 -49.330 0.000
N.gravidanze 13.012 4.342 2.997 0.003
Fumatrici -31.759 27.571 -1.152 0.249
Gestazione 32.541 3.801 8.561 0.000
Lunghezza 10.272 0.301 34.129 0.000
Cranio 10.501 0.426 24.648 0.000
Tipo.partoNat 30.296 12.098 2.504 0.012
SessoM 78.114 11.191 6.980 0.000 
tabmod3 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod3)$r.squared, summary(mod3)$adj.r.squared, summary(mod3)$fstatistic["value"]), 4)
)

kable(tabmod3, format = "html", caption = "Statistiche Mod3") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod3
Variabile Valore
Rquadro 0.7278
Rquadro_adj 0.7271
F-statistic 951.9569

Con la costruzione del terzo modello, si può osservare che l’unica variabile non significativa è Fumatrici e si abbassa di pochissimo il valore di Rquadro aggiustato 0.7271. Si può dunque eseguire un’ulteriore modifica a quest’ultimo modello, con la rimozione della variabile Fumatrici

mod4<-update(mod3,~.-Fumatrici)
kable(round(summary(mod4)$coefficients,3), format = "html", caption = "Coefficienti quarto modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti quarto modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6707.297 135.991 -49.322 0.000
N.gravidanze 12.756 4.337 2.941 0.003
Gestazione 32.271 3.794 8.506 0.000
Lunghezza 10.286 0.301 34.207 0.000
Cranio 10.506 0.426 24.659 0.000
Tipo.partoNat 30.034 12.097 2.483 0.013
SessoM 77.929 11.191 6.964 0.000 
tabmod4 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod4)$r.squared, summary(mod4)$adj.r.squared, summary(mod4)$fstatistic["value"]), 4)
)

kable(tabmod4, format = "html", caption = "Statistiche Mod4") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod4
Variabile Valore
Rquadro 0.7277
Rquadro_adj 0.7270
F-statistic 1110.2496

Si può notare che i valori dei pvalue sono tutti sotto alla soglia di 0.05, quelli più elevati sono relativi al numero di gravidanze e al tipo di parto, ma comunque si collocano inferiormanete.

Il parametro R quadro aggiustato rimane stabile, a 0.727, non è una performance malvagia, ma si potrebbe fare di più. Confrontiamo il valore BIC tra i vari modelli aggiunti.

kable(BIC(mod2,mod3,mod4), format = "html", caption = "BIC mod2/mod3/mod4") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
BIC mod2/mod3/mod4
df BIC
mod2 11 35234.64
mod3 9 35228.24
mod4 8 35221.75

Il BIC con il valore più basso risulta essere quello relativo al modello 4, costruito con 8 variabili. Quindi le prossime modifiche partiranno da quest’ultimo modello.

In particolare dall’analisi delle correlazioni, si poteva vedere che la variabile Gestazione era molto correlata con Cranio e Lunghezza, con valori di correlazione superiori a 0.4 e la forma tendente al logaritmo.

Per questo motivo si può provare a fare interagire tali variabili nel modello ed osservare se si ottengono miglioramenti.

mod5<-update(mod4,~.+Gestazione*Cranio+Gestazione*Lunghezza)
kable(round(summary(mod5)$coefficients,3), format = "html", caption = "Coefficienti quinto modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti quinto modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -274.049 1106.260 -0.248 0.804
N.gravidanze 13.416 4.310 3.113 0.002
Gestazione -139.403 29.508 -4.724 0.000
Lunghezza 9.149 3.754 2.437 0.015
Cranio -7.762 6.429 -1.207 0.227
Tipo.partoNat 28.766 12.021 2.393 0.017
SessoM 71.863 11.174 6.431 0.000
Gestazione:Cranio 0.479 0.166 2.882 0.004
Gestazione:Lunghezza 0.035 0.097 0.362 0.717 
tabmod5 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod5)$r.squared, summary(mod5)$adj.r.squared, summary(mod5)$fstatistic["value"]), 4)
)

kable(tabmod5, format = "html", caption = "Statistiche Mod5") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod5
Variabile Valore
Rquadro 0.7314
Rquadro_adj 0.7305
F-statistic 847.8896

Si osserva che Gestazione:Lunghezza non è stata un’aggiunta significativa al modello, in quanto il coefficiente è basso ed il p value molto alto. Il coefficiente relativo alla variabile Cranio ha perso di significatività. Il valore di R quadro aggiustato si è però alzato di un centesimo.

Si può provare a perfezionare il modello con la rimozione dell’interazione tra Gestazione e Lunghezza.

mod6<-update(mod5,~.-Gestazione:Lunghezza)
kable(round(summary(mod6)$coefficients,3), format = "html", caption = "Coefficienti sesto modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti sesto modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) -265.029 1105.787 -0.240 0.811
N.gravidanze 13.406 4.309 3.111 0.002
Gestazione -139.487 29.502 -4.728 0.000
Lunghezza 10.505 0.301 34.898 0.000
Cranio -9.723 3.472 -2.800 0.005
Tipo.partoNat 28.776 12.018 2.394 0.017
SessoM 72.038 11.161 6.454 0.000
Gestazione:Cranio 0.530 0.090 5.870 0.000 
tabmod6 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod6)$r.squared, summary(mod6)$adj.r.squared, summary(mod6)$fstatistic["value"]), 4)
)

kable(tabmod6, format = "html", caption = "Statistiche Mod6") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod6
Variabile Valore
Rquadro 0.7314
Rquadro_adj 0.7306
F-statistic 969.3358

Così facendo, senza perdere il valore di R quadro, la variabile Cranio è tornata ad avere significatività.

Dal grafico delle correlazioni e come rimarcato precedentemente si può osservare una dipendenza logaritmica tra le variabili Gestazione e lunghezza, per questo viene aggiunto tale effetto.

mod7<-update(mod6,~.+(log(Lunghezza))+(log(Gestazione)))
kable(round(summary(mod7)$coefficients,3), format = "html", caption = "Coefficienti settimo modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti settimo modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) 48950.452 8772.718 5.580 0.000
N.gravidanze 14.898 4.234 3.518 0.000
Gestazione -404.947 106.187 -3.814 0.000
Lunghezza 44.966 3.743 12.013 0.000
Cranio -1.428 5.814 -0.246 0.806
Tipo.partoNat 28.041 11.802 2.376 0.018
SessoM 71.244 10.964 6.498 0.000
log(Lunghezza) -16687.607 1803.304 -9.254 0.000
log(Gestazione) 13031.145 2578.921 5.053 0.000
Gestazione:Cranio 0.307 0.151 2.039 0.042 
tabmod7 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod7)$r.squared, summary(mod7)$adj.r.squared, summary(mod7)$fstatistic["value"]), 4)
)

kable(tabmod7, format = "html", caption = "Statistiche Mod7") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod7
Variabile Valore
Rquadro 0.7412
Rquadro_adj 0.7403
F-statistic 792.3806

Si può osservare un miglioramento del valore di Rquadro aggiustato, ed una perdita di significatività della variabile Cranio e dell’interazione Gestazione:Cranio. Si può provare in ultima analisi a rimuovere quest’ultima.

mod8<-update(mod7,~.-Gestazione:Cranio)
kable(round(summary(mod8)$coefficients,3), format = "html", caption = "Coefficienti ottavo modello") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti ottavo modello
Estimate Std. Error t value Pr(>|t|)
(Intercept) 60697.567 6621.173 9.167 0.000
N.gravidanze 14.792 4.237 3.491 0.000
Gestazione -219.468 54.852 -4.001 0.000
Lunghezza 47.901 3.458 13.854 0.000
Cranio 10.398 0.418 24.859 0.000
Tipo.partoNat 28.127 11.810 2.382 0.017
SessoM 71.972 10.965 6.564 0.000
log(Lunghezza) -18110.981 1663.830 -10.885 0.000
log(Gestazione) 9877.929 2065.390 4.783 0.000 
tabmod8 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod8)$r.squared, summary(mod8)$adj.r.squared, summary(mod8)$fstatistic["value"]), 4)
)

kable(tabmod8, format = "html", caption = "Statistiche Mod8") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod8
Variabile Valore
Rquadro 0.7408
Rquadro_adj 0.7399
F-statistic 889.7797

Si abbassa leggermente il valore di R quadro aggiustato, ma ora tutti i parametri sembrano essere significativi per il modello. Controlliamo tali risultati con il calcolo del BIC a partire dal modello 4.

kable(BIC(mod4,mod5,mod6,mod7,mod8), format = "html", caption = "BIC mod4/mod5/mod6/mod7/mod8") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
BIC mod4/mod5/mod6/mod7/mod8
df BIC
mod4 8 35221.75
mod5 10 35202.94
mod6 9 35195.24
mod7 11 35117.84
mod8 10 35114.19

Il caloro del BIC conferma la preferenza per l’ultimo modello costruito con 10 variabili.

Diagnostica dei residui

Può essere a questo punto eseguita una diagnostica approfondita dei residui del modello e dei potenziali valori influenti.

Con la costruzione del grafico successivo si vuole controllate che vengano rispettate le assunzioni del modello, ossia: - La linearità dei residui con i predittori; - La distrubzione normale dei residui; - L’omoschedicità: ossia la dispersione dei residui deve essere costante. Dovrebbero avere varianza costante e quindi formare una nuvola costante attorno l’orizzontale. - Devono poi essere indipendenti ed identicamente distribuiti (i.i.d)

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

Plot 1: Residuals vs. Fitted

In questo grafico i residui dovrebbero posizionarsi casualmente attorno alla media di 0, in questo caso si vede una certa concentrazione nella parte più a destra, come una sorta di asimmetria.

Plot 2: Q-Q Residuals

Viene utilizzato per capire se i residui si distribuiscono in modo normale, come si può vedere in questo caso nella parte finale ed iniziale tendono a staccarsi dalla bisettrice, rappresentata in rosso.

Plot 3: Scale-Location:

In questo caso si dovrebbe visualizzare una nuvola casuale attorno ad un valore di y, per verificare la proprietà di omoschedicità, come nel primo grafico si può notare una certa tendenza di accumulazione verso destra.

Plot 4: Residuals vs Leverage

Si può osservare che l’osservazione 1551 si trova nella soglia di avvertimento dettata dalla distanza di Cook.

Queste prime osservazioni grafiche, possono essere confermate attraverso dei test statistici.

Si può cominciare con il calcolo della normalità con il test di Shapito-WIlk:

st<-shapiro.test(residuals(mod8))$p.value

L’omoschedicità, ossia la varianza costante, può essere verificata con il test di Breusch-Pagan:

bp<-bptest(mod8)$p.value

Con il test di Darbin-Watson si può invece verificare l’indipendenza dei residui.

dw<-dwtest(mod8)$p.value
tab2 <- data.frame(
  Test = c("Shapiro-Wilk", "Breusch-Pagan", "Darbin-Watson"),
  P_value = round(c(st,bp,dw), 4)
)

kable(tab2, format = "html", caption = "Risultati Test sui residui") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Risultati Test sui residui
Test P_value
Shapiro-Wilk 0.0000
Breusch-Pagan 0.0000
Darbin-Watson 0.1018

Quindi tutte le ipotesi effettuate possono essere rifiutate a causa dei bassi valori di pvalue. Si riconfermano quindi i risultati grafici.

In questi grafici è possibile avere un’ulteriore rappresentazione grafica dei residui e dei leverage.

par(mfrow=c(2,2))

plot(density(residuals(mod8)), ylab="Densità", xlab="Residui", main="")
#leverage

lev<-hatvalues(mod8)
plot(lev, ylab="Leverage", xlab="Indici")
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
#lev[lev>soglia]
#sum(lev>soglia)
#outliers
plot(rstudent(mod8), ylab="Residui Studentizzati", xlab="Indici")
abline(h=c(-2,2),col=2)
car::outlierTest(mod8)
##       rstudent unadjusted p-value Bonferroni p
## 1551  5.538098         3.3784e-08   8.4459e-05
## 1306  4.927276         8.8857e-07   2.2214e-03
## 155   4.538064         5.9469e-06   1.4867e-02
## 1694  4.414825         1.0540e-05   2.6351e-02
## 1399 -4.400874         1.1236e-05   2.8090e-02
out<-sum(abs(rstudent(mod8))>2)
#out

#distanza di cook
cook<-cooks.distance(mod8)
plot(cook, ylim = c(0,1), ylab="Distanza di Cook", xlab="Indici")

In particolare il primo grafico rappresenta la distribuzione, che come si può osservare, non è identica alla normale, come già si poteva intuire dal test di Shapiro-Wilk. Si hanno poi 106 punti leverage. Si può altresì provare ad eliminare l’outlier 1551 e vedere come cambia il modello, rimuovendo anche la variabile log(Gestazione).

dati2<-dati[-1551,]
mod9<-lm(Peso~. ,data=dati2)
mod10<-update(mod1,~.-Anni.madre-Fumatrici-Ospedale+log(Lunghezza) )
kable(round(summary(mod10)$coefficients,2), format = "html", caption = "Coefficienti senza outlier") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Coefficienti senza outlier
Estimate Std. Error t value Pr(>|t|)
(Intercept) 60463.16 6649.99 9.09 0.00
N.gravidanze 14.43 4.25 3.39 0.00
Gestazione 42.23 3.85 10.97 0.00
Lunghezza 37.64 2.72 13.82 0.00
Cranio 10.61 0.42 25.40 0.00
Tipo.partoNat 27.76 11.86 2.34 0.02
SessoM 69.59 11.00 6.33 0.00
log(Lunghezza) -13079.21 1294.60 -10.10 0.00 
tabmod10 <- data.frame(
  Variabile = c("Rquadro", "Rquadro_adj", "F-statistic"),
  Valore = round(c(summary(mod10)$r.squared, summary(mod10)$adj.r.squared, summary(mod10)$fstatistic["value"]), 4)
)

kable(tabmod10, format = "html", caption = "Statistiche Mod10") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Statistiche Mod10
Variabile Valore
Rquadro 0.7384
Rquadro_adj 0.7377
F-statistic 1004.8040

Osservazioni: Il modello senza l’outlier non presenta grandi variazioni, ma aumenta leggermente il valore di R quadro aggiustato. Quindi ora spiega il 73,77% della variabilità del Peso.

Quanto il modello è buono per fare previsioni?

In generale alla fine dell’analisi, i modelli presentano tutti predittori significativi, considerando un valore di soglia dello 0.05. Il valore di BIC è il più basso rispetto a quello iniziale, aumentando l’efficacia del modello. Inoltre il valore di R quadro si è alzato e al momento spiega il 73,77% della variabilità del Peso. Questo modello sembra quindi essere robusto e consistente, si può quindi fare un test su un dato inserito da zero. In particolare si vuole prevedere il peso di una neonata, la cui madre è alla terza gravidanza e partorirà alla 39 esima settimana.

Prevedere il peso di una neonata alla 39esima settimana con madre alla terza gravidanza

Non avendo alcun riferimento per le varibaili Lunghezza, Cranio e Tipo.parto, si potrebbero usare i valori mediani delle variabili quantitative e fare due test nel caso di parto naturale o cesareo.

new_line1 <- data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = mean(dati$Lunghezza),
  Cranio = mean(dati$Cranio),
  Tipo.parto = "Ces",
  Sesso = "F"
)
new_line2 <- data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = mean(dati$Lunghezza),
  Cranio = mean(dati$Cranio),
  Tipo.parto = "Nat",
  Sesso = "F"
)
prev.testces<- predict(mod10,newdata = new_line1)
prev.testnat<- predict(mod10,newdata = new_line2)

tabmod11 <- data.frame(
  Variabile = c("Cesareo", "Naturale" ),
  Valore = round(c(prev.testces, prev.testnat) , 2)

)
kable(tabmod11, format = "html", caption = "Previsioni") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Previsioni
Variabile Valore
Cesareo 3239.92
Naturale 3267.67

In generale i risultati sono buoni considerando una media del peso di 3300 grammi e considerando che il neonato nato dopo 39 settimane non è considerato prematuro.

Visualizzazioni Finali

In ultima analisi Viene creato una scatterplot 3d, per non perdere informazioni e visualizzare le differenze tra maschi e femmine nel peso, a seconda del variare delle variabili Lunghezza e Cranio.

colori <- ifelse(dati$Sesso == "M", "blue", "red")
s3d <- scatterplot3d(dati$Lunghezza, dati$Cranio, dati$Peso,
                     pch = 16,
                     color = colori,
                     main = "3D Scatter Plot",
                     xlab = "Lunghezza (mm)",
                     ylab = "Cranio (mm)",
                     zlab = "Peso (g)")
legend("topright",  
       legend = c("Maschio", "Femmina"),
       col = c("blue", "red"),
       pch = 16,
       title = "Sesso",
       bty = "n",
       xpd = TRUE,  
       inset = c(-0.1, 0))

Si può osservare che non vi è molta distinzione dei valori della variabili per Sesso, ma in generale al crescere di lunghezza e Cranio aumenta anche il peso. Si osservano più valori bassi nel caso delle femmine.

In quest’ultimo grafico si vuole riportare l’interazione della variabile risposta Peso a secondo del numero di settimane di Gestazione e dell’essere o meno Fumatrici da parte delle donne.

dati2$Fumatrici=factor(dati2$Fumatrici)
suppressMessages(ggplot(data=dati2)+
  geom_point(aes(x=Gestazione, y=Peso), col="black")+
  geom_smooth(aes(x=Gestazione, y=Peso, color=Fumatrici), se=F, method="lm")+
  labs(title="Variazione del Peso secondo Gestazione e Fumatrici")
)
## `geom_smooth()` using formula = 'y ~ x'

Quello che si può vedere è che in generale nelle prime settimane non sono presenti esempi per la variabile Fumatrici. Il Peso cresce all’aumentare delle settimane di gestazione, ma in generale non sembra essere influenzato dalla variabile Fumatrici, difatti le due linee sembrano quasi sovrapporsi, si può solo notare che all’aumentare delle settimane il peso per le donne Fumatrici diminuisce leggermente.