Modello Statistico per la Previsione del Peso Neonatale

1. Raccolta dei Dati e Struttura del Dataset

dati <- read.csv("neonati.csv", stringsAsFactors = T)

2. Analisi e Modellizzazione

2.1 Analisi Preliminare

Visualizzazione dei dati osservati per ogni varibiale con distribuzione di frequenza delle variabili qualitative o quantitative discrete.

attach(dati)
library(knitr)
kable(summary(dati))
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 
indice.da.elimiare <- which(dati$Anni.madre %in% c(0, 1))
dati <- dati[-indice.da.elimiare, ]


boxplot(dati$Peso ~ dati$Sesso )  # Outlier nel peso

hist(dati$Gestazione,
     main = "Distribuzione della variabile Gestazione",
     xlab = "Numero di settimane di gestazione",               
     ylab = "Frequenza",
     col = "blue",
     xlim = c(20, 45),
     ylim = c(0, 800))

analisi_vs_peso <- function(df, variabili, risposta = "Peso") {
  risultati <- list()
  
  for (var in variabili) {
    plot(df[[var]], df[[risposta]],
         pch = 20, col = "steelblue",
         main = paste("Relazione tra", var, "e", risposta),
         xlab = var, ylab = risposta)
    
    abline(lm(df[[risposta]] ~ df[[var]]), col = "red") 
    

    corr <- cor(df[[var]], df[[risposta]], use = "complete.obs")
    risultati[[var]] <- corr
  }
  
  return(risultati)
}

ris <- analisi_vs_peso(dati, c("Gestazione", "Lunghezza", "Anni.madre", "Cranio"))

ris
## $Gestazione
## [1] 0.5919592
## 
## $Lunghezza
## [1] 0.7960415
## 
## $Anni.madre
## [1] -0.02378138
## 
## $Cranio
## [1] 0.7048438

Analisi:

  1. Dal Summary emerge che per la variabile Anni.madre ci sono valori errati che hanno un notevole impatto sugli indicatori (media ad esempio). Cibsuderando che le osservazioni errate sono 2 su 2500 si è deciso di eliminare le righe anomale dalla base dati.
  2. Il boxplot della variabile risposta rispetto, differenziato per sesso, mostra la presenza di alcuni outliers e sopratutto mostra una differenza di peso medio fra i due sessi che bisogna considerare durante l’analisi.
  3. Dall’analisi grafica e dei coefficienti di correlazione fra il Peso e le variabili Gestazione, Anni.madre, Lunghezza e Cranio si evince una correllazione positiva della variabile risposta con il numero di settimane di Gestazione e con la lunghezza e il le dimensioni del cranio. Non si evince nessuna correlazione fra Peso ed Età della madre e tale variabile può ritenersi non utile ai fini dell’analisi.
  4. Si analizzerà nelle successive fasi eventuali correlazioni fra le variabili antropometriche che possano inficiare l’analisi (lunghezza, cranio).

Ipotesi 1:in alcuni ospedali si fanno più parti cesarei

attach(dati)
newdf <- data.frame(Tipo.parto, Ospedale)
Tabella_frequenze <- table(newdf$Ospedale , newdf$Tipo.parto)


#Ipotesi nulla ((H_0)): La proporzione di parti cesarei è uguale in tutti gli ospedali.
#Ipotesi alternativa ((H_1)): La proporzione di parti cesarei è diversa tra almeno un ospedale e gli altri.

chisq.test(Tabella_frequenze)
## 
##  Pearson's Chi-squared test
## 
## data:  Tabella_frequenze
## X-squared = 1.083, df = 2, p-value = 0.5819

Esito: considerando che il p-value è maggiore di 0,05 non si rifiuta l’ipotesi nulla e non possiamo dedurre che ci siano differenze significative di tipologia di parto fra gli ospedali.

Ipotesi 2. La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

mu.Peso <-3300  #valore medio recuperato da fonti online
mu.Lunghezza <- 500 #valore medio recuperato da fonti online

t.test(dati$Peso, mu = mu.Peso)
## 
##  One Sample t-test
## 
## data:  dati$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184
t.test(dati$Lunghezza, mu = mu.Lunghezza)
## 
##  One Sample t-test
## 
## data:  dati$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6628 495.7287
## sample estimates:
## mean of x 
##  494.6958

Conclusione t test per la variabile Lunghezza

  • Il p-value è estremamente basso → rifiutiamo l’ipotesi nulla.
  • La media del campione (494.692) è significativamente diversa da 500.
  • L’intervallo di confidenza non include 500, rafforzando la conclusione

