Modello Statistico per la Previsione del Peso Neonatale

library(knitr)
knitr::opts_chunk$set(echo = TRUE)

neonati <- read.csv("C:/Users/ce168diedifi/Desktop/Profession.AI/Corso Statistica Inferenziale/Progetto/neonati.csv", sep=",", header=TRUE, encoding="UTF-8")
attach(neonati)
kable(head(neonati))

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

options(warn = -1)

#install.packages("kableExtra")
library("knitr")
library("kableExtra")
library(dplyr)

## 
## Caricamento pacchetto: 'dplyr'

## Il seguente oggetto è mascherato da 'package:kableExtra':
## 
##     group_rows

## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union

neonati.factor <- neonati %>%
  mutate(across(where(is.character), as.factor))

#str(neonati.factor)

kable(summary(neonati.factor)  , format = "html") %>%
  kable_styling(font_size = 12, full_width = FALSE, position = "center",html_font = "Consolas") %>%
  row_spec(0, bold = TRUE, font_size = 14, color = "white", background = "silver") %>%
  column_spec(1, italic = TRUE, color = "blue") %>%
  row_spec(0, bold = TRUE, font_size = 12, color = "black")  %>%
  row_spec(1:nrow(summary(neonati)), background = "white")

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

Guardando il summary del dataset noto dei possibili valori anomali per la variabile “Anni.madre” che ha dei valori minimi a zero e per la variabile “N.gravidanze” dove c’è un valore massimo di 12. Lo vedremo con l’analisi di dettaglio.

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie. Un focus particolare sarà dato alla relazione tra le variabili materne (età, fumo, gravidanze precedenti) e il peso del neonato.

Descrivo le singoli variabili presenti nel dataset:

• Peso –> variabile quantitativa continua
  
• Anni.madre –> variabile quantitativa continua
  
• N.gravidanze –> variabile quantitativa discreta
  
• Gestazione –> variabile quantitativa continua
  
• Lunghezza –> variabile quantitativa continua
  

• Cranio –> variabile quantitativa continua

• Tipo.parto –> variabile qualitativa nominale
  
• Ospedale –> variabile qualitativa nominale
  
• Sesso –> variabile qualitativa nominale
  
• Fumatrici –> variabile qualitativa nominale

Rappresento con uno scatterplot la relazione tra le variabili:

• Peso-Anni.Madre
• Peso-N.gravidanze
• Peso-Fumatrici
• Peso-Gestazione

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

Dal grafico, il peso dei neonati sembra aggirarsi intorno ai 3kg / 4kg quando gli anni della madre sono tra i 20e i 40 anni al momento della gestazione. La relazione sembra essere lineare e piuttosto orizzontale. Lungo la distribuzione e fuori noto degli ouliers.

kable(cor(Peso,Anni.madre))

x
-0.0224702

il coefficiente di correlazione di Bravais-Pearson è di poco sotto lo zero suggerisce che le due variabili non presentano una relazione lineare significativa e la variazione di una variabile non influenza apprezzabilmente l’altra.

Questa è la rappresentazione della correlazione tra Peso e il numero di gravidanze avute precedeti a quella corrente

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

kable(cor(Peso,N.gravidanze))

x
0.0024073

Qui sia grafico che indice di correlazione indica una correlazione praticamente nulla tra le due variabili

plot(Peso,Gestazione,pch=20)

kable(cor(Peso,Gestazione))

x
0.5917687

Invece visualizzando e calcolando la correlazione tra Peso-Gestazione si nota una correlazione positiva moderata-forte tra le variabili. c’è una tendenza lineare evidente, ma non perfetta perchè lo scatterplot sembra assumere una tendenza leggermente curva.

Creazione del Modello di Regressione Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato. Ad esempio, ci aspettiamo che il fumo materno e le gravidanze multiple abbiano effetti negativi sul peso alla nascita, mentre una maggiore durata della gestazione potrebbe aumentare il peso del neonato.

n <- nrow(neonati)
neonati_numeric <- neonati[sapply(neonati, is.numeric)]
kable( round(cor(neonati_numeric), 2) )

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio
Anni.madre	1.00	0.38	0.01	-0.14	-0.02	-0.06	0.02
N.gravidanze	0.38	1.00	0.05	-0.10	0.00	-0.06	0.04
Fumatrici	0.01	0.05	1.00	0.03	-0.02	-0.02	-0.01
Gestazione	-0.14	-0.10	0.03	1.00	0.59	0.62	0.46
Peso	-0.02	0.00	-0.02	0.59	1.00	0.80	0.70
Lunghezza	-0.06	-0.06	-0.02	0.62	0.80	1.00	0.60
Cranio	0.02	0.04	-0.01	0.46	0.70	0.60	1.00

