Caricamento del dataset :

Come primo passo andiamo a caricare il dataset che contiene i dati relativi ai neonati registrati in 3 ospedali diversi.

head(df)
##   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(df)
## '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 ...

1. Analisi delle variabili :

Il nostro dataset è composto da 2500 osservazioni raccolte in 3 ospedali diversi, per costruire il modello predittivo per il peso neonatale sono state raccolte le seguenti variabili :

Il nostro obiettivo è quello di individuare quali di queste variabili sono più predittive del peso del bambino alla nascita, concentrandoci maggiormente sull’impatto che il fumo, e il numero di settimane di gestazione possono avere sul peso del neonato, in quanto ad esempio per un numero di settimane di gestazione basso, si potrebbe avere una nascita prematura.

Analisi e Modellazione :

Analisi preliminare :

Nella prima fase esploreremo le variabile attraverso un’analisi descrittiva, per comprendere la distribuzione delle variabile e identificare eventuali outlier o anomalie

Calcoleremo indici di posizione,variabilità e forma per variabili quantitaive continue, mentre per variabili qualitative e quantitative discrete vedremo la distribuzione di frequenze e l’indice di eterogeneità di Gini.

Cominciamo con l’analisi del peso la nostra variabile target , cominciamo con le prime statistiche di sintesi per poi andare a visualizzare gli indici di variabilità e forma.

# Statistiche descrittive per la variabile peso :
summary(df$Peso)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2990    3300    3284    3620    4930
# Deviazione standard :
sd(df$Peso)
## [1] 525.0387
# Asimmetria :
skewness(df$Peso)
## [1] -0.6470308
# Curtosi :
kurtosis(df$Peso) - 3
## [1] 2.031532
# Visualizziamo graficamente la distribuzione del peso :
plot(density(df$Peso))

Come vediamo sia graficamente, sia grazia all’indice di asimmetria di Fisher di -0.64, la distribuzione del peso risulta essere asimmetrica negativa, presentando una coda della distribuzione più lunga sulla sinistra e con quindi una maggiore concentrazione di valori alti per le nostre osservazioni.

Visto che il peso è la nostra variabile target andiamo ad effettuare un test di Shapiro-Wilk per saggiare l’ipotesi di normalità della distribuzione, facciamo questo preventivamente perché, eventuali allontanamenti dalla normalità delle distribuzione della variabile target, possono influire sui residui del modello stimato, perché essendo un modello lineare esso non riesce bene a filtrare variabili non normali.

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

Con un livello di significatività alfa scelto a priori del 5% , il risultato del test di Shapiro-Wilk torna un p-value di praticamente 0, quindi rifiutiamo l’ipotesi nulla di normalità, la distribuzione del peso dei neonati non segue una distribuzione normale, andremo ad indagare meglio il comportamento dei residui una volta stimato il modello di regressione lineare.

Passiamo adesso alle statistiche descrittive delle altre variabili quantitative continue :

variabili_continue = c("Anni.madre","Gestazione","Lunghezza","Cranio")

for (var in variabili_continue) {
  # Statistiche descrittive
  cat("\nStatistica descrittiva per la variabile:", var, "\n")
  print(summary(df[[var]]))
  
  # Deviazione standard
  cat("\nDeviazione standard per", var, ":", sd(df[[var]]), "\n")
  
  # Asimmetria
  cat("\nAsimmetria per", var, ":", skewness(df[[var]]), "\n")
  
  # Curtosi
  cat("\nCurtosi per", var, ":", kurtosis(df[[var]]) - 3, "\n")
  
  # Visualizzazione dell'istogramma
  p = ggplot(df, aes_string(x = var)) +
    geom_histogram(fill = "blue",color = "black",bins = 20) +
    ggtitle(paste("Istogramma per la variabile:", var)) +
    theme_minimal() +
    xlab(var) +
    ylab("Frequenza")
  
  print(p)
  
}
## 
## Statistica descrittiva per la variabile: Anni.madre 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   25.00   28.00   28.16   32.00   46.00 
## 
## Deviazione standard per Anni.madre : 5.273578 
## 
## Asimmetria per Anni.madre : 0.0428115 
## 
## Curtosi per Anni.madre : 0.3804165
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## 
## Statistica descrittiva per la variabile: Gestazione 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.00   38.00   39.00   38.98   40.00   43.00 
## 
## Deviazione standard per Gestazione : 1.868639 
## 
## Asimmetria per Gestazione : -2.065313 
## 
## Curtosi per Gestazione : 8.25815

## 
## Statistica descrittiva per la variabile: Lunghezza 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   500.0   494.7   510.0   565.0 
## 
## Deviazione standard per Lunghezza : 26.31864 
## 
## Asimmetria per Lunghezza : -1.514699 
## 
## Curtosi per Lunghezza : 6.487174

