Funzioni utili

#calcolo indici di posizione e forma
calc_index<-function(x){
  x <- na.omit(as.numeric(x))
  return(c(
    Min=min(x),
    Max=max(x),
    Mean=mean(x),
    Median=median(x),
    Q1=quantile(x)[2],
    Q3=quantile(x)[4],
    Dev.st = sd(x), 
    IQR = IQR(x),
    Skewness = skewness(x),
    Kurtosis = kurtosis(x)-3
  ))
  t(x)
}

#formattazione delle tabelle
format_tbl<-function(table,variables_name,title){
  tbl_frm <-table %>%
    kable(
      format = "pipe", 
      digits = 2,
      caption = title,
      col.names = c(variables_name),
    ) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
}

Contesto Aziendale

Azienda: Neonatal Health Solutions Obiettivo: Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

1. Raccolta dei Dati e Struttura del Dataset

Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

2. Analisi e Modellizzazione

Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

Di seguito controlliamo che non ci siano valori nulli, visualiziamo una sintesi della caratterisitche dei dati ed estrapoliamo i primi 10 elementi dal daset:

dati <- read.csv("neonati.csv")
#verifico presenza eventuali dati nulli
sum(is.na(dati))
## [1] 0
kable(summary(dati))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235 Length:2500 Length:2500 Length:2500
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character Class :character Class :character
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00 Median :3300 Median :500.0 Median :340 Mode :character Mode :character Mode :character
Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98 Mean :3284 Mean :494.7 Mean :340 NA NA NA
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350 NA NA NA
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 NA NA NA 
kable(head(dati,10))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 42 3380 490 325 Nat osp3 M
21 2 0 39 3150 490 345 Nat osp1 F
34 3 0 38 3640 500 375 Nat osp2 M
28 1 0 41 3690 515 365 Nat osp2 M
20 0 0 38 3700 480 335 Nat osp3 F
32 0 0 40 3200 495 340 Nat osp2 F
26 1 0 39 3100 480 345 Nat osp3 F
25 0 0 40 3580 510 349 Nat osp1 M
22 1 0 40 3670 500 335 Ces osp2 F
23 0 0 41 3700 510 362 Ces osp2 F

Analizzando le variabili che compongno il dataset si possono individuare:

  • variabili qualitative nominali (cetegoriche): tipo.parto, ospedale, sesso, fumatrici
  • variabili quantitative discrete: numero gravidanze,gestazione
  • variabili quantitative continue: anni.madre , peso, Lunghezza, cranio

Al fine di gestire al meglio le variabili qualitative e poterle analizzare insieme alle altre variabili si procede all’importazione del dataset con l’argomento stringsAsFactors = T. In qesto modo le variabili qualitative presenti saranno considerate come dei fattori, e gestite come variabili categoriche vere e proprie e non semplice testo.

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

Ora abbiamo la ripartizione delle osservazioni per ogni possibile valore che assumono le variabili qualitative. Di seguito procederemo con un analisi descrittiva delle principali caratteristiche delle variabili presenti nel data set.

Analisi variabili quantitative

Analisi età madre

idx_eta<-calc_index(Anni.madre)
kable(idx_eta,
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Anni madre"
)
Indici di forma e variabilità - Anni madre
Valori
Min 0.0000000
Max 46.0000000
Mean 28.1640000
Median 28.0000000
Q1.25% 25.0000000
Q3.75% 32.0000000
Dev.st 5.2735783
IQR 7.0000000
Skewness 0.0428115
Kurtosis 0.3804165 
plot(density(Anni.madre),
  main= 'Distribuzione Anni madre',
  xlab="Età",
  ylab="Frequenza"
     )

Per quanto riguarda l’età delle madri si riscontra un’anomalia dovuta a dei valori errati (outliers) al di sotto del valore di 10, come si vede anche dal grafico e dal valore minimo = 0.

I valori della media e della mediana sono molto simili a confermare che i valori delle età delle madri seguono una distribuzione pressochè normale.

Si procede con l’escludere le osservazioni errate dal dataset e rifacciamo l’analisi descrittiva.

dati <- dati %>%
  filter(Anni.madre > 10)