Il risultato suggerisce che il campione non riflette fedelmente la media della popolazione e potrebbe indicare una deviazione sistematica, un bias di selezione o una caratteristica locale del campione.

Conclusione t test per la variabile Peso

  • Il p-value = 0.1296 è maggiore di 0.05 → non rifiutiamo l’ipotesi nulla
  • La media del campione (3284.08) è vicina a 3300 e non significativamente diversa
  • L’intervallo di confidenza include 3300, quindi è plausibile che la media reale sia 3300

Il test t mostra che non ci sono evidenze statisticamente significative per affermare che il peso medio del campione sia diverso da 3300 grammi.

Ipotesi 3.Le misure antropometriche sono significativamente diverse tra i due sessi

boxplot(Peso~Sesso, data= dati)

boxplot(Lunghezza~Sesso,  data= dati)

boxplot(Cranio~Sesso,  data= dati)

variabili <- c("Peso", "Lunghezza", "Cranio")

for (v in variabili) {
  cat("\n---", v, "per Maschi ---\n")
  print(summary(dati[[v]][dati$Sesso == "M"]))
  
  cat("\n---", v, "per Femmine ---\n")
  print(summary(dati[[v]][dati$Sesso == "F"]))
}
## 
## --- Peso per Maschi ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     980    3150    3430    3408    3720    4810 
## 
## --- Peso per Femmine ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2900    3160    3161    3470    4930 
## 
## --- Lunghezza per Maschi ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   320.0   490.0   500.0   499.7   515.0   560.0 
## 
## --- Lunghezza per Femmine ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   490.0   489.8   505.0   565.0 
## 
## --- Cranio per Maschi ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   265.0   334.0   343.0   342.5   352.0   390.0 
## 
## --- Cranio per Femmine ---
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   235.0   330.0   340.0   337.6   348.0   390.0
for (v in variabili) {
  cat("\nAnalisi per:", v, "\n")
    
  formula <- as.formula(paste(v, "~Sesso"))
  
  print(t.test(formula, data = dati))
}
## 
## Analisi per: Peso 
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.115, df = 2488.7, 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:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496 
## 
## 
## Analisi per: Lunghezza 
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750 
## 
## 
## Analisi per: Cranio 
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

Conclusioni: la variabile Sesso risulta essere una variabile di controllo importante nell’analisi in quanto per tutte le misure antropometriche rilevate esiste una significativa differenza di medie fra i due sessi, leggermente meno accentuato nel diametro del cranio ma comunque significativa.

2.2 Creazione del Modello di Regressione

n<- nrow(dati)
install.packages("moments")
## pacchetto 'moments' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Gianpiero\AppData\Local\Temp\RtmpakFkQg\downloaded_packages
library(moments)
skewness(dati$Peso)
## [1] -0.6474036
shapiro.test(dati$Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Peso
## W = 0.97068, p-value < 2.2e-16
plot(density(dati$Peso))

qqnorm(dati$Peso)
qqline(dati$Peso, col = "red", lwd = 2)

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[, !(names(dati) %in% c("Sesso","Fumatrici", "Tipo.parto","Ospedale"))], upper.panel= panel.smooth, lower.panel = panel.cor)

mod.full <- lm(Peso ~ ., data=dati)
summary (mod.full)
## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1123.26  -181.53   -14.45   161.05  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.7960   141.4790 -47.610  < 2e-16 ***
## Anni.madre        0.8018     1.1467   0.699   0.4845    
## N.gravidanze     11.3812     4.6686   2.438   0.0148 *  
## Fumatrici       -30.2741    27.5492  -1.099   0.2719    
## Gestazione       32.5773     3.8208   8.526  < 2e-16 ***
## Lunghezza        10.2922     0.3009  34.207  < 2e-16 ***
## Cranio           10.4722     0.4263  24.567  < 2e-16 ***
## Tipo.partoNat    29.6335    12.0905   2.451   0.0143 *  
## Ospedaleosp2    -11.0912    13.4471  -0.825   0.4096    
## Ospedaleosp3     28.2495    13.5054   2.092   0.0366 *  
## SessoM           77.5723    11.1865   6.934 5.18e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 668.7 on 10 and 2487 DF,  p-value: < 2.2e-16