Qui voglio verificare se le variabili del dataset sono indipendenti, quindi testiamo se accettare o rifiutare l’ipotesi nulla che le variabili sono indipendenti,non esiste relazione tra di esse.

test.indipendenza <- chisq.test(neonati_numeric)
test.indipendenza

## 
##  Pearson's Chi-squared test
## 
## data:  neonati_numeric
## X-squared = 43525, df = 14994, p-value < 2.2e-16

il p-value è molto basso quindi si rifiuta l’ipotesi nulla e possiamo affermare che le variabili sono correlate.

Iniziamo col visulizzare la matrice di correlazione tra le variabili del dataset per avere una situazione generale delle correlazioni più probabili,cosa che comunque va verificata.

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

L’obiettivo è quello di identificare quali delle variabili del dataset sono più predittive del Peso alla nascita. Un focus particolare lo avremo sulle variabili “Anni.Madre”, “N.gravidanze”, “Fumatrici”, “Gestazione”.

Creo il primo modello di regressione lineare multiplo:

mod1 <- lm( data = neonati, Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso )
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1141.25  -181.75   -14.87   160.54  2629.75 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6740.6584   141.3921 -47.674   <2e-16 ***
## Anni.madre        0.9543     1.1335   0.842   0.3999    
## N.gravidanze     11.5753     4.6656   2.481   0.0132 *  
## Fumatrici       -31.5916    27.5729  -1.146   0.2520    
## Gestazione       32.8814     3.8227   8.602   <2e-16 ***
## Lunghezza        10.2718     0.3010  34.128   <2e-16 ***
## Cranio           10.4828     0.4266  24.573   <2e-16 ***
## Tipo.partoNat    30.2808    12.0990   2.503   0.0124 *  
## SessoM           78.0278    11.1921   6.972    4e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7279, Adjusted R-squared:  0.727 
## F-statistic:   833 on 8 and 2491 DF,  p-value: < 2.2e-16

Il primo modello che vado a realizzare contiene tutte le sue variabili. Questo modello ha un r quadro aggiustato molto vicino a 1, sembra essere un buon modello ma nello specifico alcune delle variabili presenti hanno una bassa significatività e che quindi vanno eliminate durante lo stepwise.

Quindi analizzando questo modello per ogni singola variabile:

• “Anni.madre” per ogni incremento unitario di questa variabile l’impatto sul peso del neonato è in media di 0.95 grammi e il p-value è > del 5%. ININFLUENTE
• “N.gravidanze” per ogni incremento unitario di questa variabile l’impatto sul peso del neonato è in media di 11.57 grammi e il p-value < del 5%. INFLUENTE.
• “Fumatrici” per ogni incremento unitario di questa variabile l’impatto sul peso del neonato è in media di -31.59 grammi e il un p-value molto superiore al 5%. ININFLUENTE
• “Gestazione” per ogni incremento unitario di questa variabile l’impatto sul peso del neonato è in media di 33 grammi e il p-value molto basso. INFLUENTE.
• “Lunghezza” QUESTA VARIABILE VA TENUTA NEL MODELLO.
• “Cranio” QUESTA VARIABILE VA TENUTA NEL MODELLO.
• “Sesso” ha un impatto stimato sul peso del neonato di 78 grammi e il p-value è molto basso. QUESTA VARIABILE VA TENUTA NEL MODELLO.
• “Tipo.parto” non è un informazione predittiva per il peso del neonato. ININFLUENTE

Anche se alcune variabili come “Anni.Madre” e “Fumatrici” sono risultate di bassa influenza proverò a inserirle nei modelli futuri.

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

Inizio ad eliminare la variabile “Tipo.Parto”:

mod2 <- update(mod1,~.-Tipo.parto )
summary(mod2)

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

il risultato è un modello pressocchè invariato.

Elimino la variabile “Anni.madre”:

mod3 <- update(mod2,~.-Anni.madre )
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = neonati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1150.3  -181.3   -15.7   163.0  2636.3 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.6714   135.7178 -49.232  < 2e-16 ***
## N.gravidanze    12.7185     4.3450   2.927  0.00345 ** 
## Fumatrici      -30.4634    27.5948  -1.104  0.26972    
## Gestazione      32.5914     3.8051   8.565  < 2e-16 ***
## Lunghezza       10.2341     0.3009  34.011  < 2e-16 ***
## Cranio          10.5359     0.4262  24.718  < 2e-16 ***
## SessoM          78.1713    11.2028   6.978 3.83e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16

Il valore di r2a è leggermente cresciuto. Nonostante le credenze popolari sull’impatto del fumo, non ci sono prove sufficienti per affermare che la variabile “Fumatrici” influisce sul Peso del neonato.

Elimino la variabile “Fumatrici”:

mod4 <- update(mod3,~.-Fumatrici )
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Qui abbiamo un r2a identico al precedente, ma con tutte le variabili presenti sono molto significative perchè il singolo p-value è molto sotto il 5%.

Provo l’effetto quadratico sulla variabile “Fumatrici” per verificare se effettivamente è una variabile senza effetto:

mod5 <- update(mod4, ~.+I(Fumatrici^2))
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Fumatrici^2), data = neonati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1150.3  -181.3   -15.7   163.0  2636.3 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6681.6714   135.7178 -49.232  < 2e-16 ***
## N.gravidanze      12.7185     4.3450   2.927  0.00345 ** 
## Gestazione        32.5914     3.8051   8.565  < 2e-16 ***
## Lunghezza         10.2341     0.3009  34.011  < 2e-16 ***
## Cranio            10.5359     0.4262  24.718  < 2e-16 ***
## SessoM            78.1713    11.2028   6.978 3.83e-12 ***
## I(Fumatrici^2)   -30.4634    27.5948  -1.104  0.26972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16

Anche l’effetto quadrativo sulla variabile “Fumatrici” conferma la sua non influenza. Pertanto per adesso il modello 4 resta il migliore.

Provo con far interagire la variabile “Gestazione” con “Sesso”:

mod6 <- lm(data = neonati, Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + Gestazione*Sesso)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione * Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1144.73  -181.02   -15.21   163.38  2638.22 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -6571.7016   166.9616 -39.361  < 2e-16 ***
## N.gravidanze         12.4289     4.3396   2.864  0.00422 ** 
## Gestazione           29.4475     4.5818   6.427 1.55e-10 ***
## Lunghezza            10.2523     0.3006  34.102  < 2e-16 ***
## Cranio               10.5416     0.4262  24.732  < 2e-16 ***
## SessoM             -185.9051   234.7635  -0.792  0.42850    
## Gestazione:SessoM     6.7624     6.0090   1.125  0.26054    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16

Aggiungendo al modello l’iterazione tra Gestazione-Sesso questa non è significativa ma anche il Sesso perde di significatività.

Pertanto il modello 4 resta ancora il migliore.

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

BIC(mod1, mod2, mod3, mod4, mod5, mod6)

##      df      BIC
## mod1 10 35235.35
## mod2  9 35233.81
## mod3  8 35226.70
## mod4  7 35220.10
## mod5  8 35226.70
## mod6  8 35226.65

Il modello 4, dell’analisi sopra, sembra essere di nuovo il modello migliore.

Provo ad eseguire una procedura di selezione stepwise dei modelli in autoamatico, per verificare la mia analisi con quella automatica:

#install.packages("MASS")
library("MASS")

## 
## Caricamento pacchetto: 'MASS'

## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     select

stepwise.mod <- MASS::stepAIC(mod2, direction = "both")

## Start:  AIC=28084.7
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     53803 187973654 28083
## - Fumatrici     1     90879 188010731 28084
## <none>                      187919851 28085
## - N.gravidanze  1    439797 188359648 28089
## - Sesso         1   3662833 191582684 28131
## - Gestazione    1   5585184 193505035 28156
## - Cranio        1  45791026 233710877 28628
## - Lunghezza     1  87220887 275140738 29036
## 
## Step:  AIC=28083.42
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     91892 188065546 28083
## <none>                      187973654 28083
## + Anni.madre    1     53803 187919851 28085
## - N.gravidanze  1    646039 188619694 28090
## - Sesso         1   3671289 191644943 28130
## - Gestazione    1   5531705 193505359 28154
## - Cranio        1  46066755 234040410 28629
## - Lunghezza     1  87218857 275192512 29034
## 
## Step:  AIC=28082.64
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188065546 28083
## + Fumatrici     1     91892 187973654 28083
## + Anni.madre    1     54816 188010731 28084
## - N.gravidanze  1    623141 188688687 28089
## - Sesso         1   3655292 191720838 28129
## - Gestazione    1   5464853 193530399 28152
## - Cranio        1  46108583 234174130 28629
## - Lunghezza     1  87632762 275698308 29037