## 
## Statistica descrittiva per la variabile: Cranio 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     235     330     340     340     350     390 
## 
## Deviazione standard per Cranio : 16.42533 
## 
## Asimmetria per Cranio : -0.7850527 
## 
## Curtosi per Cranio : 2.946206

Analizziamo le variabili una alla volta, partiamo dalla variabile anni madre, vediamo come in media l’età delle madri si aggira intorno ai 28 anni e una mediana di 28, quindi i dati non sono nemmeno troppo variabili considerando anche una deviazione standard abbastanza bassa e di 5.3 circa.

Sia dal grafico che dall’indice di asimmetria di Fisher possiamo stabilire come la distribuzione delle età della madre sia molto simile a quella di una distribuzione normale, anche come indica l’indice di Fisher con un valore prossimo allo 0.

Notiamo subito dal valore del minimo del summary una anomalia, infatti viene rilevata una madre che ha una età di 0 anni, chiaramente un errore di campionamento, vediamo meglio la distribuzione di frequenze per l’età della madre per vedere meglio la situazione :

table(df$Anni.madre)
## 
##   0   1  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30 
##   1   1   1   2   6  13  18  24  45  66  74 100 115 131 180 184 197 172 174 200 
##  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46 
## 147 159 110  96  66  64  41  38  27  19  13   8   2   4   1   1

Come vediamo dalla distribuzione di frequenze, non c’è solo una anomalia, ma bensì 2, infatti abbiamo anche una madre con età 1 anno, cosa alquanto impossibile, quindi adesso andremo a liberarci di tali valori.