Il modello con tutte le variabili mostra che alcune di esse risultano non significative per predire il peso neonatale e attraverso il metodo stepwise verranno escluse dal modello utilizzando il metodo BIC. Ad esempio “Anni.madre” con un coefficiente prossimo allo 0 sarà certamente escluso dall’analisi mentre le variabili “Gestazione”, “Lunghezza” e “Cranio” sembrano contribuire significativamente sulla variabile risposta Peso.

Dal modello completo di tutte le variabili emerge che ad ogni settimana di gestazione in più corrisponde un aumento di peso di 32g per il neonato e risulta essere una variabile significativa per l’analisi.

2.3 Selezione del Modello Ottimale

mod.step <- step(mod.full,k = log(n), direction = "both", trace = FALSE)
summary(mod.step)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16
mod.inter1 <- update(mod.step, ~. + Lunghezza*Gestazione)
summary(mod.inter1)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Lunghezza, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1133.41  -179.98   -11.52   168.93  2652.65 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.991e+03  9.206e+02  -2.163 0.030631 *  
## N.gravidanze          1.303e+01  4.321e+00   3.015 0.002594 ** 
## Gestazione           -9.391e+01  2.481e+01  -3.785 0.000157 ***
## Lunghezza            -8.476e-02  2.028e+00  -0.042 0.966661    
## Cranio                1.076e+01  4.264e-01  25.234  < 2e-16 ***
## SessoM                7.225e+01  1.121e+01   6.445 1.38e-10 ***
## Gestazione:Lunghezza  2.729e-01  5.298e-02   5.151 2.79e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7299, Adjusted R-squared:  0.7292 
## F-statistic:  1122 on 6 and 2491 DF,  p-value: < 2.2e-16
mod.inter2 <- update(mod.step, ~. + Lunghezza*Cranio)
summary(mod.inter2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Lunghezza:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.65  -180.93   -13.48   165.99  2865.46 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.803e+03  1.018e+03  -1.771   0.0767 .  
## N.gravidanze      1.293e+01  4.323e+00   2.991   0.0028 ** 
## Gestazione        3.815e+01  3.967e+00   9.616  < 2e-16 ***
## Lunghezza        -3.060e-01  2.203e+00  -0.139   0.8895    
## Cranio           -4.755e+00  3.192e+00  -1.490   0.1365    
## SessoM            7.324e+01  1.120e+01   6.537 7.59e-11 ***
## Lunghezza:Cranio  3.157e-02  6.531e-03   4.835 1.41e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7296, Adjusted R-squared:  0.7289 
## F-statistic:  1120 on 6 and 2491 DF,  p-value: < 2.2e-16
mod.quad <- update(mod.step, ~. +I(Gestazione^2))
summary(mod.quad)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2), data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1144.0  -181.5   -12.9   165.8  2661.9 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4646.7158   898.6322  -5.171 2.52e-07 ***
## N.gravidanze       12.5489     4.3381   2.893  0.00385 ** 
## Gestazione        -81.2309    49.7402  -1.633  0.10257    
## Lunghezza          10.3502     0.3040  34.045  < 2e-16 ***
## Cranio             10.6376     0.4282  24.843  < 2e-16 ***
## SessoM             75.7563    11.2435   6.738 1.99e-11 ***
## I(Gestazione^2)     1.5168     0.6621   2.291  0.02206 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7269 
## F-statistic:  1109 on 6 and 2491 DF,  p-value: < 2.2e-16
mod.inter3 <- update(mod.step, ~. + Lunghezza*Gestazione + Lunghezza*Cranio)
summary(mod.inter3)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Lunghezza + Lunghezza:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1138.65  -180.94   -11.23   168.54  2735.88 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.505e+03  1.028e+03  -1.465  0.14305    
## N.gravidanze          1.304e+01  4.321e+00   3.018  0.00257 ** 
## Gestazione           -5.413e+01  4.488e+01  -1.206  0.22793    
## Lunghezza            -1.080e+00  2.233e+00  -0.484  0.62866    
## Cranio                4.799e+00  5.621e+00   0.854  0.39335    
## SessoM                7.212e+01  1.121e+01   6.434 1.49e-10 ***
## Gestazione:Lunghezza  1.917e-01  9.289e-02   2.064  0.03912 *  
## Lunghezza:Cranio      1.217e-02  1.144e-02   1.064  0.28759    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.3 on 2490 degrees of freedom
## Multiple R-squared:   0.73,  Adjusted R-squared:  0.7293 
## F-statistic: 961.8 on 7 and 2490 DF,  p-value: < 2.2e-16
print(BIC(mod.step, mod.inter1, mod.inter2, mod.inter3,mod.quad))
##            df      BIC
## mod.step    7 35193.65
## mod.inter1  8 35175.01
## mod.inter2  8 35178.14
## mod.inter3  9 35181.69
## mod.quad    8 35196.21
install.packages("car")
## pacchetto 'car' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Gianpiero\AppData\Local\Temp\RtmpakFkQg\downloaded_packages
library(car)
car::vif(mod.step)  #le variabili non presentano collinearità
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023462     1.669779     2.075747     1.624568     1.040184