summary(stepwise.mod)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Il modello selezionato in automatico è identico al mio modello 4 che si conferma essere nuovamente il migliore.

#install.packages("car")
library(car)

## Caricamento del pacchetto richiesto: carData

## 
## Caricamento pacchetto: 'car'

## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     recode

vif(mod4)

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

Infatti calcolando i VIF delle variabili del modello 4 tutti i VIF molto < di 5.

Se confronto i due modelli, modello 3 e modello calcolato con mass, questi sono identici.

BIC(mod4, stepwise.mod)

##              df     BIC
## mod4          7 35220.1
## stepwise.mod  7 35220.1

Adesso verifichiamo i test di normalità della variabile “Peso”:

shapiro.test(Peso)

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

Il test di Shapiro Wilk indica che i dati non seguono proprio una distribuzione normale, ma essendo il valore W molto vicino a 1 l’asimmetria potrebbe non essere così accentuata.

Adesso verifico la direzione dell’asimmetria:

library(moments)
skewness(Peso)

## [1] -0.6470308

Abbiamo una distribuzione asimmetrica con coda a sinistra (negativa).

Calcolando la Curtosi verifichiamo la “pesantezza delle code”:

kurtosis(Peso)

## [1] 5.031532

con un valore > 3 abbiamo una distribuzione leptocurtica.

Rappresento graficamente la distribuzione del “Peso”:

library(ggplot2)
ggplot(data = neonati, aes(x = Peso)) +
geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "blue", alpha = 0.6) +
geom_density(color = "red", lwd = 1)

Rappresentanto la variabile “Peso” col bocplot evidenzio la presenza di outliers:

ggplot(data = neonati, aes(y = Peso)) +
  geom_boxplot(fill = "lightblue")

Infatti il boxplot evidenzia degli outliers nella parte inferiore. COnsidrando che stiamo analizzando la variabile “Peso” dei neonati, questi outliers potrebbero essere dovuti a casi estremi come bambini sottopeso, nascite premature, morte prematura del feto, altro!

Di seguito l’analisi dei residui:

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

Nel primo grafico i dati si sono raccolti intorno alla media di zero ma il pattern e’ curvato e la concentrazione dei dati e’ chiaramente verso i valori piu’ alti del peso.

Nel secondo grafico i residui seguono una distribuzione normale perchè tutti raccolti intorno alla bisettrice dove si notano dei residui sul lato sinistro e sul lato destro con un osservazione 1551 molto al di fuori della tendenza.

Nel terzo grafico abbiamo una nuvola di punti casuale e quasi orizzontale intorno ad un valore di Y che indica una varianza costante. Anche qui abbiamo una curvatura ma non è troppo eccessiva.

Nel quarto grafico si evidenziano i potenziali valori influenti, cioè quei valori outliers e/o leverage o combinazioni dei due. Qui si nota solo un valore che supera la soglia di avvertimento di 0.5 che è il 1551.

Calcolo i valori di leverage:

lev <- hatvalues(mod4)
p <- sum(lev)
n <- nrow(neonati)
soglia <- 2* p/n
plot(lev)
abline(h=soglia,col=2)

Questa visualizzazione evidenzia la presenza di parecchi leverage, tutti quelli al di sopra della linea rossa.

Scrivo tutti i leverage individuati:

lev[lev>soglia]