attach(dati)
## I seguenti oggetti sono mascherati da dati (pos = 3):
## 
##     Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
##     Ospedale, Peso, Sesso, Tipo.parto
idx_eta<-calc_index(Anni.madre)
kable(idx_eta,
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Anni madre"
)
Indici di forma e variabilità - Anni madre
Valori
Min 13.0000000
Max 46.0000000
Mean 28.1861489
Median 28.0000000
Q1.25% 25.0000000
Q3.75% 32.0000000
Dev.st 5.2172061
IQR 7.0000000
Skewness 0.1510624
Kurtosis -0.1056061 
plot(density(Anni.madre),
  main= 'Distribuzione Anni madre',
  xlab="Età",
  ylab="Frequenza"
     )

Il valore minimo dell’età delle madri è 13 anni. La distribuzione presenta un’asimmetria positiva, ma molto prossima allo zero e una curtosi leggermente sopra lo zero. Media e mediana sono praticamente uguali.

Analisi numero gravidanze

kable(calc_index(N.gravidanze),
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - N° gravidanze"
)
Indici di forma e variabilità - N° gravidanze
Valori
Min 0.0000000
Max 12.0000000
Mean 0.9815853
Median 1.0000000
Q1.25% 0.0000000
Q3.75% 1.0000000
Dev.st 1.2809489
IQR 1.0000000
Skewness 2.5134123
Kurtosis 10.9816269 
barplot(table(N.gravidanze),
  main= 'Distribuzione N° gravidanze',
  xlab="N° gravidanze",
  ylab="Frequenza"
     )

Il grafico mostra come la distribuzione dei valori nelle classi di frequenza si concentri maggiormente per valori di numeri pari a 0 o 1. Tale distribuzione è confermata anche dal valore di asimmetria positivo, che indica un maggior numero di valori bassi. La distribuzione dei valori mostra un IQR di 1 con i valori che vanno da 0 a 1, rispettivamente il Q1 e il Q3. Il valore alto di curtosi è dovuto ai valori estremi che appiattiscono la distribuzione.

Analisi tempo gestazione

kable(calc_index(Gestazione),
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Gestazione"
)
Indici di forma e variabilità - Gestazione
Valori
Min 25.000000
Max 43.000000
Mean 38.979584
Median 39.000000
Q1.25% 38.000000
Q3.75% 40.000000
Dev.st 1.868950
IQR 2.000000
Skewness -2.065131
Kurtosis 8.255516 
barplot(table(Gestazione),
  main= 'Distribuzione Gestazione',
  xlab="Settimane",
  ylab="Frequenza"
     )

Dal grafico si può vedere come i valori si concentrino principalmente nelle 39 e 40 settimane di gestazione. La distribuzione è asimmetrica verso i valori più alti, confermato dal valori di skwness negativi. Media e mediana sono praticamente uguali (39 settimane), la moda è sule 40 settimane. I valori sono maggiormente distribuiti tra le 38 (Q1) e le 40 (Q3) settimane. l’indice di curtosi è alto, a causa dei valori estremi che si dicostano dal valore medio, che ne appiattiscono la curva.

Analisi Peso

kable(calc_index(Peso),
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Peso"
)
Indici di forma e variabilità - Peso
Valori
Min 830.0000000
Max 4930.0000000
Mean 3284.1841473
Median 3300.0000000
Q1.25% 2990.0000000
Q3.75% 3620.0000000
Dev.st 525.2293743
IQR 630.0000000
Skewness -0.6474036
Kurtosis 2.0287531 
plot(density(Peso),
  main= 'Distribuzione Peso',
  xlab="Peso (g)",
  ylab="Frequenza"
     )

Controlliamo che la variabile risposta peso, si ditribuisca secondo una distribuzione normale tramite lo Shapiro test:

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

La variabile di risposta Peso presenta una distribuzone leggermente asimmetrica negativa e una curtosi di tipo platicurtica. Prendendo in considerazione il valore dello Shapiro test si deve respingere l’ipotesi nulla, pertanto la distribuzione non è Normale.

Tuttavia, sia osservando il grafico della distribuzione, sia il discostamento dei valori della media (3284g) rispetto alla mediana (3300g) si ritiene che la distribuzione possa essere considerata come normale, per la realizzazione dei modelli lineari.

Analisi Lunghezza

kable(calc_index(Lunghezza),
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Lunghezza"
)
Indici di forma e variabilità - Lunghezza
Valori
Min 310.000000
Max 565.000000
Mean 494.695757
Median 500.000000
Q1.25% 480.000000
Q3.75% 510.000000
Dev.st 26.328847
IQR 30.000000
Skewness -1.514575
Kurtosis 6.480930 
plot(density(Peso),
  main= 'Distribuzione Lunghezza',
  xlab="Lunghezza (cm)",
  ylab="Frequenza"
     )

