LIBRERIE

library(moments)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(gridExtra)
library(lmtest)
library(car)
  1. RACCOLTA DEI DATI E STRUTTURE DEL DATASET

1.1 IMPORTAZIONE DATASET

dati <- read.csv("neonati.csv",
                       sep = ",",
                       encoding = "latin1",
                       stringsAsFactors = TRUE)
attach(dati)
dati$Fumatrici <- as.factor(dati$Fumatrici)
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

1.2 ANALISI DEL DATASET

n <- nrow(dati)
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   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ 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 ...

Il dataset è costituito da 2500 osservazioni e 10 variabili, suddivise in:

Iniziamo effettuiamo un summary sulle variabili quantitative, discrete e continue, del nostro dataset

summary(dati[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")])
##    Anni.madre     N.gravidanze       Gestazione         Peso     
##  Min.   : 0.00   Min.   : 0.0000   Min.   :25.00   Min.   : 830  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:38.00   1st Qu.:2990  
##  Median :28.00   Median : 1.0000   Median :39.00   Median :3300  
##  Mean   :28.16   Mean   : 0.9812   Mean   :38.98   Mean   :3284  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:40.00   3rd Qu.:3620  
##  Max.   :46.00   Max.   :12.0000   Max.   :43.00   Max.   :4930  
##    Lunghezza         Cranio   
##  Min.   :310.0   Min.   :235  
##  1st Qu.:480.0   1st Qu.:330  
##  Median :500.0   Median :340  
##  Mean   :494.7   Mean   :340  
##  3rd Qu.:510.0   3rd Qu.:350  
##  Max.   :565.0   Max.   :390

Dal summary appena lanciato emerge un’anomalia riguardante la variabile Anni.madre, ovvero un valore minimo di 0. Prendendo in considerazione le stime dell’Organizzazione Mondiale della Sanità (https://www.ospedalebambinogesu.it/gravidanza-in-minorenni-96849/), secondo le quali partoriscono ogni anno circa 16 milioni di ragazze di età compresa tra 15 e 19 anni, vediamo quante madri di età inferiore ai 15 anni sono presenti nel nostro dataset

righe_selezionate <- dati[dati$Anni.madre < 15 , ]

print(righe_selezionate)
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 138          13            0         0         38 2760       470    325
## 1075         14            1         0         39 3510       490    365
## 1152          1            1         0         41 3250       490    350
## 1380          0            0         0         39 3060       490    330
## 1532         14            0         0         39 3550       500    355
##      Tipo.parto Ospedale Sesso
## 138         Nat     osp2     F
## 1075        Nat     osp2     M
## 1152        Nat     osp2     F
## 1380        Nat     osp3     M
## 1532        Ces     osp1     M

Così facendo abbiamo identificato ben due valori anomali, 0 e 1, corrispondenti alle righe 1380 e 1152. Molto probabilmente siamo in presenza di un errore di inserimento dati, pertanto procederemo alla loro sostituzione mediante imputazione di media

media_anni <- mean(Anni.madre[Anni.madre != 0 & Anni.madre != 1],
                     na.rm = TRUE)


Anni.madre[Anni.madre == 1 | Anni.madre == 0] <- media_anni

Verifichiamo se la sostituzione è avvenuta correttamente sommando i valori della variabile Anni.madre corrispondenti a 0 e 1. Infatti, qualora la correzione non fosse andata a buon fine il risultato ottenuto sarebbe 1, in caso contrario, invece, giacchè tali valori anomali sarebbero assenti, la somma darebbe come risultato 0

sum(Anni.madre == 0 | Anni.madre == 1)
## [1] 0

I valori sono stati correttamente sostituiti, pertanto lanciamo un nuovo summary

summary(Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   25.00   28.00   28.19   32.00   46.00

2 ANALISI INDICI DI SIMMETRIA, VARIABILITà E FORMA PER VARIABILI QUANTITATIVE

Creiamo una funzione per l’analisi descrittiva delle variabili quantitative

funzione_IPVF <- function(x){
  options(scipen=999)
  x_sorted <- sort(x)
  n <- length(x_sorted)
  minimo = round(min(x),2)
  massimo = round(max(x),2)
  media = round(mean(x),2)
  mediana = round(median(x),2)
  dev_s = round(sd(x),2)
  CV =  round(sd(x)/mean(x) * 100, 2)
  simmetria = round(skewness(x), 2)
  curtosi = round(kurtosis(x)-3, 2)
  
  distr_indici <- as.data.frame(cbind(minimo, massimo, media, mediana, dev_s, CV, simmetria, curtosi))
  
  rownames(distr_indici) <- deparse(substitute(x))
  
  return(distr_indici)
}

Calcoliamo gli indici di posizione, variabilità e forma per le variabili quantitative discrete e continue

df_anni.madre <- funzione_IPVF(Anni.madre)
df_gravidanze <- funzione_IPVF(N.gravidanze)
df_gestazione <- funzione_IPVF(Gestazione)
df_peso <- funzione_IPVF(Peso)
df_lunghezza <- funzione_IPVF(Lunghezza)
df_cranio <- funzione_IPVF(Cranio)

final_df <- rbind(df_anni.madre, df_gravidanze, df_gestazione,
                  df_peso, df_lunghezza, df_cranio)

print(final_df)
##              minimo massimo   media mediana  dev_s     CV simmetria curtosi
## Anni.madre       13      46   28.19      28   5.22  18.50      0.15   -0.10
## N.gravidanze      0      12    0.98       1   1.28 130.51      2.51   10.99
## Gestazione       25      43   38.98      39   1.87   4.79     -2.07    8.26
## Peso            830    4930 3284.08    3300 525.04  15.99     -0.65    2.03
## Lunghezza       310     565  494.69     500  26.32   5.32     -1.51    6.49
## Cranio          235     390  340.03     340  16.43   4.83     -0.79    2.95

Quanto alle variabili qualitative, calcoliamo le frequenze assolute

summary(dati[, c("Tipo.parto", "Ospedale", "Sesso", "Fumatrici")])
##  Tipo.parto Ospedale   Sesso    Fumatrici
##  Ces: 728   osp1:816   F:1256   0:2396   
##  Nat:1772   osp2:849   M:1244   1: 104   
##             osp3:835

CONSIDERAZIONI

Simmetria: 0.15 -> Distribuzione quasi normale.

Curtosi: -0.10 -> Distribuzione leggermente più piatta rispetto alla normale.

Conclusione: La distribuzione dell’età materna è ben distribuita e simmetrica.

Simmetria: 2.51 -> Distribuzione asimmetrica a destra (molti valori piccoli, pochi grandi).

Curtosi: 10.99 -> La distribuzione ha molti dati vicini alla media, ma anche outlier significativi.

CV (130.51%) -> Molto variabile rispetto alla media.

Conclusione: La maggior parte delle madri ha poche gravidanze, ma ci sono alcune con molte gravidanze (outlier?).

Simmetria: -0.65 -> Più neonati con peso alto, ma con alcuni pesi bassi.

Curtosi: 2.03 -> Leggermente platicurtica, meno valori estremi.

CV (15.99%) -> Variabilità moderata nel peso dei neonati.

Conclusione: Il peso è abbastanza normale, ma ci sono neonati con peso basso che abbassano la media.

Simmetria: -1.51 -> Più neonati con lunghezza maggiore rispetto alla media.

Curtosi: 6.49 -> La distribuzione ha molti dati vicini alla media, ma anche outlier significativi.

Conclusione: La maggior parte dei neonati ha una lunghezza simile, ma ci sono casi di neonati più piccoli del normale.

-Cranio:

Simmetria: -0.79 -> Leggera asimmetria a sinistra (neonati con cranio più piccolo della media).

Curtosi: 2.95 -> Quasi normale.

Conclusione: La circonferenza cranica segue una distribuzione quasi normale, con alcuni neonati con testa più piccola.

Simmetria: -2.07 -> Distribuzione fortemente asimmetrica a sinistra (alcuni parti molto prematuri).

Curtosi: 8.26 -> La distribuzione ha molti dati vicini alla media, ma anche outlier significativi.

Conclusione: La maggior parte dei parti avviene a termine, ma ci sono casi di prematurità significativi.

3 ANALISI E MODELLIZZAZIONE

3.1 DATI SULLE MISURE ANTROPOMETRICHE

Le misure antropometriche dei neonati sono indicatori chiave della loro salute e del loro sviluppo. Queste misure includono peso, lunghezza (altezza) e circonferenza cranica e vengono confrontate con standard di riferimento internazionali per monitorare la crescita e identificare eventuali problemi di sviluppo.

Gli enti principali che forniscono queste curve di crescita sono:

Misure Antropometriche Principali e Standard WHO:

La WHO (World Health Organization) fornisce curve di crescita standardizzate basate su dati di neonati sani provenienti da diversi paesi.

Media: Tipicamente, il peso alla nascita per i neonati a termine (nati tra 37 e 42 settimane di gestazione) è compreso tra 3.2 – 3.4 kg.

Indicatori di salute: Un peso alla nascita inferiore a 2.500 grammi è classificato come basso peso alla nascita e può indicare rischi per la salute, mentre un peso superiore a 4.000 grammi può essere indicativo di macrosomia, che puòcreare complicazioni durante il parto.

Media: La lunghezza media alla nascita per i neonati a termine è di circa 49.5 – 50 cm.

Indicatori di salute: Una lunghezza significativamente inferiore può essere un indicatore di ritardi nella crescita intrauterina, mentre una lunghezza superiore alla media è generalmente meno problematica se proporzionata al peso.

Media: La circonferenza cranica media alla nascita è di circa 34.5 cm.

Indicatori di salute: Una circonferenza cranica troppo piccola può indicare microcefalia, mentre una circonferenza eccessivamente grande può essere un segno di idrocefalia. Entrambe le condizioni richiedono ulteriori valutazioni mediche.

3.2 VERIFICA IPOTESI

  1. IN ALCUNI OSPEDALI SI FANNO PIÚ PARTI CESARI

Il test di indipendenza Chi-quadrato (X^2) è un test statistico che verifica se due variabili categoriali sono indipendenti, oppure se esiste una relazione tra di esse. Il test di indipendenza Chi-quadrato confronta i valori osservati e attesi in una tabella di contingenza e calcola una statistica test.Si usa quando abbiamo due variabili categoriche e vogliamo verificare se esiste una relazione tra di esse. Nel nostro caso: il tipo di parto (Naturale/Cesareo) è indipendente dall’ospedale di nascita?

Creazione tabella di contingenza per le due variabili categoriche

tabella_contingenza <- table(Ospedale, Tipo.parto)

print(tabella_contingenza)
##         Tipo.parto
## Ospedale Ces Nat
##     osp1 242 574
##     osp2 254 595
##     osp3 232 603

Calcolo somma totale delle osservazioni

n_totale <- margin.table(tabella_contingenza)

print(n_totale)
## [1] 2500

Calcolo dei valori attesi

Ogni valore atteso si calcola come:

(totale_riga𝑖 * totale_colonna𝑗) / totale_osservazioni

attese <- tabella_contingenza

for(i in 1:nrow(tabella_contingenza)){
  for(j in 1:ncol(tabella_contingenza)){
    attese[i, j] <- (margin.table(tabella_contingenza,1)[i] * margin.table(tabella_contingenza,2)[j])/ n_totale 
  }
}

print(attese)
##         Tipo.parto
## Ospedale      Ces      Nat
##     osp1 237.6192 578.3808
##     osp2 247.2288 601.7712
##     osp3 243.1520 591.8480

Calcolo statistica Chi-Quadrato

osservate <- tabella_contingenza

x_quadro <- sum((osservate - attese)^2/attese)

print(x_quadro)
## [1] 1.097202

Analisi grafica:

Simuliamo una distribuzione chi-quadrato con due gradi di libertà e vediamo in che punto della curva di densità ricade la statistica calcolata, visualizzando la sua posizione rispetto ad un valore critico calcolato ad un livello di significatività pari a α = 0.01. Se il nostro X_quadro è maggiore del valore critico, allora rifiutiamo l’ipotesi nulla (nessuna relazione tra tipo di parto e ospedale).

plot(density(rchisq(1000000, 2)), xlim = c(0,20))

abline(v=qchisq(0.99, 2), col = 2)

points(x_quadro, 0, cex = 3, col = 4, pch = 20)

Dal grafico possiamo confermare l’impossibilità di rifiutare l’ipotesi nulla, pertanto possiamo concludere che non vi è alcuna associazione statistica tra le due variabili. Riassumendo, il test verifica se il tipo di parto è distribuito in modo diverso tra i vari ospedali. Se le percentuali di parti cesarei variano significativamente tra gli ospedali, significa che l’ospedale potrebbe influenzare la scelta del tipo di parto. Ciò si ottiene verificando se le differenze tra gli ospedali sono casuali o significative. Se il p-value < 0.05, significa che esiste una relazione significativa tra l’ospedale e il tipo di parto.

Ipotesi:

(Ipotesi nulla): La distribuzione dei parti cesarei è uguale tra i diversi ospedali (nessuna relazione tra ospedale e tipo di parto)

(Ipotesi alternativa): La distribuzione dei parti cesarei varia tra gli ospedali (esiste una relazione tra ospedale e tipo di parto)

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

Il p-value calcolato conferma l’impossibilità di rigetto dell’ipotesi nulla, pertanto non c’è dipendenza tra il tipo di parto effettuato e l’ospedale.

Possiamo renderci conto di questo vedendo la distibuzione dei parti mediante barplot:

ggplot(dati, aes(x = Ospedale, fill = Tipo.parto)) +
  geom_bar( position = "dodge",  color = "black") +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Confronto distribuzioni tipologia parto tra odspedali",
       x = "Ospedale",
       y = "Totali parti effettuati") +
  scale_y_continuous(breaks = seq(0,600,50))+
  theme_minimal()

  1. LA MEDIA DEL PESO E DELLA LUNGHEZZA DI QUESTO CAMPIONE DI NEONATI SONO SOGNIFICATIVAMENTE UGUALI A QUELLE DELLA POPOLAZIONE

Prendendo in considerazione i dati del WHO (Organizzazione Mondiale della Sanità) sulle misure antropometriche, precedentemente illustrate, verifichiamo le seguenti ipotesi:

caso n.1:

t.test(dati$Peso,
       mu = 3300,
       conf.level = 0.95,
       alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  dati$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

caso n.2:

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

CONSIDERAZIONI

Test sul peso:

Non possiamo rifiutare l’ipotesi nulla:

Un valore di t vicino a zero indica che la media campionaria è molto simile alla media ipotetica. In questo caso, -1.516 non è molto lontano da zero, suggerendo che non c’è una forte deviazione dalla media ipotetica di 3300. Inoltre, dato che l’intervallo di confidenza rappresenta il range di valori entro il quale ci aspettiamo che la vera media del campione si trovi con il 95% di probabilità, non abbiamo evidenza sufficiente per affermare che la media del campione sia significativamente diversa da 3300, poiché tale valore è incluso nell’intervallo. Anche il p_value calcolato (= 0.1296) conferma tale impossibilità, giacche maggiore del livello di significatività adottato (> di 0.05).

Test sulla lunghezza:

Non possiamo rigettare l’ipotesi nulla:

La statistica t calcolata (-0.58514) è vicinissima allo 0, pertanto non c’è una forte deviazione dalla media ipotetica di 495. Inoltre la media attesa (495) rientra nell’intervallo di confidenza calcolato e il p_value è maggiore di 0.05.

  1. LE MISURE ANTROPOMETRICHE SONO SIGNIFICATIVAMENTE DIVERSE TRA I DUE SESSI

PESO/SESSO

t.test(Peso ~ Sesso, data = dati)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

LUNGHEZZA/SESSO

t.test(Lunghezza ~ Sesso, data = dati)
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.582, df = 2459.3, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.929470  -7.876273
## sample estimates:
## mean in group F mean in group M 
##        489.7643        499.6672

CRANIO/SESSO

t.test(Cranio ~ Sesso, data = dati)
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 0.0000000000001718
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M 
##        337.6330        342.4486

CONSIDERAZIONI

  1. Test sul peso rapportato al sesso:

In pratica, si respinge l’ipotesi nulla secondo cui non ci sono differenze nei pesi medi tra i gruppi F e M, a favore dell’ipotesi alternativa.

  1. Test sulla lunghezza rapportata al sesso:

Si respinge l’ipotesi nulla secondo cui non ci sono differenze nelle lunghezze medie tra i gruppi F e M, a favore dell’ipotesi alternativa.

  1. Test sulle dimensioni del cranio rapportate al sesso:

Si respinge l’ipotesi nulla secondo cui non ci sono differenze nelle dimensioni medie del cranio tra i gruppi F e M, a favore dell’ipotesi alternativa.

3.3 RAPPRESENTAZIONE GRAFICA

Box1 <- ggplot(data=dati, aes(x = Sesso, y = Peso, fill = Sesso)) +
  geom_boxplot() +
  scale_x_discrete(labels = c("FEMMINE", "MASCHI")) +
  labs(title = "Peso in base al sesso",
       x = "Sesso",
       y = "Peso alla nascita (g)") + 
  scale_y_continuous(breaks = seq(0,5000,500))+
  theme_classic()

Box2 <- ggplot(data=dati, aes(x = Sesso, y = Lunghezza, fill = Sesso)) +
  geom_boxplot() +
  scale_x_discrete(labels = c("FEMMINE", "MASCHI")) +
  labs(title = "Lunghezza in base al sesso",
       x = "Sesso",
       y = "Lunghezza alla nascita (cm)") + 
  scale_y_continuous(breaks = seq(300,600,20))+
  theme_classic()
  
Box3 <- ggplot(data=dati, aes(x = Sesso, y = Cranio, fill = Sesso)) +
  geom_boxplot() +
  scale_x_discrete(labels = c("FEMMINE", "MASCHI")) +
  labs(title = "Cranio in base al sesso",
       x = "Sesso",
       y = "Dimensioni cranio alla nascita (cm)") + 
  scale_y_continuous(breaks = seq(100,400,10))+
  theme_classic()


grid.arrange(Box1, Box2, Box3, ncol = 3)

Dal grafico si evince chiaramente che le misure antropometriche delle femmine sono inferiori a quelle dei maschi (linea della mediana F inferiore a quella dei M). Inoltre, come visto in sede di analisi descrittiva, le distribuzioni sono simmetriche negative, giacchè presentano molti outlier nelle code di sinistra, indice del fatto che si registrano nascite di bambini/e di peso, lunghezza e cranio di molto inferiori rispetto alla media.

4.CREAZIONE DEL MODELLO DI REGRESSIONE

Verifichiamo che la variabile risposta sia approsimativamente normale

moments::skewness(Peso)
## [1] -0.6470308
moments::kurtosis(Peso) - 3
## [1] 2.031532
shapiro.test(Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 0.00000000000000022

Dal test appena effettuato risulta che la variabile risposta non segue una distribuzione normale:

Tuttavia, questo non rappresenta un limite per il nostro modello di regressione, giacchè quello che dobbiamo verificare è se i residui (le differenze tra i valori osservati e quelli predetti) siano normalmente distribuiti. Questo assicura che le inferenze statistiche (come gli intervalli di confidenza e i test di ipotesi) siano valide. Inoltre dobbiamo tener presente che il Teorema Del Limite Centrale afferma che la somma (o la media) di un gran numero di variabili indipendenti e identicamente distribuite tende a seguire una distribuzione normale, indipendentemente dalla distribuzione originale, a patto che sia finita la loro varianza. Pertanto, per campioni grandi (N > 1000), la normalità dei dati di partenza diventa meno cruciale per la validità delle inferenze, perché il TLC garantisce che la distribuzione della media campionaria sia approssimativamente normale. Il t-test e la regressione sono comunque affidabili.

4.1 RAPPRESENTAZIONE GRAFICA DELLA CURVA DI DENSITÀ DEL PESO

ggplot(data = dati, aes(x = Peso)) +
  geom_density(color = "blue") +
  geom_vline(aes(xintercept = mean(Peso)), color = "red",  size = 0.5) +
  labs(title = "Curva di densità del Peso", x = "Peso (g)", y = "Densità di probabilità") +
  theme_classic()

# Q-Q plot
qqnorm(dati$Peso, pch = 20)
qqline(dati$Peso, col = "red")

I grafici confermano la presenza di code più pesanti, nello specifico una coda di sinistra molto più lunga, indice della presenza di valori anomali più bassi rispetto alla media. Tuttavia entrambi suggeriscono che i dati del peso seguono approssimativamente una distribuzione normale, nonostante ci siano delle deviazioni per i valori più piccoli (curva verso il basso nel QQp), che potrebbero richiedere ulteriori analisi o considerazioni. Pertanto possiamo procedere con il nostro modello e, sulla base di quanto appreso, considerarlo comunque affidabile

4.2 MATRICE DI CORRELAZIONE

Creiamo una matrice di correlazione includendo solo le variabili quantitative, giacchè procederemo successivamente con un’analisi grafica di quelle qualitative. La matrice non riporterà solo i coefficenti di correlazione, bensì anche i grafici, poichè sono essenziali per consentirci di identificare eventuali relazioni non lineari (effetti quadratici), dal momento che il coefficente di correlazione lineare rileva solo correlazioni di tipo lineare.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor = 1, ...) {
  par(usr = c(0, 1, 0, 1))
  r <- (cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  text(0.5, 0.5, txt, cex = cex.cor)
}
 
dati_sb <-  subset(dati, select = c(Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio))

pairs(dati_sb, upper.panel = panel.smooth,lower.panel = panel.cor)

Grazie a questa matrice siamo in grado di identificare quali tra le variabili oggetto di studio sono significativamente correlate con la variabile Peso. Nello specifico si tratta delle variabili: Lunghezza, Diametro del Cranio e settimane di Gestazione. Quest’ultima sembrerebbe presentare un effetto quadratico.

Verifichiamo questa ipotesi mediante scatterplot

ggplot(dati, aes(x = Gestazione, y = Peso, color = Sesso)) +
  geom_point() +  
  geom_smooth(method = "loess", color = "red") +  
  theme_minimal()

Dal grafico possiamo notare come l’effetto della variabile Gestazione non sia lineare, giacchè il grafico in questione mostra un andamento simile ad una curva, la quale cresce fino a raggiunge un picco, corrispondente al momento della gestazione nel quale il feto raggiunge il suo peso ottimale, per poi appiattirsi o iniziare a scendere (questo potrebbe rappresentare una fase in cui l’estensione della gestazione oltre il termine ottimale non porta a un ulteriore aumento del peso). In sintesi il peso del neonato aumenta rapidamente fino a un certo punto della gestazione e poi si stabilizza o addirittura diminuisce leggermente se la gestazione si prolunga oltre il termine ottimale.

In realtà, anche la variabile N.gravidanze dovrebbe presentare un effetto quadratico. Infatti, oltre la terza gravidanza non ci sono differenze significative nel peso del neonato (Parity and Perinatal Outcomes, PLOS ONE, 2020). Pertanto se le cose stanno realmente così, dovremmo, anche in questo caso, visualizzare un grafico con un andamento simile ad una curva discendente.

ggplot(dati, aes(x = N.gravidanze, y = Peso, color = Sesso)) +
  geom_point() +  # Punti dello scatterplot
  geom_smooth(method = "loess", color = "red") +  # Regressione lineare
  theme_minimal()

Come ipotizzato, l’aumento si ha fino alla terza gravidanza, momento in cui la curva inizia la sua fase discendente

Per ciò che riguarda le variabili qualitative, verifichiamo se il fatto di essere una gestante Fumatrice o meno sia in qualche modo correlata al peso dei neonati. Escluderemo per il momento il Sesso, dato che verrà considerata nel nostro modello di regressione multipla in quanto variabile di controllo essenziale in studi di carattere medico e anche perchè precedentemente è stato dimostrato che c’è una differenza significativa tra le medie dei pesi dei neonati suddivise in base a questa variabile: i maschi pesano in media di più delle femmine. Inoltre, verranno esluse a priori le variabili Ospedale e Tipo di parto, in quanto è logicamente impossibile che possano avere un qualche tipo di effetto sul peso di un neonato.

Calcoliamo l’indice di correlazione

dati$Fumatrici <- as.numeric(as.factor(dati$Fumatrici))

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

Nonostante il valore sia molto piccolo, sembrerebbe esserci una qualche relazione lineare inversa tra le due variabili.

Verifichiamo graficamente quanto visto.

dati$Fumatrici <- as.factor(dati$Fumatrici)

ggplot(dati, aes(x = Fumatrici, y = Peso, fill = Fumatrici)) +
  geom_boxplot() +
  labs(title = "Boxplot del Peso per Fumatrici",
       x = "Fumatrici",
       y = "Peso") +
  theme_classic() + 
  theme(legend.position = "none") +
  scale_fill_manual(values = c("1" = "blue", "2" = "red"))+
  scale_x_discrete(labels = c("NON FUMATRICE", "FUMATRICE"))

Dal grafico emerge una piccola differenza di peso tra neonati partoriti da madri non fumatrici e fumatrici. Sembrerebbe che quest’ultime partoriscano figli di peso minore. I valori delle non fumatrici presentano una variabilità maggiore, xchè il loro numero è più ampio rispetto alle non fumatrici. Tuttavia la differenza è minima.

Verifichiamo se si tratta di una differenza significativa mediante T-student test:

t.test(Peso ~ Fumatrici)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

Il test esclude la possibilità di rigettare l’ipotesi nulla (non ci sono differenze significative nelle due medie di peso) perchè:

Molto probabilmente ciò dipenderebbe dal numero molto piccolo di madri fumatrici presenti nel campione (circa il 4%), come illustrato dal barplot riportato qui di seguito. Difatti la letteratura medica a riguardo (https://www.politicheantidroga.gov.it/media/1678/271_nicotina_prenatale.pdf) conferma che l’assimilazione di nicotina durante le settimane di gestazione ha un effetto negativo sullo sviluppo fisico del bambino. Pertanto non è da escludere che un campione con una frequenza maggiore di madri fumatrici ci consentirebbe di cogliere un effetto più significativo.

4.3 BARPLOT FUMATRICI

ggplot(dati, aes(x = Fumatrici, fill = Fumatrici)) +
  geom_bar(color = "black", width = 0.3) + 
  scale_fill_manual(values = c("2" = "red", "1" = "blue")) + 
  scale_x_discrete(labels = c("NON FUMATRICE", "FUMATRICE")) +
  labs(title = "Numero partorienti non fumatrici e fumatrici",
       x = "Partorienti",
       y = "Frequenze assolute") +
  scale_y_continuous(breaks = seq(0, 2500, 100)) +
  theme_minimal() +
  theme(legend.position = "none")

Anche qui possiamo verificare come il numero delle gestanti fumatrici sia di molto inferiore a quello delle non fumatrici

4.4 REGRESSIONE MULTIPLA

Iniziamo con un modello che includa tutte le variabili che potrebbero avere una certa significatività

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

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

Verifichiamo se sussistono problemi di multicollinearità tra le variabili, questa, infatti, potrebbe verificarsi nel caso in cui due o più variabili indipendenti sono altamente correlate tra loro.

Conseguenze:

Come rilevarla:

VIF (Variance Inflation Factor): un valore > 5 o 10 indica un problema di multicollinearità.

vif(mod1)
##    Lunghezza       Cranio   Gestazione    Fumatrici        Sesso   Anni.madre 
##     2.078644     1.628748     1.694517     1.006659     1.040359     1.186683 
## N.gravidanze 
##     1.184706

Abbiamo ottenuto un buon valore per quanto riguarda R² aggiustato (0.7264), e non è stato rilevato alcun problema di multicollinearità tra variabili, tuttavia il modello presenta al suo interno delle variabili con p-value > 0.05: Fumatrici e Anni.madre. Proviamo un nuovo modello escludendo per motivi statistici entrambe le variabili. Infatti, sebbene sia dimostrato che il fumo abbia un impatto negativo sul peso del feto, tuttavia le analisi condotte sul nostro campione ci spingono ad affermare che non c’è evidenza sufficiente per affermare che questa variabile abbia un effetto significativo sul peso del neonato in questo modello specifico. Pertanto, sulla base del Principio di Parsimonia (rasoio di Occam) manteniamo solo quelle variabili che aggiungono valore predittivo al modello.

mod2 <- update(mod1, ~. - Anni.madre - Fumatrici)

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

Il modello non ha subito variazioni significative, pertanto le variabile escluse non erano statisticamente rilevanti. Difatti, la sua bontà è rimasta alta (Adjusted R² = 0.7265) e la variabile N.gravidanze ha subito una riduzione del suo p-value e quindi un aumento di significatività, pertanto possiamo adottare questo modello come base per le nostre future predizioni.

Verifichiamo la presenza o meno di multicollinearità

vif(mod2)
##    Lunghezza       Cranio   Gestazione        Sesso N.gravidanze 
##     2.074689     1.624465     1.669189     1.040054     1.023475

Multicollinearità assente

Tuttavia, prima di procedere è bene integrare il nostro modello con qualche effetto di interazione tra variabili e con gli effetti quadratici precedentemente rilevati. Iniziamo con un modello che includa le seguenti interazioni ed effetti quadratici:

mod3 <- update(mod2, ~.  + Gestazione*Anni.madre + I(Gestazione^2) + Anni.madre*N.gravidanze)

summary(mod3)
## 
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio + Gestazione + Sesso + 
##     N.gravidanze + Anni.madre + I(Gestazione^2) + Gestazione:Anni.madre + 
##     N.gravidanze:Anni.madre, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1142.02  -181.97   -11.04   165.36  2658.11 
## 
## Coefficients:
##                           Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)             -2535.0039  1140.4888  -2.223              0.02632 *  
## Lunghezza                  10.3246     0.3038  33.990 < 0.0000000000000002 ***
## Cranio                     10.6368     0.4281  24.846 < 0.0000000000000002 ***
## Gestazione               -144.0977    53.7318  -2.682              0.00737 ** 
## SessoM                     75.6104    11.2197   6.739      0.0000000000197 ***
## N.gravidanze                1.8961    27.0294   0.070              0.94408    
## Anni.madre                -62.7658    21.2063  -2.960              0.00311 ** 
## I(Gestazione^2)             1.7306     0.6644   2.605              0.00925 ** 
## Gestazione:Anni.madre       1.6376     0.5427   3.018              0.00257 ** 
## N.gravidanze:Anni.madre     0.2888     0.8411   0.343              0.73139    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2490 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7277 
## F-statistic:   743 on 9 and 2490 DF,  p-value: < 0.00000000000000022

Il modello ottenuto possiende una bontà alta (Adjusted R² = 0.7277), e le interazioni introdotte, eccetto Anni.madre*N.gravidanze, hanno una buona significatività statistica. Nello specifico:

Proviamo ad aggiungere il seguente effetto quadratico

mod4 <- update(mod3, ~.  + I(N.gravidanze^2))

summary(mod4)
## 
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio + Gestazione + Sesso + 
##     N.gravidanze + Anni.madre + I(Gestazione^2) + I(N.gravidanze^2) + 
##     Gestazione:Anni.madre + N.gravidanze:Anni.madre, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1157.24  -182.77   -12.83   164.80  2655.33 
## 
## Coefficients:
##                           Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)             -2277.6098  1144.7706  -1.990              0.04675 *  
## Lunghezza                  10.3037     0.3036  33.937 < 0.0000000000000002 ***
## Cranio                     10.6074     0.4279  24.789 < 0.0000000000000002 ***
## Gestazione               -154.0410    53.8518  -2.860              0.00427 ** 
## SessoM                     75.6246    11.2096   6.746      0.0000000000188 ***
## N.gravidanze              -10.2752    27.5022  -0.374              0.70872    
## Anni.madre                -66.6449    21.2521  -3.136              0.00173 ** 
## I(Gestazione^2)             1.8463     0.6656   2.774              0.00558 ** 
## I(N.gravidanze^2)          -3.4237     1.4641  -2.338              0.01944 *  
## Gestazione:Anni.madre       1.7063     0.5430   3.142              0.00170 ** 
## N.gravidanze:Anni.madre     1.2563     0.9367   1.341              0.17997    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.7 on 2489 degrees of freedom
## Multiple R-squared:  0.7293, Adjusted R-squared:  0.7282 
## F-statistic: 670.4 on 10 and 2489 DF,  p-value: < 0.00000000000000022

Il modello ha una buona capacità predittiva con R² = 0.7268. L’interazione tra Fumatrici e Gestazione non ha alcuna significatività statistica, mentre il coefficente positivo di N.gravidanze e quello negativo della medesima variabile al quadrato, indicano che il peso dei neonati tende a diminuire dopo un certo numero di gravidanze.

Verifichiamo la presenza di multicollinearità

vif(mod4)
##               Lunghezza                  Cranio              Gestazione 
##                2.129315                1.647404              337.700835 
##                   Sesso            N.gravidanze              Anni.madre 
##                1.048016               41.364998              418.886021 
##         I(Gestazione^2)       I(N.gravidanze^2)   Gestazione:Anni.madre 
##              284.389198                4.428448              412.518565 
## N.gravidanze:Anni.madre 
##               53.957046

I valori sono molto alti, tuttavia, non dobbiamo preoccuparci, giacchè quando ci sono interazioni o termini quadratici, il calcolo standard del VIF potrebbe non essere accurato.

Vediamo se rimuovendo le variabili poco significative assistiamo ad un miglioramento del modello

mod5 <- update(mod4, ~. - Anni.madre*N.gravidanze)

summary(mod5)
## 
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio + Gestazione + Sesso + 
##     I(Gestazione^2) + I(N.gravidanze^2) + Gestazione:Anni.madre, 
##     data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1158.68  -182.03   -13.32   164.73  2645.90 
## 
## Coefficients:
##                          Estimate  Std. Error t value             Pr(>|t|)    
## (Intercept)           -4626.66319   900.21959  -5.139       0.000000296863 ***
## Lunghezza                10.33391     0.30410  33.982 < 0.0000000000000002 ***
## Cranio                   10.67764     0.42783  24.957 < 0.0000000000000002 ***
## Gestazione              -83.37827    49.85389  -1.672               0.0946 .  
## SessoM                   76.12662    11.24404   6.770       0.000000000016 ***
## I(Gestazione^2)           1.52546     0.66331   2.300               0.0215 *  
## I(N.gravidanze^2)         0.61456     0.72295   0.850               0.3954    
## Gestazione:Anni.madre     0.05048     0.02802   1.801               0.0717 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic:   949 on 7 and 2492 DF,  p-value: < 0.00000000000000022

Il modello ottenuto ha buona capacità predittiva con R² = 0.0.7264, tuttavia sia N.gravidanze^2, sia Gestazione, sia Gestazione*Anni.madre hanno perso significatività, pertanto non possiamo prenderlo in considerazione.

Sebbene il mod4 possieda buona capacità predittiva e sembrerebbe adeguarsi maggiormente a ciò che la letteratura scientifica ci dice sullo sviluppo fisiologico del feto, ciononostante è un modello troppo complesso che potrebbe presentare ottime performance sui dati usati per l’addestramento ma una scarsa capacità predittiva su nuovi dati e quindi una generalizzazione ridotta.

Dunque, dobbiamo concludere affermando che il miglior modello ottenuto è il mod2

Confrontiamo la procedura con la quale abbiamo selezionato il nostro modello con quella implementata in R, ovvero stepAIC, che è parte del pacchetto MASS. Questa funzione viene utilizzata per eseguire una selezione di modelli statistici, comunemente conosciuta come selezione stepwise, basata sul criterio di informazione di Akaike (AIC). Inseriamo un modello iniziale che tenga conto di tutte le variabili del dataset adottando come criterio di selezione il BIC (Bayesian Information Criterion), che penalizza maggiormente i modelli complessi rispetto all’AIC.

stepwise_mod <- MASS::stepAIC(mod1,
                direction = "both",
                k = log(n)) 
## Start:  AIC=28131.29
## Peso ~ Lunghezza + Cranio + Gestazione + Fumatrici + Sesso + 
##     Anni.madre + N.gravidanze
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     53803 187973654 28124
## - Fumatrici     1     90879 188010731 28125
## - N.gravidanze  1    439797 188359648 28129
## <none>                      187919851 28131
## - Sesso         1   3662833 191582684 28172
## - Gestazione    1   5585184 193505035 28197
## - Cranio        1  45791026 233710877 28669
## - Lunghezza     1  87220887 275140738 29077
## 
## Step:  AIC=28124.18
## Peso ~ Lunghezza + Cranio + Gestazione + Fumatrici + Sesso + 
##     N.gravidanze
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     91892 188065546 28118
## <none>                      187973654 28124
## - N.gravidanze  1    646039 188619694 28125
## + Anni.madre    1     53803 187919851 28131
## - Sesso         1   3671289 191644943 28165
## - Gestazione    1   5531705 193505359 28189
## - Cranio        1  46066755 234040410 28664
## - Lunghezza     1  87218857 275192512 29069
## 
## Step:  AIC=28117.58
## Peso ~ Lunghezza + Cranio + Gestazione + Sesso + N.gravidanze
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188065546 28118
## - N.gravidanze  1    623141 188688687 28118
## + Fumatrici     1     91892 187973654 28124
## + Anni.madre    1     54816 188010731 28125
## - Sesso         1   3655292 191720838 28158
## - Gestazione    1   5464853 193530399 28181
## - Cranio        1  46108583 234174130 28658
## - Lunghezza     1  87632762 275698308 29066
summary(stepwise_mod)
## 
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio + Gestazione + Sesso + 
##     N.gravidanze, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226 < 0.0000000000000002 ***
## Lunghezza       10.2486     0.3006  34.090 < 0.0000000000000002 ***
## Cranio          10.5402     0.4262  24.728 < 0.0000000000000002 ***
## Gestazione      32.3321     3.7980   8.513 < 0.0000000000000002 ***
## SessoM          77.9927    11.2021   6.962     0.00000000000426 ***
## N.gravidanze    12.4750     4.3396   2.875              0.00408 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 0.00000000000000022

Il modello selezionato da R corrisponde al mod2 da noi precedentemente identificato come migliore. Per quanto riguarda la variabile Fumatrici vale quanto detto in precedenza, ovvero che non possiamo adottarla nel nostro modello in base alle analisi condotte sul nostro campione di riferimento, tenendo presente che la letteratura scientifica concorda nel ritenere che il fumo abbia realmente un impatto negativo sul peso del feto (https://www.politicheantidroga.gov.it/media/1678/271_nicotina_prenatale.pdf),effetto che, molto probabilmente, non risulta in questa sede a causa del numero molto piccolo di gestanti fumatrici presenti nel nostro campione. Quanto al numero di gravidanze il modello ci dice che ogni gravidanza precedente è associata a un aumento di peso di 12.48 g, tuttavia si tratta di una semplificazione eccessiva. Infatti, è dimostrato (The Effect of Parity on Birth Weight, American Journal of Obstetrics and Gynecology, 2014) che bammbini nati da madri alla prima gravidanza (primipare) tendono ad avere un peso alla nascita leggermente inferiore rispetto ai neonati nati da madri con più gravidanze (multipare).

Il motivo è duplice:

Tuttavia, è altresì dimostrato che:

Ripetiamo, il modello 4 sembrerebbe essere quello che maggiormente si adatta alla realtà dei fatti, tuttavia il Pricipio di Parsimonia ci induce ad escluderlo a causa della sua complessità

5 ANALISI DELLA QUALITà DEL MODELLO

Confrontiamo tutti i modelli elaborati mediante L’AIC (Akaike Information Criterion) e il BIC (Bayesian Information Criterion).

Tali criteri servono per decidere quale modello è migliore, soprattutto quando si hanno più modelli che spiegano lo stesso fenomeno. Nello specifico, entrambi misurano un compromesso tra la bontà di adattamento (quanto bene il modello si adatta ai dati) e la complessità del modello

# Calcola i BIC per tutti i modelli
aic_values <- AIC(mod1, mod2, mod3, mod4, mod5)

# Ordina i risultati in base al BIC (dal più piccolo al più grande)
aic_sorted <- aic_values[order(aic_values$AIC), ]
print(aic_sorted)
##      df      AIC
## mod4 12 35168.61
## mod3 11 35172.10
## mod2  7 35179.33
## mod1  9 35181.39
## mod5  9 35181.48
# Calcola i BIC per tutti i modelli
bic_values <- BIC(mod1, mod2, mod3, mod4, mod5)

# Ordina i risultati in base al BIC (dal più piccolo al più grande)
bic_sorted <- bic_values[order(bic_values$BIC), ]
print(bic_sorted)
##      df      BIC
## mod2  7 35220.10
## mod1  9 35233.81
## mod5  9 35233.90
## mod3 11 35236.16
## mod4 12 35238.50

L’AIC e il BIC confermano quanto detto finora: il primo predilige il mod4 che è quello più complesso e che sembrerebbe maggiormente adattarsi alla realtà dei fatti, tuttavia è troppo complesso e rischia di incappare in problemi di overfitting, al contrario il BIC predilige il mod2 il quale è molto più semplice.

Eseguiamo un test ANOVA per confrontare i due modelli di regressione lineare al fine di valutare se l’aggiunta delle variabili nel mod4 migliora significativamente la spiegazione della variabilità del peso del neonato.

anova(mod2, mod4)
## Analysis of Variance Table
## 
## Model 1: Peso ~ Lunghezza + Cranio + Gestazione + Sesso + N.gravidanze
## Model 2: Peso ~ Lunghezza + Cranio + Gestazione + Sesso + N.gravidanze + 
##     Anni.madre + I(Gestazione^2) + I(N.gravidanze^2) + Gestazione:Anni.madre + 
##     N.gravidanze:Anni.madre
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2494 188065546                                  
## 2   2489 186513240  5   1552306 4.1431 0.0009464 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RSS (Residual Sum of Squares): misura la variabilità residua non spiegata dal modello.

Quando si aggiungono nuove variabili al modello, il RSS dovrebbe diminuire perché il modello spiega meglio i dati.

La differenza tra i due RSS è il “Sum of Squares (Differenza)”, ovvero quanto in più il modello complesso riesce a spiegare rispetto al modello semplice.

Il modello complesso ha ridotto la variabilità residua di 1552306 unità rispetto al modello semplice.

La statistica F viene calcolata proprio utilizzando questa differenza

Questo valore di F viene poi utilizzato per calcolare il p-value, che ci dice se la riduzione della variabilità è significativa.

In sintesi:

Il valore 1552306 indica quanta variabilità aggiuntiva è spiegata dal modello più complesso. È significativo? Sì, perché il p-value = 0.0009464, molto inferiore a 0.05. Pertanto l’aggiunta delle nuove variabili ha effettivamente migliorato la capacità del modello di spiegare il peso del neonato.

Tuttavia il rischio di overfitting resta alto, pertanto, sempre sulla base del principio di Parsimonia, verrà privilegiato il mod2

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

MOD2:

5.1 ANALISI RESIDUI MOD2

L’analisi dei residui è una fase fondamentale per valutare la qualità e l’affidabilità di un modello di regressione. I residui rappresentano la differenza tra i valori osservati e quelli predetti dal modello

Un modello di regressione lineare si basa su alcune ipotesi fondamentali. L’analisi dei residui serve proprio per controllare se queste ipotesi sono rispettate.

Linearità:

L’associazione tra le variabili indipendenti e la variabile dipendente deve essere lineare. Come verificarlo: Un grafico dei residui vs valori predetti dovrebbe mostrare un andamento casuale senza schemi evidenti.

Omoscedasticità (Varianza Costante):

I residui devono avere una varianza costante (non devono “allargarsi” o “restringersi” al variare dei valori predetti). Come verificarlo: Un funnel o un ventaglio nei residui indica un problema di eteroscedasticità.

Indipendenza dei residui:

I residui devono essere indipendenti tra loro, soprattutto in modelli con dati temporali o spaziali. Test: Il test di Durbin-Watson viene spesso utilizzato per controllare l’autocorrelazione.

Normalità dei residui:

I residui dovrebbero seguire una distribuzione normale, soprattutto per la validità dei test statistici sul modello. Come verificarlo: Un Q-Q plot o il test di Shapiro-Wilk può aiutare a verificarlo.

Inoltre, l’analisi sui residui serve ad identificare:

Outlier: Valori anomali che si discostano molto dal resto dei dati, che hanno un impatto sproporzionato sui coefficienti del modello. Cosa usare: Il grafico dei residui standardizzati o la misura del leverage e del Cook’s distance

Iniziamo effettuando un Breusch-Pagan test, un Durbin-Watson test e uno Shapiro-Wilk normality test per verificare rispettivamente l’Omoscedasticità, l’Indipendenza e la Normalità dei residui

bptest(mod2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 90.253, df = 5, p-value < 0.00000000000000022
dwtest(mod2) 
## 
##  Durbin-Watson test
## 
## data:  mod2
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(residuals(mod2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod2)
## W = 0.97408, p-value < 0.00000000000000022

Risultato:

Breusch-Pagan: p_value < 2.2e-16 -> Eteroschedasticità presente

Durbin-Watson: p_value = 0.1224 -> Nessuna autocorrelazione

Shapiro-Wilk: p_value < 2.2e-16 -> Non normalità dei residui

Sembrerebbe che dati i risultati appena ottenuti, non siamo nella condizione di ritenere il mod2 un modello affidabile. Tuttavia è necessaria un’analisi grafica ulteriore giacchè si tratta di test molto rigidi e molto sensibili.

Infatti, bisogna fare le seguenti considerazioni:

Pertanto procediamo all’analisi grafica

par(mfrow = c(2,2))
plot(mod2, pch = 20)

Come volevasi dimostrare, l’analisi grafica conferma quanto detto:

Analizziamo mendiante una curva di densità la distribuzione dei residui

plot(density(residuals(mod2)))

Simmetria della Curva:

La curva sembra abbastanza simmetrica attorno allo 0, il che è positivo e indica che non ci sono forti bias sistematici nelle previsioni del modello.

Forma della Distribuzione:

La curva è simile a una forma a campana (gaussiana), suggerendo una certa normalità. Tuttavia, sembra leggermente più appuntita e con code leggere rispetto a una distribuzione normale perfetta (leptocurtica).

Code della Distribuzione:

Ci sono code leggere verso destra, che indicano la presenza di pochi residui molto grandi (outlier positivi). Ciò è coerente con i risultati precedenti dove abbiamo osservato alcuni outlier nei grafici dei residui. In particolare c’è un valore nella coda di destra molto grande che corrisponde all’osservazione 1551

La presenza di code leggere e la forma appuntita indicano una deviazione dalla normalità perfetta, tuttavia, la leggera non normalità non dovrebbe essere un problema serio dato il grande numero di osservazioni (N = 2500), perché il teorema del limite centrale aiuta a mantenere valide le stime.

5.2 ANALISI GRAFICA VALORI LEVERAGE

Procediamo ad un’analisi grafica dei valori leverage, essenziali nell’analisi di regressione per identificare quanto una particolare osservazione influisce sulla stima dei parametri del modello di regressione. Ricordiamo che il leverage di un’osservazione è una misura di quanto un punto dati si discosti dalle altre osservazioni in termini di valori delle variabili indipendenti.

Alcuni punti chiave sui valori leverage:

Intervallo di valori: I valori leverage variano tra 0 e 1. Un valore leverage vicino a 1 indica che il punto ha un’influenza elevata sul modello, mentre un valore vicino a 0 indica che il punto ha poca influenza.

Influenza sui residui: Un punto con un leverage elevato può avere un grande impatto sui residui del modello di regressione. Ciò significa che la rimozione o la modifica di tale punto potrebbe cambiare significativamente i risultati del modello.

Identificazione di outlier: Anche se un valore leverage elevato non implica necessariamente che un’osservazione sia un outlier, può indicare che l’osservazione è un “outlier in X”, cioè che è distante dagli altri punti nel dominio delle variabili indipendenti.

# valori leverage

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

print(length(lev[lev>soglia]))
## [1] 152

Stampiamoci i valori leverage

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

Distribuzione Generale dei Leverage:

La maggior parte dei punti ha valori di leverage molto bassi, vicino a 0, il che è normale. Significa che la maggior parte delle osservazioni non ha un’influenza eccessiva sulla regressione.

Punti sopra la Linea Rossa:

Ci sono diversi punti con leverage più alto rispetto alla soglia critica. Questi punti hanno un impatto significativo sulla regressione, perché sono molto distanti dalla media delle variabili indipendenti.

Punto Estremo con Leverage Alto (> 0.05):

Un’osservazione ha un valore di leverage nettamente superiore a tutti gli altri. Questo punto è potenzialmente influenzale e potrebbe alterare la stima dei coefficienti.

L’obiettivo adesso è capire se i punti con leverage alto hanno anche residui molto grandi. Se un’osservazione ha entrambe le caratteristiche, allora è un punto influenzale critico, che potrebbe distorcere il modello di regressione, altrimenti, se il leverage è alto ma il residuo è basso, il punto ha influenza ma non distorce il modello.

Dobbiamo verificare:

#outliers

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

Questo grafico mostra i residui studentizzati del modello mod2, un’analisi utile per individuare outlier e punti influenzali.

Linee Rosse a ±2 (Soglia comune per identificare outlier nei residui):

Residui tra -2 e +2: Normali

Residui oltre ±2: Possibili outlier

Residui oltre ±3: Outlier critici

Distribuzione generale dei residui:

La maggior parte dei punti è compresa tra -2 e +2, il che è buono. Indica che il modello spiega abbastanza bene i dati senza deviazioni eccessive.

Outlier con residui elevati:

Alcuni punti superano +3 o -3, segnalando outlier critici. Il punto più evidente è oltre 10 (potenziale valore anomalo o errore nei dati) e corrisponde all’osservazione 1551.

Distribuzione simmetrica attorno allo zero:

I residui sono bilanciati tra valori negativi e positivi (assenza di bias sistematico). Se ci fosse una concentrazione di residui positivi o negativi, potrebbe indicare problemi nel modello.

Identifichiamo quali tra questi residui studentizzati siano outlier significativi mediante il comando outlierTest()

outlierTest(mod2)
##       rstudent            unadjusted p-value               Bonferroni p
## 1551 10.051908 0.000000000000000000000024906 0.000000000000000000062265
## 155   5.027798 0.000000531379999999999984382 0.001328500000000000036762
## 1306  4.827238 0.000001468099999999999970756 0.003670200000000000198352

Residuo Studentizzato:

Indica quanto un punto si discosta dalla previsione del modello, normalizzando l’errore.

Regola generale:

1551:

Residuo = 10.05 -> Outlier critico p-value Bonferroni = 6.23e-20 -> Molto significativo anche dopo la correzione. Potrebbe essere un errore nei dati o un caso molto anomalo.

155:

Residuo = 5.02 -> Outlier significativo p-value Bonferroni = 1.33e-03 -> Ancora significativo dopo Bonferroni.

1306:

Residuo = 4.82 -> Outlier significativo p-value Bonferroni = 3.67e-03 -> Ancora significativo dopo la correzione.

In sintesi abbiamo una singola osservazione che sembrerebbe essere un outlier critico. Pertanto procederemo controllando la Cook’s Distance: in questo modo se il punto ha Cook’s Distance alta, significa che ha un’influenza sproporzionata sul modello.

5.3 ANALISI COOK’S DISTANCE

La distanza di Cook è una misura utilizzata nell’analisi di regressione per identificare osservazioni influenti, cioè dati che hanno un impatto sproporzionato sulle stime dei parametri del modello di regressione. Quando un outlier supera un valore soglia della distanza di Cook, implica che quell’osservazione potrebbe avere un’influenza significativa sui risultati del modello, distorcendone i risultati, portando così a inferenze errate. È importante identificare e valutare queste osservazioni per assicurarsi che il modello sia robusto e rappresenti accuratamente i dati sottostanti.

cook <- cooks.distance(mod2)
plot(cook)

La maggior parte dei punti ha Cook’s Distance prossima allo 0, il che significa che non stanno influenzando molto la regressione. C’è un punto fuori scala, con Cook’s Distance molto alta:

max(cook)
## [1] 0.8300965

Identifichiamo il punto in questione

which(cooks.distance(mod2) > 0.5)
## 1551 
## 1551

Dato che il punto 1551 ha Cook’s Distance alta, leverage alto e residuo molto grande, allora possiamo considerarlo altamente influenzale e quindi deve essere analizzato con attenzione.

Ricordiamo:

Analizziamo con attenzione il punto individuato in modo da comprendere se:

dati[1551, ]
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551         35            1         1         38 4370       315    374
##      Tipo.parto Ospedale Sesso
## 1551        Nat     osp3     F

Stando ai valori riportanti dall’OMS, secondo i quali, tipicamente:

Possiamo facilmente concludere che ci troviamo dinanzi ad una neonata con un peso ed una lunghezza rispettivamente sopra e sotto la media.

Nonostante ciò, non abbiamo elementi sufficienti per affermare che si tratti di un errore di inserimento dati o di un’anomalia, imputabile ad un errore umano, pertanto verrà considerato un valore legittimo rappresentante un caso raro e quindi importante ai fini dell’indagine.

  1. PREVISIONI E RISULTATI

Giunti a questo punto, possiamo utilizzare il nostro modello per fare previsioni pratiche: ad esempio, stimeremo il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana:

nuovo_caso <- data.frame(
  Sesso='F', 
  N.gravidanze=3,
  Gestazione=39, 
  Lunghezza=round(mean(Lunghezza),2), 
  Cranio=round(mean(Cranio),2)
  )

peso_previsto <- predict(mod2, newdata = nuovo_caso)

cat("Il peso stimato del neonato è:", round(peso_previsto, 2), "grammi\n")
## Il peso stimato del neonato è: 3271.08 grammi
  1. RAPPRESENTAZIONE GRAFICA

Adesso, mediante grafici visualizzeremo i risultati del modello, mostrando le relazioni più significative tra le variabili

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 = "red",se=F, method = "lm")

Il peso alla nascita cresce all’aumentare delle settimane di gestazione, il che è coerente con la fisiologia neonatale.

Differenza tra maschi e femmine:

I maschi mostrano un incremento del peso più pronunciato all’aumentare della gestazione.

ggplot(data=dati)+
  geom_point(aes(x=N.gravidanze,
                 y=Peso,
                 col = Sesso), position = "jitter")+
  geom_smooth(aes(x=N.gravidanze,
                  y=Peso,
                  col =Sesso), se=F, method = "loess")+
  geom_smooth(aes(x=N.gravidanze,
                y=Peso), col = "red",se=F, method = "loess")

Gravidanze successive alla prima favorirebbero una crescita fetale leggermente superiore, fenomeno che sembrerebbe arrestarsi dopo la terza, il che è coerente con studi medici che indicano che un alto numero di gravidanze può essere associato a un maggior rischio di basso peso alla nascita. Tuttavia dato il basso numero di nascite oltre la quarta gravidanza del nostro campione, non è da escludersi che la causa sia da ricercarsi altrove.

ggplot(data=dati)+
  geom_point(aes(x=Lunghezza,
                 y=Peso,
                 col = Sesso), position = "jitter")+
  geom_smooth(aes(x=Lunghezza,
                  y=Peso,
                  col =Sesso), se=F, method = "lm")+
  geom_smooth(aes(x=Lunghezza,
                y=Peso), col = "red",se=F, method = "lm")

Il peso alla nascita cresce all’aumentare della lunghezza, il che è coerente con la fisiologia neonatale.

Differenza tra maschi e femmine:

I maschi mostrano un incremento del peso più pronunciato all’aumentare della lunghezza.

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")

Infine visualizziamo l’impatto del fumo sullo sviluppo del peso: le linee di regressione per fumatrici e non fumatrici sono diverse, pertanto il fumo potrebbe avere un impatto sul peso del neonato. Tuttavia dato il numero piccolissimo di fumatrici nel nostro campione non possiamo affermare con assoluta certezza una simile ipotesi