##          13          15          34          67          89          96 
## 0.005630918 0.007050422 0.006744143 0.005891590 0.012816470 0.005351586 
##         101         106         131         134         151         155 
## 0.007526751 0.014479871 0.007229339 0.007552819 0.010883937 0.007207682 
##         161         189         190         204         205         206 
## 0.020335483 0.004893297 0.005366557 0.014490604 0.005351634 0.009476652 
##         220         294         305         310         312         315 
## 0.007393997 0.005912765 0.005442061 0.028812123 0.013169272 0.005385800 
##         378         440         442         445         486         492 
## 0.015934061 0.005404736 0.007723662 0.007509382 0.005164446 0.008274018 
##         497         516         582         587         592         614 
## 0.005166306 0.013079851 0.011665555 0.008412325 0.006384116 0.005299262 
##         638         656         657         684         697         702 
## 0.006688287 0.005927777 0.005322685 0.008818987 0.005863826 0.005202259 
##         729         748         750         757         765         805 
## 0.005023115 0.008565543 0.006942097 0.008145491 0.006070298 0.014356657 
##         828         893         895         913         928         946 
## 0.007179817 0.005075205 0.005295896 0.005571144 0.022742332 0.006909044 
##         947         956         985        1008        1014        1049 
## 0.008409465 0.007784123 0.007039416 0.005343037 0.008470133 0.004956169 
##        1067        1091        1106        1130        1166        1181 
## 0.008465430 0.008933360 0.005967317 0.031728597 0.005513559 0.005677676 
##        1188        1200        1219        1238        1248        1273 
## 0.006477203 0.005492370 0.030694311 0.005908078 0.014622914 0.007085831 
##        1291        1293        1311        1321        1325        1356 
## 0.006117497 0.006073639 0.009625908 0.009293111 0.004857169 0.005303442 
##        1357        1385        1395        1400        1402        1411 
## 0.006965051 0.012636943 0.005126697 0.005925069 0.004811441 0.008048184 
##        1420        1428        1429        1450        1505        1551 
## 0.005155654 0.008192811 0.021757172 0.015104831 0.013330439 0.048769569 
##        1553        1556        1573        1593        1606        1610 
## 0.008504889 0.005919673 0.005047204 0.005623758 0.005001812 0.008722184 
##        1617        1619        1628        1686        1693        1701 
## 0.004866796 0.015067498 0.005069731 0.009349313 0.005077858 0.010842957 
##        1712        1718        1727        1735        1780        1781 
## 0.006992084 0.006958857 0.013300523 0.004884846 0.025538678 0.016832361 
##        1809        1827        1868        1892        1962        1967 
## 0.008707504 0.006065698 0.005205637 0.005332812 0.005540442 0.005337356 
##        1977        2037        2040        2046        2086        2089 
## 0.006927281 0.004889127 0.011494872 0.005471670 0.013193090 0.006293550 
##        2098        2114        2115        2120        2140        2146 
## 0.005094455 0.013316875 0.011772090 0.018659995 0.006244232 0.005802168 
##        2148        2149        2157        2175        2200        2215 
## 0.007926839 0.013583436 0.005907225 0.032527273 0.011670024 0.004892265 
##        2216        2220        2221        2224        2225        2244 
## 0.008117864 0.005414040 0.021628717 0.005838076 0.005591261 0.006929217 
##        2257        2307        2317        2318        2337        2359 
## 0.006170254 0.013965608 0.007673614 0.004831118 0.005230450 0.010067364 
##        2408        2422        2436        2437        2452        2458 
## 0.009696691 0.021532808 0.004986522 0.023943328 0.023838489 0.008506087 
##        2471        2478 
## 0.020903740 0.005775173

Queste sono le osservazioni che si trovano lontano dalla soglia dei regressori.

Vediamo gli outliers:

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

Anche qui gli outliers sono parecchi, ma andiamo a verificare con un test t quali sono gli effettivi outliers presenti nel modello:

car::outlierTest(mod4)

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

Con questo test si evidenzia che gli effettivi outliers sono soltanto 3 perchè gli altri non sono considerati eccessivamente lontani dal resto della nuvola di punti.

Verifichiamo la distanza di cook

cook <- cooks.distance(mod4)
plot(cook, main = "Distanza di Cook", type = "h", ylab = "Distanza di Cook")
abline(h = 1, col = "red", lty = 2)

Il grafico mostra un valore massimo della distanza di 0.8

kable(max(cook))

x
0.8300965

Il valore max mostra un valore di 0.83 che come abbiamo visto sopra supera la soglia di allarme ma non quella di pericolo per poter essere considerato un valore influente per la bontà del modello.

Altri test:

library(lmtest)

## Caricamento del pacchetto richiesto: zoo

## 
## Caricamento pacchetto: 'zoo'

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     as.Date, as.Date.numeric