df = df[df$Anni.madre >= 13,]
summary(df$Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   25.00   28.00   28.19   32.00   46.00

Perfetto passiamo adesso alla variabile Gestazione, possiamo vedere dal summary una media di 39 circa e una mediana di 38, con una deviazione standard di 1.87 circa , come ben visibile dal grafico la distribuzione risulta essere asimmetrica negativa con un indice di asimmetria di Fisher di -2 circa, che conferma appunto una forte asimmetria negativa.

Per la variabile lunghezza abbiamo una media di 494.7 mm e una mediana di 500, con una deviazione standard di 26.32 circa, dal grafico ma anche dall’indice di asimmetria di Fisher notiamo una forte asimmetria negativa con valore dell’indice di -1.51.

Infine per la variabile cranio abbiamo media di 340, una mediana anche essa di 340, una deviazione standard di 16.42, la distribuzione risulta essere come dal grafico moderatamente asimmetrica negativa , con un indice di asimmetria di Fisher di -0.78.

Vediamo ora le variabili qualitative e quantitative discrete :

Partiamo dalla variabile numero di gravidanze, di essa vediamo gli indici di posizione variabilità e forma,ed anche la distribuzione di frequenze.

summary(df$N.gravidanze)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  1.0000  0.9816  1.0000 12.0000
sd(df$N.gravidanze)
## [1] 1.280949
skewness(df$N.gravidanze)
## [1] 2.513412
kurtosis(df$N.gravidanze)
## [1] 13.98163
table(df$N.gravidanze)
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12 
## 1095  817  340  150   48   21   11    1    8    2    3    1    1

Allora abbiamo una media di 1 gravidanza circa, quindi in media le madri hanno già avuto una gravidanza precedente, una mediana di 1, deviazione standard di 1.28, inoltre la distribuzione della variabile numero di gravidanze è asimmetrica positiva con un indice di asimmetria di Fisher di 2.51.

Vediamo ora la distribuzione di frequenze e l’indici di eterogeneità di Gini per le variabili qualitative, fumo materno, tipo parto, ospedale di nascita e sesso del neonato.

Partiamo dal fumo materno :

table(df$Fumatrici)
## 
##    0    1 
## 2394  104
gini.index = function(x){
  x = sort(x)
  N = length(x)
  ni = table(x)
  fi = ni/N
  fi2 = fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  
  return(gini.normalizzato) 
}

gini.index(df$Fumatrici)
## [1] 0.1595999
# bar chart per visualizzazione grafica :
ggplot(df,aes(x=Fumatrici))+
  geom_bar(fill="lightblue",color="black")+
  labs(title = "Distribuzione delle madri fumatrici",
       x="Madri",
       y="Frequenze")

Come vediamo dalla distribuzione delle frequenze assolute, la maggior parte delle madri non sono fumatrici, la moda infatti risulta essere quella delle madri non fumatrici, dall’indice di eterogeneità di Gini abbiamo appunto come la distribuzione sia per lo più omogenea con un valore dell’indice di 0.16 circa, quindi molto vicino allo 0 che come sappiamo indica massima omogeneità nei dati.

Vediamo ora la variabile tipo parto :

table(df$Tipo.parto)
## 
##  Ces  Nat 
##  728 1770
gini.index = function(x){
  x = sort(x)
  N = length(x)
  ni = table(x)
  fi = ni/N
  fi2 = fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  
  return(gini.normalizzato) 
}

gini.index(df$Tipo.parto)
## [1] 0.8259995
# bar chart per visualizzazione grafica :
ggplot(df,aes(x=Tipo.parto))+
  geom_bar(fill="lightgreen",color="black")+
  labs(title = "Distribuzione del tipo di parto",
       x="Tipo parto",
       y="Frequenze")

Come vediamo dalla tabella delle frequenze assolute la maggior parte dei parti è avvenuta naturalmente con ben 1770 parti naturali contro 728 cesarei, quindi la moda è parti naturali, l’indici di Gini di 0.83 circa indica una buona eterogeneità dei dati.

Analizziamo adesso la variabile ospedale di nascita :

table(df$Ospedale)
## 
## osp1 osp2 osp3 
##  816  848  834
gini.index = function(x){
  x = sort(x)
  N = length(x)
  ni = table(x)
  fi = ni/N
  fi2 = fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  
  return(gini.normalizzato) 
}

gini.index(df$Ospedale)
## [1] 0.9998763
# bar chart per visualizzazione grafica :
ggplot(df,aes(x=Ospedale))+
  geom_bar(fill="lightyellow",color="black")+
  labs(title = "Distribuzione del tipo di ospedale",
       x="Tipo di Ospedale",
       y="Frequenze")

Come vediamo dalla tabella di frequenze abbiamo che l’ospedale dove sono nati la maggior parte di bambini è l’ospedale 2, ma comunque in media le distribuzioni sembrano essere abbastanza uguali, questo è molto visibile dal grafico a barre tra i rispettivi ospedali dove non si vede una netta differenza fra i vari ospedali, inoltre l’indice di Ginii di 0.99 circa sta a significare una forte eterogeneità, tutte le modalità della variabile sono equamente distribuite.

Infine vediamo la variabile sesso del neonato :

table(df$Sesso)
## 
##    F    M 
## 1255 1243
gini.index = function(x){
  x = sort(x)
  N = length(x)
  ni = table(x)
  fi = ni/N
  fi2 = fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  
  return(gini.normalizzato) 
}

gini.index(df$Sesso)
## [1] 0.9999769
# bar chart per visualizzazione grafica :
ggplot(df,aes(x=Sesso))+
  geom_bar(fill="lightpink",color="black")+
  labs(title = "Distribuzione del sesso del neonato",
       x="Sesso del neonato",
       y="Frequenze")

Per la variabile sesso abbiamo un maggior numero di bambine femmine con una frequenza assoluta di 1255, come vediamo dal grafico e dall’indice di Gini con valore di 0.99 , abbiamo una distribuzione equa per le due modalità della variabile.

Matrice di correlazione e test statistici :

Andiamo adesso a vedere come le variabili sono correlate con la nostra variabile oggetto di studio Peso, questo ci darà un’idea su come costruire il nostro modello, in quanto dobbiamo cercare quelle che sono le variabili che riescono maggiormente a spiegare, dare informazioni sul peso del neonato.

Cominciamo con la matrice di correlazione ma considerando solamente le variabili quantitative, per le variabili qualitative valuteremo l’associazione con il peso del neonato in maniera diversa.

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)
}
df_copia <- df[, c("Peso", "Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")]
pairs(df_copia, upper.panel = panel.smooth, lower.panel = panel.cor)

Osservando la matrice di correlazione che mostra le correlazioni a due a due fra le variabili, nel nostro caso anche con il relativo scatterplot, possiamo vedere come le variabili con una correlazioni maggiore con la nostra variabile target Peso sono :

  • Gestazione con un valore di 0.59, come vediamo dallo scatterplot le settimane di gestazione sembrano avere un effetto abbastanza lineare positivo sul peso del neonato, infatti la pendenza della nuvola dei punti sembra essere verso l’alto.

  • Lunghezza con un valore di 0.80 , lo scatterplot mostra infatti una nuvola dei punti che ha una pendenza verso l’alto, con una relazione che sembra essere molto lineare fra le due variabili, quindi all’aumentare della lunghezza del bambino aumenta anche il suo peso

  • Cranio con un valore di 0.70, e come vediamo anche dallo scatterplot della nuvola di punti mostra un effetto lineare positivo sul peso del bambino, quindi maggiori sono i mm del diametro del cranio, più saranno i kg del bambino

