Il progetto ha come finalità quella di prevedere il peso dei neonati alla nascita e individuare i fattori maggiormente a rischio. Questo permetterà al personale medico di intervenire in tempo in caso di anomalie e di prevedere in anticipo quali neonati potrebbero avere bisogno di cure intensive.

dati <- read.csv("neonati.csv",stringsAsFactors = T)
head(dati)
##   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
## 6         32            0         0         40 3200       495    340        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
## 6     osp2     F
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  : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Ospedale    : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
##  $ Sesso       : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...

Abbiamo un dataframe di 2500 osservazioni di 10 variabili.

Le variabili sono: gli anni della madre: variabile quantitativa continua; il numero di gravidanze avute da ogni madre: variabile quantitativa discreta; se la madre è fumatrice o meno: variabile qualitativa nominale; le settimane di gestazione: variabile quantitativa continua; il peso in grammi del neonato: variabile quantitativa continua; lunghezza in centimetri nel neonato: variabile quantitativa continua; diametro del cranio del neonato (in cm): variabile quantitativa continua; il tipo di parto, naturale o cesareo: variabile qualitativa nominale; l’ospedale nel quale deve avvenire il parto: variabile qualitativa nominale; sesso nel neonato: variabile qualitativa nominale.

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

Notiamo subito una anomalia negli anni della madre che non può essere 0. Allora eseguiamo un imputazione di quei dati per qui il valore di Anni.madre sia, per esempio, inferiore a 18.

attach(dati)

dati <- dati[Anni.madre >= 18, ]

summary(dati)
##    Anni.madre     N.gravidanze       Fumatrici         Gestazione   
##  Min.   :18.00   Min.   : 0.0000   Min.   :0.00000   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.00000   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000   Median :0.00000   Median :39.00  
##  Mean   :28.38   Mean   : 0.9947   Mean   :0.04231   Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.00000   3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000   Max.   :1.00000   Max.   :43.00  
##       Peso        Lunghezza         Cranio      Tipo.parto Ospedale   Sesso   
##  Min.   : 830   Min.   :310.0   Min.   :235.0   Ces: 712   osp1:803   F:1237  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330.0   Nat:1746   osp2:831   M:1221  
##  Median :3300   Median :500.0   Median :340.0              osp3:824           
##  Mean   :3285   Mean   :494.8   Mean   :340.1                                 
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350.0                                 
##  Max.   :4930   Max.   :565.0   Max.   :390.0

Ora vedremo tramite gli indici di correlazione e i grafici quali variabili sono maggiormente ppredittive del peso dei neonati alla nascita.

cor(Anni.madre,Peso)
## [1] -0.02247017
plot(Peso~Anni.madre)

Come possiamo notare le due variabili hanno un indice di correlazione negativo molto basso, ed infatti il relativo scatter plot presenta punti molto sparsi.

Vediamo le altre variabili.

cor(N.gravidanze,Peso)
## [1] 0.0024073

Non significativa.

cor(Fumatrici,Peso)
## [1] -0.01894534

Anche il fumo non è una buona indicazione del peso del nascituro, tuttavia è noto che non è salutare per il feto, e quindi richiederebbe ulterioni indagini statistiche, per esempio quanto neonati nascono senza problemi respiratori.

cor(Gestazione,Peso)
## [1] 0.5917687
plot(Peso~Gestazione)+
  abline(h=3000,col=2)

## integer(0)

Con un inddice di correlazione superiore a 0,5, Le settimane di gestazione, risultano significative.

cor(Lunghezza,Peso)
## [1] 0.7960368
plot(Peso~Lunghezza)+
  abline(h=3000,col=2)

## integer(0)

Assolutamente significativa.

cor(Cranio,Peso)
## [1] 0.7048015
plot(Peso~Cranio)+
  abline(h=3000,col=2)

## integer(0)

Significativa.

Usiamo il t-test per confrontare le medie della variabile Peso rispetto ai due gruppi della variabile Tipo.parto.

t.test(Peso~Tipo.parto,data=dati)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = -0.36254, df = 1446.5, p-value = 0.717
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -52.16367  35.88993
## sample estimates:
## mean in group Ces mean in group Nat 
##          3279.705          3287.842

Il p-value indica che non vi è una differenza significativa tra le medie dei due gruppi, e quindi la variabile Tipo.parto non è significatia.

Per la significatività dell’ospedale in cui avviene il parto possiamo usare il test ANOVA.