SUMMARY.RESIDUALS <- data.frame(
  "Breusch-Pagan test"  = sprintf("%.3e", unname(bptest(mod4)$p.value)),
  "Durbin-Watson test"  = sprintf("%.3e", unname(dwtest(mod4)$p.value)),
  "Shapiro-Wilk normality test" = sprintf("%.3e", unname(shapiro.test(residuals(mod4))$p.value)),
  check.names = FALSE
)

kable(SUMMARY.RESIDUALS)

Breusch-Pagan test	Durbin-Watson test	Shapiro-Wilk normality test
5.946e-18	1.224e-01	6.782e-21

I test statistici sui residui indicano:

• Presenza di eteroschedasticità
  
• Non ci sono prove di autocorrelazione
• I residui potrebbero non seguire un andamento normale

plot(density(residuals(mod4)))

anche se i test statistici sembrano suggerire un andamento non normale dei residui, la rappresentazione grafica si avvicina moltissimo a una normale, quindi non mi sento di bocciare il modello 4 nell’utilizzarlo per le previsioni.

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

Lunghezza.media <- mean(neonati$Lunghezza, na.rm = TRUE)
Cranio.media <- mean(neonati$Cranio, na.rm = TRUE)

predict(mod4
        , newdata = data.frame(N.gravidanze = 3
                               , Fumatrici = 0
                               , Gestazione = 39
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               , Sesso = "M"
                               )
        ,interval = "confidence" )

##        fit      lwr      upr
## 1 3349.083 3326.274 3371.892

Coi parametri forniti possiamo fare due stime, la prima è nel caso il bambino sia un maschio. In questo caso la previsione del Peso sarebbe intorno ai 3Kg e 350g.

predict(mod4
        , newdata = data.frame(N.gravidanze = 3
                               , Fumatrici = 0
                               , Gestazione = 39
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               , Sesso = "F"
                               )
        ,interval = "confidence" )

##       fit      lwr      upr
## 1 3271.09 3247.763 3294.416

Nel caso di una bambino Femmina la previsione del Peso sarebbe intorno ai 3Kg e 270g.

4. Visualizzazioni Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

• Impatto del numero di settimane di gestazione sul peso del neonato:

library(ggplot2)

ggplot(data=neonati)+
 geom_point( aes( x=Gestazione, y=Peso, col=Sesso ), position = "jitter" )+
geom_smooth( aes( x=Gestazione, y=Peso, col=Sesso ), se=F, method = "lm", formula=y ~ poly(x, 2) )+
geom_smooth( aes( x=Gestazione, y=Peso, col="Aggregato" ), se=F, method = "lm", formula=y ~ poly(x, 2) )

Il numero di settimane di “Grestazione” ha un impatto positivo sul peso del neonato per entrambi i sessi.

• Impatto di una madre fumatrice sul peso del neonato:

library(ggplot2)
ggplot(data=neonati)+
 geom_point( aes( x=Fumatrici, y=Peso, col=Sesso ), position = "jitter" )+
geom_smooth( aes( x=Fumatrici, y=Peso, col=Sesso ), se=F, method = "lm", formula=y ~ x )+
geom_smooth( aes( x=Fumatrici, y=Peso, col="Aggregato" ), se=F, method = "lm", formula=y ~ x )

La maggior parte delle osservazioni in questo dataset contiene informazioni di madri non fumatrici. Interpretando questo grafico mi viene da supporre che l’essere una fumatrice non ha un impatto significativo sul Peso del neonato per entrambi i sessi.

• Impatto del numero di gravidanze sul peso del neonato:

library(ggplot2)
ggplot(data=neonati)+
 geom_point( aes( x=N.gravidanze, y=Peso, col=Sesso ), position = "jitter" )+
geom_smooth( aes( x=N.gravidanze, y=Peso, col=Sesso ), se=F, method = "lm", formula=y ~ x )+
geom_smooth( aes( x=N.gravidanze, y=Peso, col="Aggregato" ), se=F, method = "lm", formula=y ~ x )

Leggendo il grafico da sx a dx, la nuvola di punti sembra formare un triangolo, questo fa supporre che con l’aumentare delle gravidanze il peso del neonato non riesca ad arrivare ai 3kg. Quindi si il numero di gravidanze ha un impatto sul Peso del neonato per entrambi i sessi.