Dopo aver valutato il modello contenete tutte le variabili, utilizzando la funzione step, con metodo BIC, si è definito il modello migliore. Per selezionare il miglior modello si è considerato quanto emerso dall’analisi grafica condotta in precedenza, si è QUINDI valutata l’eventuale interazione fra LunghezzaVSGestazione e fra LunghezzaVSGrandezza del Cranio ed entrambe contemporaneamente. Successivamente si è valutato se il modello rispondesse meglio alle esigenze di analisi considerando la non linearità della variabile Gestazione. Seppur dai risultati del test BIC si riterrebbe più attendibile il modello “mod.inter1”, considerata la minima differenza di R^2 si decide di continuare l’analisi con il modello mod.step che ha un minore numero di variabili e sembra spiegare bene la variabile risposta considerando i coefficienti delle variabili utilizzate.

summary(mod.step)$r.squared  #percentuale di variabile risposta spiegata dal modello 
## [1] 0.7270152
rmse <- sqrt(mean(residuals(mod.step)^2))
rmse  #basso rispetto alla media della variabile risposta
## [1] 274.3666
mean(dati$Peso)
## [1] 3284.184
rmse/mean(dati$Peso)*100
## [1] 8.354179
par(mfrow=c(2,2))
plot(mod.step)

shapiro.test(residuals(mod.step))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod.step)
## W = 0.97414, p-value < 2.2e-16
plot(density(residuals(mod.step)))

install.packages("lmtest")
## pacchetto 'lmtest' aperto con successo con controllo somme MD5
## 
## I pacchetti binari scaricati sono in
##  C:\Users\Gianpiero\AppData\Local\Temp\RtmpakFkQg\downloaded_packages
library(lmtest)
lmtest::bptest(mod.step)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod.step
## BP = 90.297, df = 5, p-value < 2.2e-16
lmtest::dwtest(mod.step)
## 
##  Durbin-Watson test
## 
## data:  mod.step
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0
lev<-hatvalues(mod.step)
plot(lev)
p=length(coefficients(mod.step)-1)
soglia.lev <- 2 * (p+1) / n
abline(h=soglia.lev, col=2)
lev.alti <- which(lev > soglia.lev)
length(lev.alti)
## [1] 108
plot(rstudent(mod.step))
residui <- rstudent(mod.step)
abline(h=c(-4.5,4.5), col=2)
sum(abs(residui)> 4.5)
## [1] 3
soglia.cook <- 0.5
cook <- cooks.distance(mod.step)
plot(cook, main = "Cook's Distance")
abline(h = 0.5, col = "red", lty = 2)

# Indici con Cook’s Distance alta
cook_alti <- which(cook > soglia.cook)
cook_alti
## 1551 
## 1549
library(car)
outlierTest(mod.step)
##       rstudent unadjusted p-value Bonferroni p
## 1551 10.046230         2.6345e-23   6.5810e-20
## 155   5.025345         5.3818e-07   1.3444e-03
## 1306  4.824963         1.4848e-06   3.7092e-03
max(cook)
## [1] 0.8297645
# Indici con Cook’s Distance alta
cook.alti <- which(cook > soglia.cook)
influenti <- intersect(lev.alti, cook.alti)

influenti
## [1] 1549
plot(cook, type = "h", main = "Cook's Distance")
abline(h = soglia.cook, col = "red", lty = 2)
points(influenti, cook[influenti], col = "blue", pch = 19)


dati[1549, ]
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551         35            1         0         38 4370       315    374
##      Tipo.parto Ospedale Sesso
## 1551        Nat     osp3     F
dati[dati$Lunghezza < dati$Cranio, ]
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551         35            1         0         38 4370       315    374
##      Tipo.parto Ospedale Sesso
## 1551        Nat     osp3     F
mod.nooutlier <- update(mod.step, data = dati[-influenti, ])

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