aov(Peso~Ospedale,data=dati)
## Call:
##    aov(formula = Peso ~ Ospedale, data = dati)
## 
## Terms:
##                  Ospedale Residuals
## Sum of Squares     902352 678537074
## Deg. of Freedom         2      2455
## 
## Residual standard error: 525.7279
## Estimated effects may be unbalanced
pairwise.t.test(Peso, Ospedale, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Peso and Ospedale 
## 
##      osp1 osp2
## osp2 1.00 -   
## osp3 0.33 0.33
## 
## P value adjustment method: bonferroni

Il p-value aggiustato è superiore in tutti e 3 i casi al valore α/3, e quindi non vi è una differenza significativa tra le medie del peso nei tre ospedali.

Concludiamo col Sesso.

t.test(Peso~Sesso,data=dati)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.07, df = 2450.8, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -289.1489 -208.3266
## sample estimates:
## mean in group F mean in group M 
##        3161.926        3410.663

Con un p-value estremamente piccolo rifiutiamo l’ipotesi di uuguaglianza tra i due sessi, quindi il sesso è significativo al fine del peso del neonato. Vediamolo anche graficamente.

boxplot(Peso~Sesso)

Dal boxplot osserviamo che il peso medio delle femmine supera di poco i 3 kg, Quello dei maschi si trova piu in alto sulla scala del peso.

Dal dataset risulta che in alcuni ospedali si conducono più parti cesarei che in altri, verifichiamolo usando il test chi-quadro.

creaiamo la tabella di contingenza.

tabella_cont = table(Tipo.parto,Ospedale)

eseguiamo il test chi-quadrato.

chisq.test(tabella_cont)
## 
##  Pearson's Chi-squared test
## 
## data:  tabella_cont
## X-squared = 1.0972, df = 2, p-value = 0.5778

Con un p-value che supera 0,5 non si rifiuta l’ipotesi nulla, e quindi non ci sono differenze sostanziali tra i tre ospedali per quanto riguarda il tipo di parto.

Adesso saggiamo l’ipotesi che le medie del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle dell’intera popolazione

Secondo il sito ufficiale dell’ospedale Bambinello Gesu, il peso medio dei neonati è di circa 3300 gr, con una differenza media di circa 150 gr tra maschi e femmine.

Applichiamo un t-test per la media di 3300.

t.test(Peso,
       mu = 3300,
       conf.level = 0.95,
       alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

Con un t-test che ricade nell’intervallo di accettazione, un p-value di circa il 12% ed un intervallo di confidenza che contiene il valore 3300 non si rifiuta l’ipotesi nulla, e quindi il peso medio, stabilito in 3284 gr non è significativamente diverso da quello dell’intera popolazione.

La lunghezza media dei maschi alla nascita è 50,5 cm, e quella delle femmine 49,5 cm, quindi possiamo usare la media di 500 mm nel t-test.

t.test(Lunghezza,
       mu = 500,
       conf.level = 0.95,
       alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

Questa volta si rifiuta l’ipotesi nulla, quindi sono necessarie ulteriori indagini statistiche.

Verifichiamo se le misure antropometriche, peso, lunghezza e dimensioni del cranio coincidono per i due sessi.

Abbiamo già verificato in precedenza che il peso per i due sessi e significativamente diverso. Rimangono da verificare la lunghezza e le dimensioni craniali.

boxplot(Lunghezza~Sesso)

boxplot(Cranio~Sesso)

Come si vede benissimo dai boxplot sopra, esistono differenza significative nelle misure antropometriche tra i due sessi.

Come abbiamo potuto notare le uniche variabili rilevanti per il peso finale dei neonati sono il periodo di gestazione, la lunghezza, le dimensioni del cranio nonché il sesso.Quindi costruiamo un modello di regressione multipla usando solo queste variabili.

mod_lin <- lm(Peso ~ Gestazione+Lunghezza+Cranio+Sesso)
summary(mod_lin)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16

Dalle caratteristiche del modello risulta un buon R^2, Quindi le variabili esplicative spiegano bene l’andamento della variabile risposta (ossia il peso).

Vediamo se il modello è migliorabile inserendo le interazioni tra variabili esplicative. Osserviamo la seguente matrice di correlazione per possibili rilevanti tra le variabili.

x11()

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

Notiamo che tutte le variabili di interesse sono correlate a due a due, quindi possiamo costruire altri modelli considerando le loro interazioni.

mod_lin2 <- update(mod_lin, ~.+Gestazione*Lunghezza)
summary(mod_lin2)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Gestazione:Lunghezza)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1122.0  -181.5   -14.0   166.3  2640.2 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -2.030e+03  9.216e+02  -2.202  0.02774 *  
## Gestazione           -9.318e+01  2.484e+01  -3.751  0.00018 ***
## Lunghezza             2.870e-02  2.030e+00   0.014  0.98872    
## Cranio                1.089e+01  4.246e-01  25.651  < 2e-16 ***
## SessoM                7.353e+01  1.121e+01   6.559 6.57e-11 ***
## Gestazione:Lunghezza  2.688e-01  5.302e-02   5.069 4.29e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.7 on 2494 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7283 
## F-statistic:  1341 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin3 <- update(mod_lin, ~.+Lunghezza*Cranio)
summary(mod_lin3)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Lunghezza:Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1139.05  -181.44   -15.46   163.32  2850.10 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.839e+03  1.019e+03  -1.804   0.0714 .  
## Gestazione        3.692e+01  3.951e+00   9.344  < 2e-16 ***
## Lunghezza        -2.013e-01  2.205e+00  -0.091   0.9273    
## Cranio           -4.409e+00  3.194e+00  -1.380   0.1676    
## SessoM            7.449e+01  1.121e+01   6.648 3.64e-11 ***
## Lunghezza:Cranio  3.113e-02  6.537e-03   4.763 2.01e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7286, Adjusted R-squared:  0.728 
## F-statistic:  1339 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin4 <- update(mod_lin, ~.+Gestazione*Cranio)
summary(mod_lin4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Gestazione:Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1125.48  -181.04   -13.99   163.92  2682.19 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -249.1238  1108.1125  -0.225  0.82214    
## Gestazione        -139.4557    29.5725  -4.716 2.54e-06 ***
## Lunghezza           10.4220     0.3010  34.620  < 2e-16 ***
## Cranio              -9.4246     3.4781  -2.710  0.00678 ** 
## SessoM              73.3078    11.1830   6.555 6.72e-11 ***
## Gestazione:Cranio    0.5262     0.0904   5.821 6.62e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7292 
## F-statistic:  1347 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin5 <- update(mod_lin, ~.+Sesso*Lunghezza)
summary(mod_lin5)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Lunghezza:Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.94  -183.34   -18.14   163.27  2564.38 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -6472.4840   162.8684 -39.741   <2e-16 ***
## Gestazione          31.5210     3.7854   8.327   <2e-16 ***
## Lunghezza            9.8390     0.3532  27.857   <2e-16 ***
## Cranio              10.6445     0.4244  25.081   <2e-16 ***
## SessoM            -341.5613   213.3117  -1.601   0.1095    
## Lunghezza:SessoM     0.8492     0.4300   1.975   0.0484 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7265, Adjusted R-squared:  0.726 
## F-statistic:  1325 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin6 <- update(mod_lin, ~.+Sesso*Cranio)
summary(mod_lin6)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Cranio:Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1140.34  -184.92   -16.89   163.91  2639.72 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6545.4176   170.3532 -38.423   <2e-16 ***
## Gestazione       31.4577     3.7898   8.301   <2e-16 ***
## Lunghezza        10.2060     0.3007  33.942   <2e-16 ***
## Cranio           10.3353     0.5359  19.284   <2e-16 ***
## SessoM         -157.7067   231.5384  -0.681    0.496    
## Cranio:SessoM     0.6960     0.6797   1.024    0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2494 degrees of freedom
## Multiple R-squared:  0.7262, Adjusted R-squared:  0.7257 
## F-statistic:  1323 on 5 and 2494 DF,  p-value: < 2.2e-16

#E gli affetti quadratici.

mod_lin7 <- update(mod_lin, ~.+I(Gestazione^2))
summary(mod_lin7)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1132.84  -181.45   -15.99   162.62  2649.78 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4643.6094   899.4989  -5.162 2.63e-07 ***
## Gestazione        -80.7957    49.7863  -1.623   0.1047    
## Lunghezza          10.3087     0.3039  33.920  < 2e-16 ***
## Cranio             10.7663     0.4262  25.259  < 2e-16 ***
## SessoM             76.9359    11.2436   6.843 9.75e-12 ***
## I(Gestazione^2)     1.4960     0.6627   2.258   0.0241 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7267, Adjusted R-squared:  0.7261 
## F-statistic:  1326 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin8 <- update(mod_lin, ~.+I(Lunghezza^2))
summary(mod_lin4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Gestazione:Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1125.48  -181.04   -13.99   163.92  2682.19 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -249.1238  1108.1125  -0.225  0.82214    
## Gestazione        -139.4557    29.5725  -4.716 2.54e-06 ***
## Lunghezza           10.4220     0.3010  34.620  < 2e-16 ***
## Cranio              -9.4246     3.4781  -2.710  0.00678 ** 
## SessoM              73.3078    11.1830   6.555 6.72e-11 ***
## Gestazione:Cranio    0.5262     0.0904   5.821 6.62e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7292 
## F-statistic:  1347 on 5 and 2494 DF,  p-value: < 2.2e-16
mod_lin9 <- update(mod_lin, ~.+I(Cranio^2))
summary(mod_lin9)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1127.21  -183.53   -14.88   163.98  2610.62 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   76.70852 1153.12221   0.067    0.947    
## Gestazione    37.73144    3.91778   9.631  < 2e-16 ***
## Lunghezza     10.44511    0.30147  34.647  < 2e-16 ***
## Cranio       -31.41831    7.17690  -4.378 1.25e-05 ***
## SessoM        74.30205   11.16712   6.654 3.50e-11 ***
## I(Cranio^2)    0.06226    0.01060   5.875 4.80e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7293 
## F-statistic:  1347 on 5 and 2494 DF,  p-value: < 2.2e-16

Come possiamo osservare non ci sono differenze significative, quindi effettuiamo ulteriori analisi statistiche al fine di scegliere il modello più conveniente, Come il test ANOVA, per l’analisi della varianza o il criterio di informazione Bayesiano (BIC). Per applicare il test ANOVA è necessario che le variabili rispettino alcune condizioni, come la distribuzione normale della variabile risposte, la similitudine tra le varianze e l’indipendenza delle osservazioni. Verifichiamo la normalità con un test di Shapiro-Wilk.

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

Come si nota si rifiuta l’ipotesi di normalità della variabile peso, quindi non utilizzeremo il test ANOVA.

Procediamo allora col BIC.

BIC(mod_lin,mod_lin2,mod_lin3,mod_lin4,mod_lin5,mod_lin6,mod_lin7,mod_lin8,mod_lin9)
##          df      BIC
## mod_lin   6 35220.54
## mod_lin2  7 35202.74
## mod_lin3  7 35205.73
## mod_lin4  7 35194.64
## mod_lin5  7 35224.46
## mod_lin6  7 35227.32
## mod_lin7  7 35223.27
## mod_lin8  7 35138.48
## mod_lin9  7 35194.01

Seguendo la logica del metodo BIC quindi sceglieremo l’ottavo modello.

Analiziamo ora la parte erratica, per individuare eventuali valori leverage o outliers che potrebbero distorcerne le previsioni.

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

I residui seguono una distribuzione approssimativamente normale e non ci sono valori leverage o outliers che superano la soglia di allarme.

shapiro.test(residuals(mod_lin8))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_lin8)
## W = 0.98552, p-value = 2.86e-15

Il fatto che il test di shapiro-Wilk bocci l’ipotesi di normalità può essere dovuto al fatto che la loro distribuzione di densità ha code più allungate, o è asimmetrica o ancora lepto/platicurtica. Verifichiamolo.

plot(density(residuals(mod_lin8)))

Le altre condizioni fondamentali che i residui di un buon modello di regressione devono rispettare sono quella di omoschedasticità e quella di non correlazione, verifichiamoli tramite i test di Breusch-Pagan e quello di Durbin-Watson, rispettivamente.

lmtest::bptest(mod_lin8)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_lin8
## BP = 126.16, df = 5, p-value < 2.2e-16
lmtest::dwtest(mod_lin8)
## 
##  Durbin-Watson test
## 
## data:  mod_lin8
## DW = 1.9508, p-value = 0.1093
## alternative hypothesis: true autocorrelation is greater than 0

Purtroppo viene rifiutata anche l’ipotesi di omoschedasticità, quindi è un modello da usare con cautela.

Accertiamoci in fine che il nostro modello non presenti probleni di multicolinearita tra mite il fattore di inflazione della varianza (VIF).

car::vif(mod_lin8)
##     Gestazione      Lunghezza         Cranio          Sesso I(Lunghezza^2) 
##       1.781859     237.696630       1.607724       1.044426     229.653852

Osserviamo che compaiono dei VIF molto alti, ma è plausibile dati gli effetti di interazione e l’effetto quadratico della variabile Lunghezza.

Facciamo dunque delle previsioni col nostro modello. Per esempio di una femmina nata dopo la 39°settimana di gestazione.

predict(mod_lin8, newdata = data.frame(Gestazione=39,Sesso="F",
                                Lunghezza=mean(Lunghezza),Cranio=mean(Cranio)))
##        1 
## 3227.744

Per finire diamo una rappresentazione grafica dell’andamento del peso col periodo di gestazione.

library(ggplot2)

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso),position="jitter")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso),se=F, method="lm")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso), col="black",se=F, method="lm")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

come abbiamo già notato prima, la retta che rappresenta i maschi si trova più in alto rispetto a quella che rappresenta le femmine, proprio perché pesano mediamente di più. Inoltre la retta di regressione passa molto bene in mezzo alla nuvola di punti, il che ci conferma che il nostro è un buon modello predittivo.

Il progetto di previsione del peso neonatale è un’iniziativa fondamentale per Neonatal Health Solutions. Attraverso l’uso di dati clinici dettagliati e strumenti di analisi statistica avanzati, possiamo contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere. Questo progetto rappresenta un punto di svolta per l’azienda, consentendo non solo un miglioramento della pratica clinica ma anche l’implementazione di politiche sanitarie più informate e proattive.