Dalla matrice di correlazione è possibile vedere anche quelle che sono le variabili che hanno una bassa correlazione con la variabile target, in particolare le variabili Anni Madre e Numero di gravidanze, presentano una correlazione quasi nulla con la variabile peso, questo è visibile anche dagli scatterplot delle due variabili con la variabile target, dove vediamo come la nuvola dei punti è sparsa orizzontalmente, non presentando alcun pattern particolare, denotando quindi un effetto probabilmente nullo delle variabili sul peso del neonato.

In generale poi i regressori fra loro non sembrano avere valori elevati di correlazioni fra loro, quindi non dovrebbero esserci problemi di multicollinearità, ma ne avremo conferma più avanti dopo aver stimato il modello e analizzati i residui.

Specifichiamo anche che solo dopo aver stimato un primo modello potremmo dire che fermezza se alcune variabili sono significativamente influenti sulla variabile target, questo verrà effettuato analizzando i coefficienti di regressione beta delle varie variabili, in quanto sono essi stessi a stabilire la significatività della relazione con la variabile target.

Per le variabili qualitative invece non ci è possibile usare degli scatterplot in quanto non si prestano bene a cogliere relazioni fra una variabile quantitativa come il peso, e variabili qualitative come ad esempio il sesso, per questo è molto meglio usare dei boxplot e verificare l’associazione tra le variabili tramite test statistici.

Cominciamo mostrando il boxplot fra peso e madri fumatrici :

par(mfrow=c(1,2))
boxplot(df$Peso)
boxplot(df$Peso ~ df$Fumatrici)

Dall’analisi del boxplot vediamo come il fumo ha un impatto negativo sul peso del neonato, andiamo anche a vedere tramite un test t per confronto fra medie di gruppi indipendenti, se l’impatto del fumo è significativo sul peso del bambino.

t.test(df$Peso ~ df$Fumatrici)
## 
##  Welch Two Sample t-test
## 
## data:  df$Peso by df$Fumatrici
## t = 1.0362, df = 114.12, p-value = 0.3023
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.5076 145.3399
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.262        3236.346

Il test t per il confronto fra medie ha dato un p-value di 0.30, che ci porta a non rifiutare l’ipotesi nulla sull’uguaglianza fra medie, quindi l’impatto del fumo sul peso del bambino sembra non essere abbastanza significativo, questo probabilmente è dovuto alla dimensione del campione, visto anche come il numero della madri fumatrici è nettamente inferiori a quello delle non fumatrici, in generale sappiamo invece che il fumo ha un impatto molto negativo sul peso del bambino.

Vediamo ora il boxplot fra il peso e il tipo di parto :

boxplot(df$Peso ~ df$Tipo.parto)

Come è evidente dal boxplot, non c’è nessuna associazione fra il peso e il tipo del parto del neonato, quindi il tipo di parto non ha alcun effetto sul peso del bambino alla nascita, possiamo dire la stessa cosa a priori per il tipo di ospedale, in quanto di certo il peso di un bambino non dipende da quello.

Andiamo invece ad analizzare la relazione fra peso e sesso del neonato tramite boxplot :

boxplot(df$Peso ~ df$Sesso)

Dal boxplot possiamo vedere come il sesso del bambino abbia un effetto sul peso, infatti dal grafico notiamo come ci sia una differenza nella distribuzione del peso per bambini maschi che pesano di più delle femmine. Andiamo a saggiare l’ipotesi di uguaglianza fra medie usando un test t, per verificare la significatività del sesso sul peso dei bambini.