L’analisi dei residui porta alle seguenti conclusioni:

  1. Un RMSE basso rispetto alla media della variabile risposta indica una buona capacità predittiva del modello.
  2. Il test Shapiro-Wilk segnala una deviazione dalla normalità sulle code, ma il valore di W è alto e il modello può essere considerato statisticamente affidabile
  3. i test di Breusch-Pagan e Durbin-Watson evidenziano eteroschedasticità e assenza di autocorrelazione il che induce a valutare l’ipotesi di utilizzare un modello generalizzato o la trasformazioni delle variabili per maggiore accuratezza delle analisi 4.Dall’analisi degli outliers (3) emerge che uno di essi risulta essere influnete in quanto è anche un valore leverage. Il campione in esame presenta un valore della variabile lunghezza inferiore rispetto alla variabile cranio che porta a dedurre che tale campione è non consistente per l’analisi. Quindi si è ricalcolato il modello escludendo l’outlier.

3. Previsioni e Risultati

df.predict1 <- data.frame(N.gravidanze = 3,
                        Gestazione = 39,
                        Lunghezza = mean(dati$Lunghezza[-1549]),
                        Cranio = mean(dati$Cranio[-1549]),
                        Tipo.parto = 1,
  Ospedale = 1,
  Sesso = "M"
                         )

df.predict2 <- data.frame(N.gravidanze = 3,
                        Gestazione = 39,
                        Lunghezza = mean(dati$Lunghezza[-1549]),
                        Cranio = mean(dati$Cranio[-1549]),
                        Tipo.parto = 1,
  Ospedale = 1,
  Sesso = "F"
                         )
predict(mod.step, newdata = df.predict1, interval = "prediction")
##        fit      lwr      upr
## 1 3349.757 2810.617 3888.897
predict(mod.step, newdata = df.predict2, interval = "prediction")
##        fit      lwr      upr
## 1 3271.776 2732.613 3810.939

Dalle previsioni eseguite risulta che un bambino nato alla 39esita settimana di gestazione da una madre alla terza gravidanza peserà circa 3314g se Maschio e 3236g se Femmina. Tipo di parto e Modalità di parto non hanno un impatto significativo sull’analisi.

4. Visualizzazioni

dati$Fumatrici <- factor(dati$Fumatrici,
                                   levels = c(0, 1),
                                   labels = c("Non fumatrice", "Fumatrice"))

                             
table(dati$Fumatrici)
## 
## Non fumatrice     Fumatrice 
##          2394           104
library("ggplot2")
ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
             y=Peso,
             col=Fumatrici), position = 'jitter')+
  geom_smooth(aes(x=Gestazione,
             y=Peso,
             col=Fumatrici), se=F, method = "lm")+
  labs(title = "Peso vs Gestazione e Fumo Materno")

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")+
  labs(title = "Peso rispetto alle settimane di Gestazione e al Sesso")

ggplot(dati, aes(x = Fumatrici, y = Peso, fill = Fumatrici)) +
  geom_boxplot() +
  labs(title = "Distribuzione del peso per fumo materno",
       x = "Fumatrici",
       y = "Peso del neonato (g)",
       fill = "Fumo materno") +
  scale_fill_manual(values = c("Non fumatrice" = "lightblue", "Fumatrice" = "blue")) +
  theme_minimal()

Di seguito le considerazioni desumibili dall’analisi grafica:

  1. Il primo grafico mostra che con l’aumentare delle settimane di gestazione aumenta il peso del neonato e sembra esserci una tendenza positiva più significativa all’aumento di peso per i bambini nati da madri non fumatrici con l’avanzare delle settimane di gestazione, visibile dalle rette che si incrociano nel grafico.

  2. Dal secondo grafico emerge che il peso del neonato aumenta in manieta identica fra i due sessi con l’avanzare delle settimane di gestazione considerando un peso di poco inferiore del gruppo di neonati di sesso femminile che resta costante all’aumentare del numero di settimane di gestazione.

  3. Dal primo grafico emerge visivamente che il campione analizzato ha una percentuale di neonati con madri non fumatrici nettamente maggiore rispetto ai nati da madri fumatrici. La differenza fra campioni dei due gruppi potrebbe spiegare il motivo per cui la variabile Fumatrici non viene considerata significativa nei vari modelli ma suggerisce che sarebbe utile analizzare l’incidenza del fumo sulla natalità. Anche il peso del neonato alla nascita sembra essere più alto per le madri non fumatrici rispetto alle madri fumatrici.