La variabile lunghezza presenta una asimmetria negativa, che indica che la distribuzione è spostata più versoi valori alti. la curtosi positiva è indice che la variabile si ditribuisce in maniera più ampia rispetto a una distribuzione normale, il che è dovuto alla presenza di valori estremi (minimo = 310, massimo = 565). La media e mediana differiscoono tra loro, ma di poco, e i valori hanno un range IQR di 30 cm, che indica che il 50% dei valori si aggira tra i 480 cm e i 510 cm.

Analisi diametro cranio

kable(calc_index(Cranio),
  col.names = c("Valori"),
  caption = "Indici di forma e variabilità - Cranio"
)
Indici di forma e variabilità - Cranio
Valori
Min 235.0000000
Max 390.0000000
Mean 340.0292234
Median 340.0000000
Q1.25% 330.0000000
Q3.75% 350.0000000
Dev.st 16.4294692
IQR 20.0000000
Skewness -0.7850906
Kurtosis 2.9448704 
plot(density(Peso),
  main= 'Distribuzione Dimensione Cranio',
  xlab="Cranio (mm)",
  ylab="Frequenza"
     )

I valori della dimensione del cranio si distribuiscono secondo una distribuzione che presenta un’asimmetria negativa. Media e mediana sono uguali, e si assestano su un valore di 340 mm. la variabilità dei valori non è molto elevata e il 50% di essi si concetra tra un range IQR che va da 330 mm a 350 mm, e un valore basso di deviazione standard pari a 16 mm.

Analisi variabili qualitative

di seguito si analizzeranno le variabili qualitative, al fine di capire come si distribuiscono i valori di ciascuna di esse:

par(mfrow=c(1, 2))
barplot(table(as.factor(Fumatrici)), xlab = 'Fumatrici [0 = no 1 = si]', ylab = 'Freq.')
barplot(table(Tipo.parto), xlab = 'Tipo parto', ylab = 'Freq.')

t_fum<-summary(as.factor(Fumatrici))
t_par<-summary(Tipo.parto)

kable(summary(as.factor(Fumatrici)),
  col.names = c("Totale"),
  caption = 'Distribuzione valori Fumatrici [0 = no 1 = si]',
  )
Distribuzione valori Fumatrici [0 = no 1 = si]
Totale
0 2394
1 104 
kable(summary(Tipo.parto),
  col.names = c("Totale"),
  caption = 'Distribuzione valori Tipo parto')
Distribuzione valori Tipo parto
Totale
Ces 728
Nat 1770 
par(mfrow=c(1, 2))
barplot(table(Ospedale), xlab = 'Ospedale',ylab = 'Freq.')
barplot(table(Sesso), xlab = 'Sesso',ylab = 'Freq.')

kable(summary(Ospedale),
  col.names = c("Totale"),
  caption = 'Distribuzione valori tra i vari ospedai')
Distribuzione valori tra i vari ospedai
Totale
osp1 816
osp2 848
osp3 834 
kable(summary(Sesso),
  col.names = c("Totale"),
  caption = 'Distribuzione valori per sesso')
Distribuzione valori per sesso
Totale
F 1255
M 1243

Verifica di alcune ipotesi

Si saggeranno le seguenti ipotesi con i test adatti:

1. in alcuni ospedali si fanno più parti cesarei

Osservizmo come si ditribuiscono i tipi di parto nei vari ospedali tramite una tabella di contingenza

cfr_ospedali<-table(Tipo.parto,Ospedale)
kable(cfr_ospedali)
osp1 osp2 osp3
Ces 242 254 232
Nat 574 594 602

Ora tramite il test chi-quadro verifichiamo se ci sono differenze significative tra i tipi di parto nei vari ospedali

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

dal testi del chi-quadro il p-value è al di sopra del valore di significatività del 5%, pertanto non si può rifiutare l’ipotesi nulla, e pertanto tra i vari ospedali non si ha una differenza significativa per quanto riguarda il tipo di parto effettuato. Possiamo verificare anche graficamente se ci sono eventuali associazioni significative tra i tipi di parto e i vari ospedali:

ggpubr::ggballoonplot(data=as.data.frame(cfr_ospedali))

Anche graficamente si può vedere che non ci sono significative associazioni tra tipo di parto e i vari ospedali.

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

Al fine di testare se le lunghezze presenti in questo campione sono significativamente diverse da quelle della popolazione si effettuerà il Test-T, utilizzando come parametro medio della popolazione il valore di 50 cm