t.test(df$Peso ~ df$Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  df$Peso by df$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

Con un p-value molto vicino allo 0, si rifiuta l’ipotesi nulla di uguaglianza e non si riufiuta l’ipotesi alternativa, le medie del peso per i due sessi sono significativamente differenti, quindi il sesso sarà sicuramente una variabile significativa per determinare il peso del bambino.

Prima di passare alla modellazione di un primo modello di regressione lineare multipla, ci è stato chiesto di verificare tramite test statistici le seguenti ipotesi :

  • In alcuni ospedali si fanno più parti cesarei

  • La media del peso e della lunghezza del nostro campione di neonati sono significativamente uguali a quelle della popolazione

  • Le misure antropometriche sono significativamente diverse fra i due sessi

Per la prima ipotesi dobbiamo in pratica verificare l’associazione fra il tipo di parto e l’ospedale, possiamo fare questo usando il test chi-quadrato che misura l’associazione fra variabili qualitative su scala nominale :

# creaiamo la tabella di contingenza :
tabella_cont = table(df$Tipo.parto,df$Ospedale)

# eseguiamo il test chi-quadrato :
chisq.test(tabella_cont)
## 
##  Pearson's Chi-squared test
## 
## data:  tabella_cont
## X-squared = 1.083, df = 2, p-value = 0.5819

Il test chi-quadrato ha tornato un p-value di 0.58 circa quindi non si riufiuta l’ipotesi nulla, pertanto non c’è una relazione significativa fra il tipo di parto e i vari ospedali.

Andiamo adesso a saggiare l’ipotesi di uguaglianza fra medie del peso e della lunghezza del nostro campione di neonati tramite test t, consideriamo come valori di riferimento i dati presi dall’OSM che indica come peso medio dei neonati un peso di 3400 grammi, e una lunghezza in mm di 500, per verificare che le due medie sono uguali andiamo ad effettuare un test t :

# valori medi di peso e lunghezza della popolazione :
media_peso_pop = 3400
media_lunghezza_pop = 500

# test t per il peso :
t.test(df$Peso,mu = media_peso_pop)
## 
##  One Sample t-test
## 
## data:  df$Peso
## t = -11.021, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 3400
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184
# test t per la lunghezza :
t.test(df$Lunghezza, mu = media_lunghezza_pop)
## 
##  One Sample t-test
## 
## data:  df$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

Con i valori del p-value dei rispettivi test t, possiamo rifiutare l’ipotesi nulla di uguaglianza fra medie, e non rifiutare l’ipotesi alternativa, la media di peso e lunghezza del campione di bambini è significativamente diversa da quella della popolazione.

Come ultimo test dobbiamo andare a verificare che le misure antropometriche quindi lunghezza e cranio, sono significativamente diverse fra i due sessi, useremo come fatto in precedenza anche per il peso un test t dove mettiamo a confronto lunghezza e diametro del cranio con il sesso del bambino :

# test t per lunghezza e sesso :
t.test(df$Lunghezza ~ df$Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  df$Lunghezza by df$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
# test t per diametro cranio e sesso :
t.test(df$Cranio ~ df$Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  df$Cranio by df$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

Dai valori del p-value dei rispettivi test t, i quali sono molto vicini fra loro, rifiutiamo l’ipotesi nulla di uguaglianza fra medie, e non riufitiamo l’ipotesi alternativa di differenza fra medie, le misure antropometriche di lunghezza e diametro del cranio sono significativamente diverse fra loro.

Model8lazione di un modello di regressione :

Creiamo un primo modello di regressione lineare multipla inserendo al suo interno tutte le variabili esplicative a nostro disposizione :

# primo modello di regressione :
mod1 = lm(df$Peso ~ ., data = df)

# summary del modello :
summary(mod1)
## 
## Call:
## lm(formula = df$Peso ~ ., data = df)
## 
## 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

Abbiamo un primo modello che presenta un R quadro aggiustato di 0.73 circa che non è male, ma stiamo appunto considerando un modello abbastanza complesso e che fa uso di tutte le variabili, analizziamo i coefficienti beta per capire meglio quali variabili sono realmente significative all’interno del modello.

Dei coefficienti beta di ogni variabile dobbiamo vedere la stima, quindi il valore in valore assoluto, il segno, e infine il valore del p-value che ci dice quanto quel coefficiente sia significativo, e quindi quanto è significativa la variabile esplicativa sulla variabile risposta.

Analizziamo prima i coefficienti più significativi per le varie variabili :

  • La variabile gestazione con una stima del coefficiente beta di 32.58 circa indica un effetto di incremento marginale sulla variabile peso, quindi per ogni incremento unitario delle settimane di gestazione un bambino peserà in media 32.58 grammi in più, questo è confermato anche dal p-value prossimo allo 0 che denota la forte significatività del coefficiente beta per la gestazione.

  • La variabile lunghezza con una stima del coefficiente beta di 10.3 circa , ha anche essa un effetto positivo marginale sulla variabile peso, all’incremento unitario dei mm di lunghezza del neonato, si avrà in media un aumento del peso di 10.3 grammi, anche qui con una alta significatività.

  • La variabile diametro del cranio con una stima del beta di regressione di 10.47 , presenta un effetto di incremento positivo sulla variabile peso, anche qui con un p-value praticamente di 0 , indicatore di una forte significatività

  • La variabile sesso, che prende in considerazione la baseline di maschio come modalità principale, con una stima del beta di regressione di 77.57 , essa rappresenta l’effetto marginale della variabile sesso tenendo fissate le altre variabili, in quanto variabile di controllo, questo significa che a parità delle altre variabili, il peso dei bambini è in media di 77.57 grammi in più rispetto all’altra modalità della variabile sesso, ovvero per le femmine, questo lo vediamo anche con una forte significatività data da un p-value di quasi 0.

Come vediamo dalla tabella dei coefficienti altre variabili invece risultano essere totalmente non influenti sul peso del bambino, come la variabile Anni madre che con un p-value di 0.48 e quindi maggiore del livello di significatività alfa del 5% è praticamente non significativa, come avevamo visto in precendenza anche dalla matrice di correlazione.

La variabile Numero di Gravidanze nonostante abbia un coefficiente con significatività sotto il 5%, quindi considerabile significativo, sappiamo come questa variabile ha un effetto nullo sul peso del bambino quindi trascurabile.

Stessa cosa vale per il tipo di parto e l’ospedale di nascita, che non incidono minimamente sul peso del neonato.

Per quanto riguarda la variabile Fumatrici essa presenta un beta di regressione di -30 circa, ma con un p-value di 0.27 non rifiutiamo l’ipotesi nulla quindi la variabile non è statisticamente significativa, questo ripetiamo è dovuto dalla grandezza del campione, che presenta un numero di madri fumatrici nettamente inferiore a quello della madri non fumatrici, come è stato evidente notare in precendenza il fumo ha un impatto negativo sul peso, e questo ci viene anche segnalato dal segno del coefficiente beta, infatti potremmo dire che in media il peso di un neonato decresce in media di 30 grammi per madri fumatrici. Questo però non è il caso dato il nostro campione di dati.

Calcoliamo i VIF del modello per vedere se ci sono problemi di multicollinearità fra i vari regressori :

# VIF per modello 1 :
vif(mod1)
##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.190241  1        1.090982
## N.gravidanze 1.189278  1        1.090540
## Fumatrici    1.007426  1        1.003706
## Gestazione   1.695675  1        1.302181
## Lunghezza    2.086879  1        1.444604
## Cranio       1.631049  1        1.277125
## Tipo.parto   1.004227  1        1.002111
## Ospedale     1.004267  2        1.001065
## Sesso        1.040743  1        1.020168

Tutti i valori del VIF sono minori della soglia di 5, quindi non ci sono problemi di multicollinearità fra i regressori.

Selezione del modello migliore :

Dopo aver visto come si comportano i regressori nel primo modello, andiamo a costruire un secondo modello andando a togliere dal modello precedente quelle variabile che sappiamo essere non influenti sulla risposta, andiamo a rimuovere dal modello le variabili : Anni Madre, Tipo parto e Ospedale dal nostro modello.

# update del modello precendente, creazione di un secondo modello di regressione :
mod2 = update(mod1,~. - Anni.madre - Tipo.parto - Ospedale)

# summary modello 2 :
summary(mod2)
## 
## Call:
## lm(formula = df$Peso ~ N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.24  -181.32   -15.73   162.95  2635.69 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6682.2637   135.7983 -49.207  < 2e-16 ***
## N.gravidanze    12.6996     4.3470   2.921  0.00352 ** 
## Fumatrici      -30.5728    27.6048  -1.108  0.26818    
## Gestazione      32.6437     3.8079   8.573  < 2e-16 ***
## Lunghezza       10.2309     0.3011  33.979  < 2e-16 ***
## Cranio          10.5366     0.4265  24.707  < 2e-16 ***
## SessoM          78.1596    11.2117   6.971 4.01e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1106 on 6 and 2491 DF,  p-value: < 2.2e-16

Il nostro nuovo modello presenta un R quadro aggiustato praticamente identico a quello precendente, quindi a parità di varianza spiegata attualmente preferiamo un modello più parsimonioso e il rasoio di Occam, possiamo vedere come i beta di regressioni delle nostre variabili sono rimasti per lo più invariati, abbiamo ancora che il coefficiente della variabile fumatrici risulta ancora statisticamente poco significativo per cui ha un effetto marginale praticamente nullo sul peso.

Vediamo tramite il criterio di informazione Bayesiano quale dei due modelli fra il primo e il secondo risulta migliore :

BIC(mod1,mod2)
##      df      BIC
## mod1 12 35215.45
## mod2  8 35200.24

Con un BIC di 35200 preferiamo il modello 2, andiamo adesso ad usare la funzione stepAIC che usa la procedura stepwise , per andare ad individuare il modello migliore e più semplice, imposteremo la funzione con il parametro k = log(n) che equivale al criterio di informazione Bayesiano BIC che predilige modelli più parsimoniosi :

# procedura stepwise per selezione modello migliore :
mod_step = stepAIC(mod1,direction = "both", k=log(nrow(df)))
## Start:  AIC=28118.61
## df$Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36710 186779904 28111
## - Fumatrici     1     90677 186833870 28112
## - Ospedale      2    687555 187430749 28112
## - N.gravidanze  1    446244 187189438 28117
## - Tipo.parto    1    451073 187194266 28117
## <none>                      186743194 28119
## - Sesso         1   3610705 190353899 28159
## - Gestazione    1   5458852 192202046 28183
## - Cranio        1  45318506 232061700 28654
## - Lunghezza     1  87861708 274604902 29074
## 
## Step:  AIC=28111.28
## df$Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     91599 186871503 28105
## - Ospedale      2    693914 187473818 28105
## - Tipo.parto    1    452049 187231953 28110
## <none>                      186779904 28111
## - N.gravidanze  1    631082 187410986 28112
## + Anni.madre    1     36710 186743194 28119
## - Sesso         1   3617809 190397713 28151
## - Gestazione    1   5424800 192204704 28175
## - Cranio        1  45569477 232349381 28649
## - Lunghezza     1  87852027 274631931 29066
## 
## Step:  AIC=28104.68
## df$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    702925 187574428 28098
## - Tipo.parto    1    444404 187315907 28103
## <none>                      186871503 28105
## - N.gravidanze  1    608136 187479640 28105
## + Fumatrici     1     91599 186779904 28111
## + Anni.madre    1     37633 186833870 28112
## - Sesso         1   3601860 190473363 28144
## - Gestazione    1   5358199 192229702 28168
## - Cranio        1  45613331 232484834 28642
## - Lunghezza     1  88259386 275130889 29063
## 
## Step:  AIC=28098.41
## df$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    467626 188042054 28097
## <none>                      187574428 28098
## - N.gravidanze  1    648873 188223301 28099
## + Ospedale      2    702925 186871503 28105
## + Fumatrici     1    100610 187473818 28105
## + Anni.madre    1     44184 187530244 28106
## - Sesso         1   3644818 191219246 28139
## - Gestazione    1   5457887 193032315 28162
## - Cranio        1  45747094 233321522 28636
## - Lunghezza     1  87955701 275530129 29051
## 
## Step:  AIC=28096.81
## df$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188042054 28097
## - N.gravidanze  1    621053 188663107 28097
## + Tipo.parto    1    467626 187574428 28098
## + Ospedale      2    726146 187315907 28103
## + Fumatrici     1     92548 187949505 28103
## + Anni.madre    1     45366 187996688 28104
## - Sesso         1   3650790 191692844 28137
## - Gestazione    1   5477493 193519547 28161
## - Cranio        1  46098547 234140601 28637
## - Lunghezza     1  87532691 275574744 29044
# summary del modello scelto con stepAIC :
summary(mod_step)
## 
## Call:
## lm(formula = df$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = df)
## 
## 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
# BIC del modello con stepAIC :
BIC(mod_step)
## [1] 35193.65

La funzione stepAIC ha scelto un modello che contiene solamente le variabili, N.gravidanze, Gestazione, Lunghezza , Cranio e Sesso che come avevamo visto anche in precendeza sono quelle più significative, il modello scelto dalla funzione stepAIC è anche molto simile al modello 2 , che però conteneva anche la variabile fumatrici, come possiamo vedere anche dal BIC del modello scelto dalla funzione stepAIC esso presenta anche un BIC minore rispetto al modello 2 creato da noi, quindi è più preferibile.

Scegliamo quindi di tenere il modello preso grazie alla funzione stepAIC come modello finale, che riesce a spiegare un 72.65% di variabilità della risposta, controlliamo nuovamente i VIF , per controllare non ci siano problemi di multicollinearità fra i regressori :

# vif modello stepAIC :
vif(mod_step)
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023462     1.669779     2.075747     1.624568     1.040184

Valori ancora sotto la soglia di 5, quindi siamo tranquilli per quanto riguarda la multicollinearità.

Analisi dei residui :

Affinchè un modello di regressione sia buono, dobbiamo verificare delle assunzioni sui residui ovvero :

  • I residui devono seguire una distribuzione normale con media 0

  • La varianza della distribuzione dei residui deve essere costante (omoschedasticità)

  • I residui non devono essere correlati fra loro

Andiamo prima di tutto a fare una analisi grafica dei residui, che ci servirà anche ad individuare eventuali valori di leva o outlier, o una combinazione dei due.

# grafici dei residui :
par(mfrow=c(2,2))

plot(mod_step)

Dal primo grafico possiamo notare come la nuvola dei punti dei residui sembra essere sparsa casualmente attorno alla media di 0 senza seguire un pattern preciso.

Il secondo grafico mostra i punti seguire circa la bisettrice del grafico, quindi i residui seguono una distribuzione normale.

Il terzo grafico mostra anche una nuvola dei punti abbastanza sparsa casualmente, senza un pattern o trend specifico, quindi i residui hanno varianza costante.

L’ultimo grafico mostra quelli che sono i potenziali valori influenti siano essi leverage, outlier o una combinazione dei due, nel nostro caso individuiamo come potenziale valore molto influente l’osservazione 1551.

Andiamo adesso ad effettuare i vari test statistici per saggiare le ipotesi delle assunzioni sui residui, ovvero il test di Shapiro-Wilk per l’ipotesi di normalità, il test di Breusch-Pagan per l’ipotesi di omoschedasticità (varianza costante), e il test di Durbin-Watson per l’ipotesi di incorrelazione.

# Shapiro test :
shapiro.test(residuals(mod_step))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_step)
## W = 0.97414, p-value < 2.2e-16
# Breusch-Pagan test :
bptest(mod_step)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_step
## BP = 90.297, df = 5, p-value < 2.2e-16
# Durbin-Watson test :
dwtest(mod_step)
## 
##  Durbin-Watson test
## 
## data:  mod_step
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0

Dai risultati dei test statistici effettuati sui residui possiamo dire :

  • Shapiro test : si rifiuta l’ipotesi nulla di normalità, i residui non seguono una distribuzione normale

  • Breusch-Pagan test : si rifiuta l’ipotesi nulla di omoschedasticità (varianza costante), i residui non hanno varianza costante

  • Durbin-Watson test : non si rifiuta l’ipotesi nulla di incorrelazione, i residui non sono correlati fra di loro

Andiamo anche a visualizzare la distribuzione dei residui :

plot(density(residuals(mod_step)))

Come vediamo dal grafico infatti la distribuzione dei residui è asimmetrica positiva presentando una lunga coda destra della distribuzione.

Andiamo adesso a fare una verifica dei valori influenti e a calcolare la distanza di Cook :

# valori leverage :
lev = hatvalues(mod_step)
p = sum(lev)
soglia = 2*p / nrow(df)
plot(lev)
abline(h=soglia,col=2)

Vediamo ora gli outliers :

# valori outliers :
plot(rstudent(mod_step))
abline(h=c(-2,2),col=2)

Facciamo un outlier test per verificare quali valori sono significativamente più outlier :

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

Vengono individuati 3 valori outlier, ovvero le osservazioni 1551, già vista in precedenza quindi era da aspettarselo, l’osservazione 155 e l’osservazione 1306.

Vediamo adesso la distanza di Cook :

cook = cooks.distance(mod_step)
plot(cook)

# massima distanza di Cook :
round(max(cook),2)
## [1] 0.83

Con una distanza di Cook che supera la soglia di 0.5 l’osservazione 1551 che è propria quella che è a limite anche col valore di 1, non sarebbe da considerarsi una vera e propria anomalia in quanto non si trova oltre la soglia di allarme, ma sicuramente è un valore molto influente.

Nel complesso non dovremmo sentirci nella posizione di rifiutare il modello, nonostante esso non rispetti le assunzioni di normalità e varianza costante, anche perchè il modello attuale presenta una buona variabilità spiegata e con variabili significative.

3. Previsioni e risultati

Andiamo adesso a fare delle previsioni usando il modello da noi scelto, ci viene chiesto in particolare di fare una previsione del peso per una neonata, quindi di sesso femminile, considerando una madre alla terza gravidanza e alla 39esima settimana di gestazione, si assume anche che la madre non sia fumatrice in quanto la variabile fumo non è inclusa nel modello.

Oltre i dati che ci sono stati forniti il modello richiede anche la lunghezza e il diametro del cranio del bambino, per la nostra previsione useremo la media di entrambe le variabili :

# creaiamo un dataframe con le caratteristiche richieste :
neonata = data.frame(N.gravidanze = 3, Gestazione = 39, Sesso = "F",Lunghezza = round(mean(df$Lunghezza),2), Cranio = round(mean(df$Cranio)))

# facciamo la previsione :
prev = predict(mod_step, newdata = neonata)
prev
##        1 
## 3270.918

Il peso stimato per la neonata è di 3270 grammi circa.

4. Visualizzazioni :

Andiamo a visualizzare l’impatto del numero di settimane di gestazione e del sesso del neonato sul peso previsto :

ggplot(data=df)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso),position='jitter')+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso),
              se=F, #standardError=False
              method = "lm")+
    labs(title='Impatto delle settimane di gestazione e del sesso sul peso del neonato')+
  theme_classic()
## `geom_smooth()` using formula = 'y ~ x'

Visualizziamo anche l’impatto delle settimane di gestazione e del fumo sul peso del neonato :

ggplot(data=df) +
  geom_point(aes(x=Gestazione, y=Peso, col=as.factor(Fumatrici)), position='jitter') +
  geom_smooth(aes(x=Gestazione, y=Peso, col=as.factor(Fumatrici), group = Fumatrici),
              se=F, method = "lm") +
  labs(title='Impatto delle settimane di gestazione e del fumo sul peso del neonato',
       col='Fumatrici') +
  theme_classic()
## `geom_smooth()` using formula = 'y ~ x'