t.test(Lunghezza, 
       mu = 500,
       conf.level = 0.95, 
       alternative = "two")
## 
##  One Sample t-test
## 
## data:  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

Il test restituisce un p-value inferiore alla soglia di significatività, pertanto, si rifiuta l’ipotesi nulla. La media delle lunghezze del campione è significativamente diversa dalla media della popolazione.

Allo stesso modo si verifica se il peso medio dei neonati di questo campione è uguale al valore medio della popolazione, fissato sui 3300 grammi

t.test(Peso,
       mu = 3300,
       conf.level = 0.95, 
       alternative = "two"
       )
## 
##  One Sample t-test
## 
## data:  Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

Per quanto riguarda il peso si è ottenuto un p-value di 0.1296 che superiore all’intervallo di confidenza del 5%, pertanto in questo caso si accetta l’ipotesi nulla e il peso medio del compione non è statisticamente diverso dal peso medio della popolazione.

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

Analisi correlazioni tra dimensione cranio e sesso

boxplot(
  Cranio~Sesso, 
  ylab ='Diametro (mm)', 
  col = c('red','royalblue2'),
  boxwex=.2
  )

t.test(Cranio~Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

Dai test effettuati si concluide che, tra tra maschi e femmine, i valori medi del diametro del cranio hanno differenze statisticamente rilevanti. L’ipotesi nulla viene rifiutata e le differenze sono significative.

Analisi correlazioni tra peso e sesso

boxplot(
  Peso~Sesso, 
  ylab ='Peso (g)', 
  col = c('red','royalblue2'),
  boxwex=.2)

t.test(Peso~Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496

Sia dal grafico che dal risultato del T-Test si può affermare che le medie dei pesi, dei due gruppi maschi e femmine, sono statisticamente differenti, si rifiuta l’ipotesi nulla.

Confronto tra lunghezza e sesso

boxplot(
  Lunghezza~Sesso, 
  ylab ='Peso (g)', 
  col = c('red','royalblue2'),
  boxwex=.2)

t.test(Lunghezza~Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750

Anche per la lunghezza, le differenze tra i valori medi, tra maschi e femmine, sono significative.

CREAZIONE DEL MODELLO DI REGRESSIONE LINEARE

Analisi correlazioni tra variabili quantitative

Ora analizziamo le correlazioni tra variabili quantitative, escludendo i fattori:

dati_num<- dati[,c("Anni.madre","N.gravidanze","Gestazione","Peso","Lunghezza","Cranio")]
shapiro.test(dati_num$Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  dati_num$Peso
## W = 0.97068, p-value < 2.2e-16
#?pairs()

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr")
  par(usr = c(0, 1, 0, 1))
  r <- (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 = 1.5)
}

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

Osservando i valori di correlazione tra le variabili si può osservare: - la variabile numero di gravidanze ha una bassa correlazione positiva con l’età delle madri. Maggiore è l’età della madre, maggiore è il numero di gravidanze. Si evidenziano valori leverages con l’età della madre - la durata della gestazione ha una correlaizone positiva, anche se non del tutto lineare, con il peso, la lunghezza e la dimensione del cranio del neonato. La correlazione maggiore la si ha con la lunghezza del neonato. Maggiore è il numero di settimane di gestazione, maggiori sono le dimensioni del neonato. Si possono osservare dei valori leverages per quanto riguarda la lunghezza e la dimensione del cranio. - la lunghezza del neonato è correlata positivamente con il peso, la dimensione del cranio e la durata della gestazione. - la variabile risposta peso è correlata positivamente con la durata della gestazione, la lungheza del neonato e la dimennsione del cranio, per i quali fa registrare una forte correlazione lineare. Si osservano valori outliner.

Le altre correlazioni sono basse e prossime allo zero, pertanto sono da ritenersi poco significative.

Analisi correlazioni tra variabili qualitative e variabile risposta Peso

Boxplot realtivi alla variabile peso e le variabili qualitative: fumatrici, sesso, tipo parto e ospedale:

par(mfrow=c(1,2))
boxplot(Peso~as.factor(Fumatrici), xlab='Fumatrici [0 = no 1 = si]')
boxplot(Peso~Sesso, xlab='Sesso')

par(mfrow=c(1,2))
boxplot(Peso~Tipo.parto, xlab='Tipo parto')
boxplot(Peso~Ospedale, xlab='Ospedale')

graficamente si può osservare che il tipo di parto e l’ospedale non influenzano la variabile peso. Invece sembrerebbe esserci uno scostamento per quanto riguarda le variabili sesso e fumatrici.

Per queste ultime, verifichiamo tramite il T-test se c’è una differenza statisticamente significativa dei valori medi del peso, e se pertanto hanno una reale influenza su di esso.

t.test(Peso~Fumatrici)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by 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
t.test(Peso~Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496

I valori dei test indicano che c’è una significativa differenza tra i valori medi del peso tra neonati maschi e femmine, metre per le madri fumatrici e non fumatrici, non si può rifiutare ’ipotesi nulla, pertanto non ci sono differenze significative per quanto riguarda i pesi medi.

Modello Lineare 1

Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

Le variabili ritenute rilevanti sono: Anni della madre, fumetrici, periodo gestazione, sesso,lunghezza e dimensione cranio.

Si decide di escludere ospedale e tipo parto in quanto irrilevanti per quanto riguarda la crescita del feto o comunque il peso del neonato.

mod_1<-lm(Peso~ Anni.madre+as.factor(Fumatrici)+Gestazione+Sesso+Lunghezza+Cranio, data=dati)
summary(mod_1)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + as.factor(Fumatrici) + Gestazione + 
##     Sesso + Lunghezza + Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.76  -184.44   -14.07   162.78  2615.42 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -6723.6155   141.3960 -47.552  < 2e-16 ***
## Anni.madre                1.9110     1.0692   1.787    0.074 .  
## as.factor(Fumatrici)1   -27.0218    27.6004  -0.979    0.328    
## Gestazione               32.4434     3.8270   8.478  < 2e-16 ***
## SessoM                   78.8488    11.2198   7.028  2.7e-12 ***
## Lunghezza                10.2001     0.3011  33.872  < 2e-16 ***
## Cranio                   10.6028     0.4261  24.882  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2491 degrees of freedom
## Multiple R-squared:  0.7266, Adjusted R-squared:  0.7259 
## F-statistic:  1103 on 6 and 2491 DF,  p-value: < 2.2e-16
kable(vif(mod_1), caption = 'Verifica collinearità')
Verifica collinearità
x
Anni.madre 1.027671
as.factor(Fumatrici) 1.004152
Gestazione 1.689377
Sesso 1.039664
Lunghezza 2.075902
Cranio 1.618613

Per il modello 1 la variabile Fumatrici è quella meno significativa, seguita da anni madre.

Si decide di procedere con un sistema step by step ed eliminare una variabile alla volta e vedere come viene influenzato il valore di R-Quadro aggiustato.

Modello 2

mod_2<-lm(Peso~ Anni.madre+Gestazione+Sesso+Lunghezza+Cranio, data=dati)
summary(mod_2)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Sesso + Lunghezza + 
##     Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.03  -184.20   -14.07   163.24  2618.69 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6723.2290   141.3943 -47.550  < 2e-16 ***
## Anni.madre      1.8995     1.0692   1.777   0.0758 .  
## Gestazione     32.2256     3.8205   8.435  < 2e-16 ***
## SessoM         78.6738    11.2182   7.013 2.99e-12 ***
## Lunghezza      10.2137     0.3008  33.954  < 2e-16 ***
## Cranio         10.6047     0.4261  24.887  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.7265, Adjusted R-squared:  0.7259 
## F-statistic:  1324 on 5 and 2492 DF,  p-value: < 2.2e-16
kable(vif(mod_2), caption = 'Verifica collinearità')
Verifica collinearità
x
Anni.madre 1.027548
Gestazione 1.683669
Sesso 1.039400
Lunghezza 2.071502
Cranio 1.618576

Il modello 2 fa resistrare un R quadro aggiustato uguale a quello del modello 1, la variabile meno significativa ora è anni madre. Si procede pertanto a eliminarla creando un terzo modello.

Modello 3

mod_3<-lm(Peso~ Gestazione+Sesso+Lunghezza+Cranio, data=dati)
summary(mod_3)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Sesso + Lunghezza + Cranio, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.1  -184.4   -17.4   163.3  2626.7 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.6732   135.5952 -49.055  < 2e-16 ***
## Gestazione     31.3262     3.7884   8.269  < 2e-16 ***
## SessoM         79.1027    11.2205   7.050 2.31e-12 ***
## Lunghezza      10.2024     0.3009  33.909  < 2e-16 ***
## Cranio         10.6706     0.4247  25.126  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275.1 on 2493 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1652 on 4 and 2493 DF,  p-value: < 2.2e-16
kable(vif(mod_3), caption = 'Verifica collinearità')
Verifica collinearità
x
Gestazione 1.654101
Sesso 1.038918
Lunghezza 2.070582
Cranio 1.606316

il modello 3 fa registrare una diminuzione del R quadro aggiustato, ma la diminuzione è minima 0.0002, in compenso si utilizzano meno variabili.

Selezione per step wise

Giusto per prova, procediamo con la creazione e selezione del modello lineare tramite la funzione stepwise di R, a partire dal modello 1.

library(MASS)
## 
## Caricamento pacchetto: 'MASS'
## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     select
mod_stp <- stepAIC(mod_1, direction = "both")
## Start:  AIC=28067.98
## Peso ~ Anni.madre + as.factor(Fumatrici) + Gestazione + Sesso + 
##     Lunghezza + Cranio
## 
##                        Df Sum of Sq       RSS   AIC
## - as.factor(Fumatrici)  1     72476 188424446 28067
## <none>                              188351970 28068
## - Anni.madre            1    241517 188593487 28069
## - Sesso                 1   3734388 192086358 28115
## - Gestazione            1   5434264 193786234 28137
## - Cranio                1  46812401 235164371 28621
## - Lunghezza             1  86752952 275104922 29012
## 
## Step:  AIC=28066.94
## Peso ~ Anni.madre + Gestazione + Sesso + Lunghezza + Cranio
## 
##                        Df Sum of Sq       RSS   AIC
## <none>                              188424446 28067
## + as.factor(Fumatrici)  1     72476 188351970 28068
## - Anni.madre            1    238661 188663107 28068
## - Sesso                 1   3718775 192143221 28114
## - Gestazione            1   5379726 193804172 28135
## - Cranio                1  46830959 235255405 28619
## - Lunghezza             1  87168676 275593122 29015
summary(mod_stp)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Sesso + Lunghezza + 
##     Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.03  -184.20   -14.07   163.24  2618.69 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6723.2290   141.3943 -47.550  < 2e-16 ***
## Anni.madre      1.8995     1.0692   1.777   0.0758 .  
## Gestazione     32.2256     3.8205   8.435  < 2e-16 ***
## SessoM         78.6738    11.2182   7.013 2.99e-12 ***
## Lunghezza      10.2137     0.3008  33.954  < 2e-16 ***
## Cranio         10.6047     0.4261  24.887  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.7265, Adjusted R-squared:  0.7259 
## F-statistic:  1324 on 5 and 2492 DF,  p-value: < 2.2e-16

La modelizzazione automatica indicherebbe il modello 2 come modello ottimale, ma andremo a indagare meglio la scelta del modello migliore.

Selezione del Modello Ottimale

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

Attraverso gli indici statistici AIC e BIC si valuterà quale modello è il migliore:

kable(AIC(mod_1,mod_2,mod_3))
df AIC
mod_1 8 35159.00
mod_2 7 35157.96
mod_3 6 35159.12 
kable(BIC(mod_1,mod_2,mod_3))
df BIC
mod_1 8 35205.58
mod_2 7 35198.72
mod_3 6 35194.06

I modelli sono pressochè identici dal punto di vista dei test. l’indice AIC indicherebbe il modello 2 come il migliore, mentre l’indice BIC il modello 3.

Proviamo ad effettuare un’analisi ANOVA sui modelli 2 e 3, per capire se la variabile Anni madre è utile e aumenta l’informazione del modello.

anova(mod_2,mod_3)
## Analysis of Variance Table
## 
## Model 1: Peso ~ Anni.madre + Gestazione + Sesso + Lunghezza + Cranio
## Model 2: Peso ~ Gestazione + Sesso + Lunghezza + Cranio
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1   2492 188424446                              
## 2   2493 188663107 -1   -238661 3.1564 0.07575 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L’Anova restituisce un p-value di 0.076, pertanto la variabile anni madre non è così significativa.

Pertanto si decide di utilizazre il modello 3 come meodello finale, anche tenendo conto che è preferible un modello con meno variabili a uno con più variabili.

Analisi della Qualità del Modello

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

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

Procediamo con l’analisi visiva dei residui del modello lineare.

Residuals vs Fitted: Per capire se la distribuzione dei residui segue la linearità, essi si devono disporre casualmente attorno al valore zero (linea trateggiata). Dal grafico si può osservare che la linea di tendenza (rossa) conferma la linearità dei residui, tranne per i dati di peso più bassi (la linea si curva). Questo indica che il modello non sarà in grado, o predirrà con minor precisione, i valori corrispondenti a pesi piccoli. inoltre si osserva la presenza di una valore di lavarage, il punto 1549

Q-Q Residuals: i dati si distribuiscono lungo la retta, tranne l’ulitmo tratto che si incurva verso l’alto. Anche in questo caso il valore 1549 è isolato dagli altri, e potrebbre essere la causa della curvatura dei valori, che causa la non normalità della distribuzione dei residui.

Scale-Location: i residui sono disposti in maniera costante, tranne nella prima parte del grafico con i valori più bassi di peso. Dove si registra un’alta eteroschedasticità

Residuals vs Leverage: il valore 15490 super la soglia della distanza di Cook dello 0,5 che rappresenta una soglia di allerta.

Dato che il valore 1549 potrebbe essere un valore influente, procediano con l’analisi dei valori di levarage e degli outliers. Di seguito individuiamo i punti di leverage

#leverage
lev<-hatvalues(mod_3)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
plot(lev)
abline(h=soglia,col=2)

lev[lev>soglia]
##          15          34          42          61          67          96 
## 0.006151861 0.006243278 0.004255297 0.004587531 0.005349019 0.004802045 
##         101         106         117         131         151         155 
## 0.007186795 0.012822506 0.004750378 0.006972348 0.010852877 0.006671488 
##         190         205         206         220         249         295 
## 0.005288680 0.005297151 0.009408332 0.007386096 0.004655804 0.004010796 
##         304         305         310         312         315         348 
## 0.004421683 0.005385519 0.028760605 0.013130123 0.005343752 0.004188947 
##         378         383         445         471         486         492 
## 0.015374067 0.004306002 0.007095381 0.004291248 0.004741540 0.008175742 
##         565         587         592         615         629         638 
## 0.004639536 0.008378669 0.006314805 0.004576590 0.004041835 0.006220559 
##         656         666         684         697         702         726 
## 0.004742368 0.004345956 0.008757517 0.005810826 0.004720571 0.004039128 
##         748         750         765         805         821         895 
## 0.008237128 0.006713862 0.006055540 0.014305924 0.004041356 0.005205549 
##         928         947         956         964         968         991 
## 0.022097377 0.007803683 0.007677273 0.004618560 0.004079825 0.004259154 
##        1014        1067        1091        1096        1130        1156 
## 0.008225437 0.007870792 0.008934430 0.004268185 0.006568161 0.004084916 
##        1165        1180        1187        1199        1237        1247 
## 0.004020912 0.005586745 0.006408532 0.005446218 0.005378187 0.014581054 
##        1272        1282        1292        1293        1355        1356 
## 0.007083865 0.004056888 0.005556170 0.004742521 0.005289297 0.006528860 
##        1360        1383        1393        1398        1400        1418 
## 0.004081173 0.012021412 0.004649879 0.005565737 0.004793538 0.005132039 
##        1426        1427        1549        1551        1554        1558 
## 0.008194247 0.021037931 0.048555718 0.006812048 0.005871766 0.004592075 
##        1591        1604        1608        1617        1626        1632 
## 0.004857727 0.004754528 0.007765891 0.014505134 0.004664490 0.004529497 
##        1684        1690        1691        1699        1710        1733 
## 0.008721997 0.004262114 0.005042673 0.010199722 0.006947104 0.004371431 
##        1778        1800        1804        1807        1825        1856 
## 0.025535302 0.004008990 0.004091435 0.008391698 0.006011201 0.004329297 
##        1866        1891        1909        1918        1975        2035 
## 0.004746934 0.004246662 0.004485587 0.004131744 0.005384719 0.004235451 
##        2038        2064        2087        2112        2113        2118 
## 0.011439028 0.004017066 0.006252500 0.012851847 0.011770955 0.018008805 
##        2138        2147        2173        2198        2213        2214 
## 0.006091129 0.013041062 0.021175260 0.011037915 0.004835341 0.007550907 
##        2218        2222        2223        2255        2305        2316 
## 0.004026782 0.005782813 0.005410507 0.005768404 0.013904929 0.004773600 
##        2335        2357        2389        2406        2435        2450 
## 0.004700877 0.004750715 0.004387051 0.009638225 0.023764159 0.023235828 
##        2456        2474        2476 
## 0.008005109 0.004070657 0.005747348

si procede con l’analisi degli outliners

plot(rstudent(mod_3))
abline(h=c(-2,2))

outlierTest(mod_3)
##      rstudent unadjusted p-value Bonferroni p
## 1549 9.980701         4.9791e-23   1.2438e-19
## 155  4.948947         7.9596e-07   1.9883e-03
## 1305 4.779033         1.8638e-06   4.6558e-03

Valutiamo la distanza di Cook per evidenziare i valori più influente:

N=nrow(dati) #salvo il numero di osservazioni
cooks<-cooks.distance(mod_3) 

plot(cooks,ylim = c(0,1)) 
max(cooks) 
## [1] 0.9780498
text(x=2:length(cooks), y=cooks-0.05, labels=ifelse(cooks== max(cooks), names(cooks),""), col="red")  # add labels

L’analisi sui leverages e suigli outliners conferma che il valore 1549 è un valore influente.

Eliminiamo il valore influente dai dati e riverifichiamo il modello

influential <- as.numeric(names(cooks)[(cooks== max(cooks))])
dati<-dati[-influential, ]
mod_3<-lm(Peso~Gestazione+Sesso+Lunghezza+Cranio, data=dati)
summary(mod_3)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Sesso + Lunghezza + Cranio, 
##     data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1153.75  -182.01   -14.88   163.78  1391.18 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6652.115    132.990 -50.020  < 2e-16 ***
## Gestazione     28.533      3.726   7.657 2.69e-14 ***
## SessoM         79.321     11.005   7.208 7.51e-13 ***
## Lunghezza      10.841      0.302  35.903  < 2e-16 ***
## Cranio         10.059      0.421  23.893  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.8 on 2492 degrees of freedom
## Multiple R-squared:  0.7362, Adjusted R-squared:  0.7358 
## F-statistic:  1739 on 4 and 2492 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2)) 
plot(mod_3)

Eliminando il valore outliers si è amentato il valore di R quadro corretto e nessun valore supera la soglia di 0,5 e di 1 della distanza di Cook. Pertanto si decide di mantenere la correzione effettuata al fine di migliorare il modello predittivo

effettuando i test sui residui parametrici

#test sui residui numerici
lmtest::bptest(mod_3) #omoschelasticità varianza costante non si rifiuta lì'ipotesi nulla
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_3
## BP = 9.3293, df = 4, p-value = 0.05338
lmtest::dwtest(mod_3) #darwin watson i residui non sono auto correlate, non si rifiuta l'ipotesi nulla
## 
##  Durbin-Watson test
## 
## data:  mod_3
## DW = 1.9555, p-value = 0.1325
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(mod_3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod_3$residuals
## W = 0.98868, p-value = 3.6e-13
plot(density(residuals(mod_3))) #grafico per vedere la distribuzione dei residui. Problema su coda sinistra.

Dall’analisi dei residui sul modello 3 senza outliner risulta che: - Il test di Breusch-Pagan indica la presenza di eteroschedasticità. - si rifiuta l’ipotesi nulla per lo Shapiro test, pertanto i residui non si distribuiscono normalmente.
- il test di Durbin-Watson indica che i residui non sono correlati.

3. Previsioni e Risultati

Una volta validato il modello, lo useremo per fare previsioni pratiche.

Stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

Come dati di lunghezza e dimensione del cranio usiamo i dati medi presi da letteratura.

parametri <- data.frame(N.gravidanze = 3, Gestazione = 39, Sesso = "F", Lunghezza=500, Cranio=350)
predict(mod_3,parametri)
##        1 
## 3401.994

4. Visualizzazioni

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili.

Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

ggplot(data = dati, aes(x = Gestazione, y = Peso),position='jitter') +
  geom_point(position='jitter') + 
  geom_smooth(aes(x=Gestazione, y=Peso), method='lm')+
  ggtitle("Settimane gestazione vs peso")
## `geom_smooth()` using formula = 'y ~ x'

La retta di regressione presenta un relazione lineare positiva. Per quanto riguarda le osservazioni con valori più bassi di gestazione, si posizonano al di sotto della linea di regressione. La maggior parte delle osservazioni si posizionano al disopra delle 36 settimane e sono ben distribuite sopra e sotto la linea di regressione.

boxplot(Peso~Fumatrici)

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

non ci sono differenze sginificative sui valori medi del peso per la condizione di madre fumatrice o non fumatrice.

Conclusioni

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

Il modello risulta essere un buono strumento per predire il peso dei neonati. Esso si basa su parametri facilmente rilevabili e di facile ottenimento sia durante il periodo di gestazione che dopo il parto.