Un modello statistico per prevedere il peso dei neonati

1) Importa il dataset “neonati.csv” e controlla che sia stato letto correttamente dal software

Viene riprotato di seguito lo script in R che crea il data.frame “dataset” a partire dal file “neonati.csv”

Il dataset contiene 2500 osservazioni di 10 variabili.

if (!require("dplyr")) install.packages("dplyr")
if (!require("ineq")) install.packages("ineq")
if (!require("e1071")) install.packages("e1071")
if (!require("ggplot2")) install.packages("ggplot2")

library(dplyr)
library(ineq)
library(ggplot2)
library(e1071)
library(ggplot2)

csv_neonati_path <- "neonati.csv"
dataset <- read.csv(csv_neonati_path, stringsAsFactor = T)

#elimina 2 osservazioni errate
dataset = subset(dataset, Anni.madre > 1)

dataset$Fumatrici <- ifelse(dataset$Fumatrici == 1, "SI", "NO")
dataset <- dataset[, c("Peso", "Lunghezza", "Cranio", "Anni.madre", "N.gravidanze", "Gestazione", "Fumatrici", "Tipo.parto", "Ospedale", "Sesso")]
attach(dataset)

# Selezionare solo le colonne numeriche
numeric_cols <- sapply(dataset, is.numeric)
dataset_numeric <- dataset[, numeric_cols]

n = nrow(dataset)

2) Descrivi il dataset, la sua composizione, il tipo di variabili e l’obiettivo dello studio

Nome variabile	Descrizione	Tipo
Anni.madre	Età della madre	quantitativa discreta
N.gravidanze	Numero di gravidanze precedenti	quantitativa discreta
Fumatrici	Indica se la madre è fumatrice (0=NO, SI=1)	qualitativa dicotomica
Gestazione	Durata della gestazione in settimane	quantitativa discreta
Peso	Peso del neonato alla nascita in grammi	quantitativa continua
Lunghezza	Lunghezza del neonato alla nascita in millimetri	quantitativa continua
Cranio	Circonferenza cranica del neonato in millimetri	quantitativa continua
Tipo.parto	Tipo di parto, ad esempio naturale o cesareo	qualitativa dicotomica
Ospedale	Ospedale in cui è avvenuta la nascita	qualitativa nominale
Sesso	Sesso del neonato, indicato come “M” (maschio) o “F” (femmina)	qualitativa dicotomica

3) Indaga le variabili effettuando una breve analisi descrittiva, utilizzando indici e strumenti grafici che conosci

Calcolo dei principali indici

Allo scopo sono state implementate in R due funzioni che calcolano, dato un set di osservazioni di una variabile, gli indici di posizione, variabilità e forma, e le distribuzioni di frequenze.

Le funzioni generiche implementate verranno utilizzate in funzione della tipologia di variabili trattate presenti nel dataset.

Di seguito viene riportato lo script in R.

# definizione della funzione custom_summary per il calcolo degli indici di posizione, variabilità e forma
# e della funzione frequency_distribution per il calcolo delle frequenze di distribuzione

custom_summary <- function(x, variable_name) {
  stat_df <- data.frame(
    Category = variable_name,
    Count = length(x),
    Mean = mean(x, na.rm = TRUE),
    Var = var(x, na.rm = TRUE),
    Std = sd(x, na.rm = TRUE),
    Coeff_var = sd(x) / mean(x),
    Min = min(x, na.rm = TRUE),
    Q25 = quantile(x, 0.25, na.rm = TRUE),
    Median = median(x, na.rm = TRUE),
    Q75 = quantile(x, 0.75, na.rm = TRUE),
    Max = max(x, na.rm = TRUE),
    Skewness = skewness(x, na.rm = TRUE),
    Kurtosis = kurtosis(x, na.rm = TRUE)
  )
  
  return(stat_df)
}

frequency_distribution <- function(x, variable_name) {
  freq_table <- table(x)
  rel_freq_table <- prop.table(freq_table)
  
  freq_df <- data.frame(
    Category = variable_name,
    Category_Value = names(freq_table),
    Freq_ass = as.integer(freq_table),
    Freq_rel = as.numeric(rel_freq_table)
  )
  
  return(freq_df)
}


# Calcolo degli indici di posizione, variabilità e forma per le variabili quantitative continue e discrete del dataset

summary_stats <- data.frame()

summary_stats <- rbind(summary_stats, custom_summary(dataset$Anni.madre, "anni_madre"))
summary_stats <- rbind(summary_stats, custom_summary(dataset$N.gravidanze, "numero_gravidanze"))
summary_stats <- rbind(summary_stats, custom_summary(dataset$Gestazione, "gestazione"))
summary_stats <- rbind(summary_stats, custom_summary(dataset$Peso, "peso"))
summary_stats <- rbind(summary_stats, custom_summary(dataset$Lunghezza, "lunghezza"))
summary_stats <- rbind(summary_stats, custom_summary(dataset$Cranio, "cranio"))


# Calcolo delle distribuzioni di frequenze per le variabili qualitative

freq_distr_fumatrici <- frequency_distribution(dataset$Fumatrici, "fumatrice")
freq_distr_tipo_parto <- frequency_distribution(dataset$Tipo.parto, "tipo_parto")
freq_distr_ospedale <- frequency_distribution(dataset$Ospedale, "ospedale")
freq_distr_sesso <- frequency_distribution(dataset$Sesso, "sesso")

print(freq_distr_fumatrici)
print(freq_distr_tipo_parto)
print(freq_distr_ospedale)
print(freq_distr_sesso)

---

La seguente tabella mostra le statistiche principali calcolate per le variabili quantitative continue e discrete del dataset

Category	Count	Mean	Var	Std	C_var	Min	Q25	Median	Q75	Max	Skewness	Kurtosis
anni_madre	2500	28.16	27.81	5.27	0.19	0	25	28	32	46	0.04	0.38
numero_gravidanze	2500	0.98	1.64	1.28	1.31	0	0	1	1	12	2.51	10.98
gestazione	2500	38.98	3.49	1.87	0.05	25	38	39	40	43	-2.06	8.25
peso	2500	3284.08	275665.68	525.04	0.16	830	2990	3300	3620	4930	-0.65	2.03
lunghezza	2500	494.69	692.67	26.32	0.05	310	480	500	510	565	-1.51	6.48
cranio	2500	340.03	267.79	16.43	0.05	235	330	340	350	390	-0.78	2.94

---

Output delle distribuzioni di frequenze per la variabile dicotomica Fumatrici freq_distr_fumatrici

Category	Category_Value	Freq_ass	Freq_rel
Fumatrici	NO	2396	0.9584
Fumatrici	SI	104	0.0416

Output delle distribuzioni di frequenze per la variabile dicotomica tipo_parto freq_distr_tipo_parto

Category	Category_Value	Freq_ass	Freq_rel
tipo_parto	Ces	728	0.2912
tipo_parto	Nat	1772	0.7088

Output delle distribuzioni di frequenze per la variabile qualitativa ospedale freq_distr_ospedale

Category	Category_Value	Freq_ass	Freq_rel
ospedale	osp1	816	0.3264
ospedale	osp2	849	0.3396
ospedale	osp3	835	0.3340

Output delle distribuzioni di frequenze per la variabile dicotomica sesso freq_distr_sesso

Category	Category_Value	Freq_ass	Freq_rel
sesso	F	1256	0.5024
sesso	M	1244	0.4976

A partire dai dati di cui sopra, possiamo effettuare le seguenti considerazioni:

• Anni.madre: L’età media delle madri è di circa 28 anni, con una deviazione standard di circa 5,3 anni. L’età minima registrata è 0 (probabilmente un errore), mentre la massima è 46 anni.

• N.gravidanze: La media del numero di gravidanze precedenti è di circa 1, con la maggior parte delle madri che hanno avuto tra 0 e 1 gravidanze precedenti. Il massimo è 12 gravidanze.

• Fumatrici: Solo il 4,2% delle madri è fumatrice, mentre il 95,8% non lo è.

• Gestazione: La durata media della gestazione è di circa 39 settimane, con una deviazione standard di 1,87 settimane. La gestazione più breve è stata di 25 settimane, mentre la più lunga di 43 settimane.

• Peso: Il peso medio dei neonati è di circa 3284 grammi, con una deviazione standard di 525 grammi. Il peso varia da un minimo di 830 grammi a un massimo di 4930 grammi.

• Lunghezza: La lunghezza media alla nascita è di circa 495 mm, con una deviazione standard di 26,3 mm.

• Cranio: La circonferenza cranica media è di circa 340 mm, con una deviazione standard di 16,4 mm. -

• Tipo.parto: La maggior parte dei parti è di tipo naturale (Nat).

• Ospedale: Gli ospedali sono ben distribuiti, con una leggera predominanza di alcuni rispetto ad altri.

• Sesso: La distribuzione tra maschi e femmine è quasi equamente bilanciata.

Grafici

Di seguito sono riportati alcuni grafici che aiutano a comprendere la distribuzione dei dati delle singole variabili, nonchè la correlazione tra le variabili stesse del dataset, con particolare focus rispetto la variabile Peso.

Dai grafici di cui sopra si riportano le seguenti considerazioni:

1. Distribuzione dell’età della madre

   • La distribuzione dell’età della madre ha una forma normale, come si evince anche dal valore di Skewness = 0.04.

   • La maggior parte dei parti avviene in un’eta che oscilla tra i 25 ed i 32 anni.

   • Si osservano 2 osservazioni con ogni probabilità errate, in quanto la madre risulta avere 0 ed 1 anno.

2. Distribuzione del numero di gravidanze

   • La distribuzione del numero di gravidanze è fortemente asimmetrica con una coda più lunga verso destra.

   • La maggior parte delle gravidanze rappresenta per le donne il primo parto.

3. Distribuzione del numero di settimane di gestazione:

   • La distribuzione delle settimane di gestazione è fortemente asimmetrica con una coda più lunga verso sinistra.

   • La maggior parte dei dati sono concentrati tra le 37 e le 40 settimane. Questo suggerisce che la maggior parte delle gravidanze in questo dataset raggiunge il termine (circa 37-40 settimane), mentre le gravidanze premature o post-termine sono meno frequenti.

   • Si nota un calo significativo della frequenza dopo le 40 settimane, in linea con la pratica clinica di indurre il parto in caso di gravidanza post-termine.

4. Distribuzione della lunghezza in mm:

   • La distribuzione delle lunghezze alla nascita sembra avere una forma normale, leggermente asimmetrica verso sinistra. La maggior parte dei neonati ha una lunghezza compresa tra 450 e 550 mm, con un picco intorno ai 500 mm.

   • Ci sono alcuni outlier (valori anomali) su entrambi i lati, che indicano che ci sono neonati significativamente più piccoli o più grandi rispetto alla media.

5. Distribuzione della circonferenza cranica:

   • La distribuzione della circonferenza cranica appare approssimativamente normale, centrata intorno ai 340-350 mm. Ciò indica che la maggior parte dei neonati ha una circonferenza cranica entro un intervallo relativamente ristretto.

   • Sono presenti alcuni outlier con circonferenze craniche molto piccole o molto grandi, che potrebbero indicare condizioni mediche particolari o una popolazione di neonati più diversificata.

---

Per visualizzare graficamente la distribuzione del peso alla nascita, la variabile Peso è stata divisa in classi da 150 grammi.

Analizzando l’istogramma sulla distribuzione del peso alla nascita, possiamo fare alcune osservazioni:

1. Distribuzione approssimativamente normale: Il grafico mostra una distribuzione che tende ad avere una forma normale. C’è una chiara concentrazione di nascite intorno a un peso centrale, con variazioni che oscillano tra 2930 e 3380 grammi, suggerendo che il peso medio alla nascita si trova in questo intervallo.

2. Asimmetria leggera verso sinistra: C’è una leggera asimmetria a sinistra (distribuzione con “coda” più lunga verso pesi inferiori). Questo indica che c’è una piccola proporzione di neonati con peso significativamente inferiore rispetto alla media, ma non tale da alterare drasticamente la simmetria complessiva.

3. Frequenza massima: L’intervallo di peso più comune (moda) si trova tra 2930 e 3080 grammi, con un picco di 354 neonati in quell’intervallo. Gli intervalli vicini, come quelli tra 3080-3230 grammi e 2780-2930 grammi, presentano anch’essi frequenze relativamente alte (315 e 267 rispettivamente).

4. Pesi estremi: Sono presenti pochi casi di neonati con pesi estremamente bassi (meno di 1000 grammi) e pesi estremamente alti (più di 4880 grammi), ma questi sono rari, con pochi valori ai margini. Questo potrebbe rappresentare neonati prematuri o con complicazioni alla nascita

---

Di seguito verranno visualizzate le correlazioni tra le variabili, con particolare riferimento alla variabile Peso, tramite la matrice di correlazione per le variabili quantitative.

La matrice mostra una evidente correlazione tra la variabile Peso e le variabili Lunghezza (0.8), Cranio (0.7) ed, in misura ridotta, con la variabile Gestazione (0.59).

Tali correlazioni sono evidenziati anche dai relativi scatterplot in cui si evidenzia un pattern riconoscibile che indica una chiara relazione lineare positiva, il che indica che l’aumento della lunghezza del neonato, piuttosto che la circonferenza cranica od il numero di settimane di gestazione, porta ad un aumento anche il peso.

Per la variabile Gestazione, in particolare, si evidenzia ad un certo punto un abbassamento della positività della relazione il che suggerisce che il numero di settimane di gestazione ad un certo punto influisce di meno sul peso del neonato.

---

I seguenti boxplot mostrano invece la correlazione tra la variabile Peso e le variabili qualitative del dataset

I grafici mostrano che i neonati maschi pesano di più delle femmine, ed il peso alla nascita dei neonati partoriti da donne che dichiarano di fumare, sia minore rispetto le donne non fumatrici.

Non si evidenziano relazioni particolari tra il peso e l’ospedale oggetto del parto.

Stessa valutazione per la tipologia di parto.

Pulizia del dataset

Sull base di quanto osservato, si è deciso di eliminare le osservazioni 1380 ed 1152 i cui valori della variabile Anni.madre sono palesemente errati (0 e 1 anno).

Tutti i test ed i modelli di regressione di seguito riportati faranno uso del dataset pulito.

4) Saggia l’ipotesi che la media del peso e della lunghezza di questo campione di neonati siano significativamente uguali a quelle della popolazione

Allo scopo verrà eseguito un test t per il peso ed uno per la lunghezza, con il seguente comune sistema di ipotesi:

• Per il peso / lunghezza:

  • H0 (ipotesi nulla) = H0 ​: μ=μ0

  • H1 (ipotesi alternativa)​ = H1 ​: μ!= μ0​

Dove:

• μ è la media del peso / lunghezza dei neonati nel campione.

• μ0​ è la media del peso / lunghezza dei neonati nella popolazione.

Fissiamo, per entrambi i test t, un livello di significativà α = 0.05.

Ciò significa che c’è una probabilità del 5% di rifiutare l’ipotesi nulla quando essa è vera.

Test t per la variabile peso

Per eseguire il test t supponiamo che la media del peso alla nascita dei bambini della popolazione sia 3300 grammi (fonte Crescita bimbi in Italia)

Il risultato del test è il seguente:

One Sample t-test

data:  dataset$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 

Con un valore di p-value = 0.1324, essendo maggiore di alfa (0.05), non rifiutiamo l’ipotesi nulla, quindi possiamo continuare ad affermare, fino a prova contraria, che la media del campione sia significativamente uguale a quella della popolazione.

Ciò è evidenziato anche dal 95% di confidenza dell’intervallo di accettazione dell’ipotesi, pari a 3263.577 - 3304.791 per cui il valore di media della popolazione (3300) ricade al suo interno.

Test t per la variabile lunghezza

Per eseguire il test t supponiamo che la media della lunghezza alla nascita dei bambini della popolazione sia 50 cm (fonte Crescita bimbi in Italia)

IL risultato del test è il seguente:

One Sample t-test

data:  dataset$Lunghezza

t = -10.069, df = 2497, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 500

95 percent confidence interval:
 493.6628 495.7287

sample estimates:
mean of x 
 494.6958

Con un valore di p-value < 2.2e-16 (ovvero prossimo a 0), essendo minore di alfa (0.05), rifiutiamo l’ipotesi nulla, ovvero che la media del campione sia significativamente uguale a quella della popolazione.

Ciò è evidenziato anche dal 95% di confidenza dell’intervallo di accettazione dell’ipotesi, pari a 493.6628 - 495.7287 per cui il valore di media della popolazione (500 mm) non ricade al suo interno.

Script in R dei test t

Di seguito viene riportato lo script in R dei test t sopra citati.

# Test t per il peso

media_peso_popolazione <- 3300

test_peso <- t.test(dataset$Peso, mu = media_peso_popolazione)

# Estrazione del p-valore e della statistica t
p_value <- test_peso$p.value
t_stat <- as.numeric(test_peso$statistic)

# Livello di significatività
alpha <- 0.05

print("Risultati del test t per il peso:")
print(test_peso)

# Decisione
if (p_value < alpha) {
  print("Rifiuta l'ipotesi nulla H0")
} else {
  print("Non rifiutare l'ipotesi nulla H0")
}


# Test t per la lunghezza

media_lunghezza_popolazione <- 500

test_lunghezza <- t.test(dataset$Lunghezza, mu = media_lunghezza_popolazione)

# Estrazione del p-valore e della statistica t
p_value <- test_lunghezza$p.value
t_stat <- as.numeric(test_lunghezza$statistic)

# Stampa i risultati dei test
print("Risultati del test t per la lunghezza:")
print(test_lunghezza)

# Decisione
if (p_value < alpha) {
  print("Rifiuta l'ipotesi nulla H0")
} else {
  print("Non rifiutare l'ipotesi nulla H0")
}

5) Per le stesse variabili, o per altre per le quali ha senso farlo, verifica differenze significative tra i due sessi

Allo scopo verranno effettuati 3 test t di Student, ciascuno per verificare differenze significative tra i due sessi per le variabili Peso, Lunghezza e Cranio, con il seguente sistema delle ipotesi:

• H0 (ipotesi nulla) = H0 ​: μ-μ0=0

• H1 (ipotesi alternativa)​ = H1 ​: μ-μ0 != 0​

Dove:

• μ è la media del gruppo delle femmine (F)

• μ0​ è la media del gruppo delle maschi (M)

Fissiamo, per tutti i test t, un livello di significativà α = 0.05.

Di seguito lo script in R

#test t per Peso, Lunghezza, Cranio
test_peso <- t.test(Peso ~ Sesso, data = dataset)
test_lunghezza <- t.test(Lunghezza ~ Sesso, data = dataset)
test_cranio <- t.test(Peso ~ Sesso, data = dataset)

# Mostra i risultati
print(test_peso)
print(test_lunghezza)
print(test_cranio)

Di seguito sono mostrati i risultati dei test:

    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 



    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 



    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 

Per tutti i test t il valore p-value < α (fissato a 0.05), per cui si rifiuta l’ipotesi nulla e concludere che esiste una differenza significativa tra le medie del gruppo delle femmine ed il gruppo dei maschi.

6) Si vocifera che in alcuni ospedali si facciano più parti cesarei, sai verificare questa ipotesi?

Per verificare l’ipotesi che in alcuni ospedali si facciano più parti cesarei rispetto ad altri, si eseguirà un test del chi-quadrato (χ2) di associazione o indipendenza tra la variabile “Tipo.parto” e la variabile “Ospedale”. Si andrà a verificare in questo modo se esiste una relazione significativa tra due variabili categoriali.

Di seguito il sistema di ipotesi:

• H0 (ipotesi nulla): La proporzione di parti cesarei è la stessa in tutti gli ospedali (non c’è associazione tra l’ospedale e il tipo di parto).
• H1 (ipotesi alternativa): La proporzione di parti cesarei varia tra gli ospedali (c’è un’associazione tra l’ospedale e il tipo di parto).

Di seguito lo script in R

# Creazione della tabella di contingenza
tabella_contingenza <- table(dataset$Ospedale, dataset$Tipo.parto)

# Visualizza la tabella di contingenza
print("Tabella di contingenza tra Ospedale e Tipo di parto:")
print(tabella_contingenza)

# Esegui il test del chi-quadrato
test_chi_quadrato <- chisq.test(tabella_contingenza)

# Stampa i risultati del test
print("Risultati del test del chi-quadrato:")
print(test_chi_quadrato)

# Decisione
if (test_chi_quadrato$p.value < alpha) {
  print("Rifiuta l'ipotesi nulla H0")
} else {
  print("Non rifiutare l'ipotesi nulla H0")
}

Il test χ2 è stato eseguito su una tabella di contingenza che incrocia le variabili “Ospedale” e “Tipo.parto”.

Viene riportato di seguito l’ouptut dello script che stampa il risultato del test, che include il valore del chi-quadrato, i gradi di libertà e il valore p.

    Pearson's Chi-squared test

data:  tabella_contingenza
X-squared = 1.083, df = 2, p-value = 0.5819

Con un valore di p-value = 0.5819, essendo minore del valore di alfa (0.05), non si può rifiutare l’ipotesi nulla, suggerendo che non c’è evidenza sufficiente per affermare che la proporzione di parti cesarei varia significativamente tra gli ospedali.

7) Crea un modello di regressione lineare multipla con tutte le variabili e commenta i coefficienti e il risultato ottenuto

Di seguito lo script in R di creazione del modello:

# Modello completo (tutte le variabili)
modello_completo <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
                         Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dataset)

# Visualizza il sommario del modello
summary(modello_completo)

Il modello di regressione lineare multipla viene creato utilizzando la funzione lm(), dove il peso è la variabile dipendente e tutte le altre variabili sono indipendenti.

La funzione summary fornisce una panoramica completa del modello, inclusi i coefficienti stimati, i valori p, l’R-squared, e altri indicatori.

Di seguito l’output del modello:

Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
    Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1123.26  -181.53   -14.45   161.05  2611.89 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6735.7960   141.4790 -47.610  < 2e-16 ***
Anni.madre        0.8018     1.1467   0.699   0.4845    
N.gravidanze     11.3812     4.6686   2.438   0.0148 *  
FumatriciSI     -30.2741    27.5492  -1.099   0.2719    
Gestazione       32.5773     3.8208   8.526  < 2e-16 ***
Lunghezza        10.2922     0.3009  34.207  < 2e-16 ***
Cranio           10.4722     0.4263  24.567  < 2e-16 ***
Tipo.partoNat    29.6335    12.0905   2.451   0.0143 *  
Ospedaleosp2    -11.0912    13.4471  -0.825   0.4096    
Ospedaleosp3     28.2495    13.5054   2.092   0.0366 *  
SessoM           77.5723    11.1865   6.934 5.18e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 274 on 2487 degrees of freedom
Multiple R-squared:  0.7289,    Adjusted R-squared:  0.7278 
F-statistic: 668.7 on 10 and 2487 DF,  p-value: < 2.2e-16

• Coefficiente: Ogni variabile indipendente ha un coefficiente associato, che rappresenta l’effetto stimato di quella variabile sul peso del neonato, tenendo costanti tutte le altre variabili.

  • Un coefficiente positivo indica che un aumento della variabile indipendente è associato a un aumento del peso del neonato.

  • Un coefficiente negativo indica che un aumento della variabile indipendente è associato a una diminuzione del peso del neonato.

• Valore p: Il valore p per ciascun coefficiente indica se l’effetto della variabile indipendente è statisticamente significativo.

  • Un valore p inferiore a 0,05 (tipicamente) suggerisce che la variabile è significativamente associata al peso del neonato. Nel modello completo le variabili N.gravidanze, Gestazione, lunghezza, cranio, Tipo.partoNat, ospedaleosp3 e sessoM vanno in questa direzione.

  • Un valore p maggiore di 0,05 suggerisce che non c’è evidenza sufficiente per affermare che la variabile abbia un effetto significativo. Le variabili sono Anni.madre, FumatriciSI, ospedaleosp2

• R-squared: Questo valore rappresenta la proporzione della varianza nella variabile dipendente (peso) spiegata dal modello. Un valore più alto indica che il modello spiega meglio la variabilità dei dati. Nel modello di riferimento il valore è 0.7289

• Adjusted R-squared: È una versione modificata dell’R-squared che tiene conto del numero di predittori nel modello. È utile per confrontare modelli con un diverso numero di predittori. Nel modello completo il valore è 0.7278

8) Cerca il modello “migliore”, utilizzando tutti i criteri di selezione che conosci e spiegali.

Definizione del modello migliore tramite procedura stepwise manuale

Al modello completo è stato applicato in modo manuale la procedura stepwise al fine di valutare la possibilità di semplificarlo, togliendo le variabili poco significative.

Togliendo uno alla volta le variabili con Valore p > 0,05 e si è osservato un valore di Adjusted R-squared pressocchè simile a quello del modello completo, ma con il vantaggio di aver semplificato il modello stesso.

Di seguito sono riportati gli output dei vari modelli semplificati:

Modello completo senza la variabile Fumatrice

Si è deciso di togliere la variabile Fumatrice avendo un p-value molto alto(0.2719) anche se sono conclamati gli effetti del fumo sulla salute e sul peso alla nascita dei neonati.

La poca correlazione rispetto al peso potrebbe spiegarsi dal fatto che nel campione solo meno del 5% delle donne dichiara di fumare.

Tuttavia i valori della variabile potrebbero non rappresentare a pieno la popolazione, ragion per cui il modello potrebbe includere in futuro anche questa variabile.

Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
    Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1122.63  -181.96   -14.91   161.39  2615.07 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6735.4444   141.4845 -47.606  < 2e-16 ***
Anni.madre        0.8118     1.1467   0.708   0.4791    
N.gravidanze     11.1201     4.6627   2.385   0.0172 *  
Gestazione       32.3210     3.8138   8.475  < 2e-16 ***
Lunghezza        10.3064     0.3006  34.285  < 2e-16 ***
Cranio           10.4766     0.4263  24.577  < 2e-16 ***
Tipo.partoNat    29.3770    12.0888   2.430   0.0152 *  
Ospedaleosp2    -11.0363    13.4475  -0.821   0.4119    
Ospedaleosp3     28.5194    13.5038   2.112   0.0348 *  
SessoM           77.3928    11.1858   6.919 5.77e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 274 on 2488 degrees of freedom
Multiple R-squared:  0.7288,    Adjusted R-squared:  0.7278 
F-statistic: 742.8 on 9 and 2488 DF,  p-value: < 2.2e-16

• Multiple R-squared modello completo: 0.7289

• Multiple R-squared modello semplificato: 0.7288

• R-squared modello completo: 0.7278

• R-squared modello semplificato: 0.7278

Modello completo senza le variabili Fumatrice e Anni.madre

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Tipo.parto + Ospedale + Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1113.07  -181.71   -16.66   161.08  2619.57 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6707.9252   136.0257 -49.314  < 2e-16 ***
N.gravidanze     12.3360     4.3344   2.846  0.00446 ** 
Gestazione       32.0386     3.7925   8.448  < 2e-16 ***
Lunghezza        10.3059     0.3006  34.286  < 2e-16 ***
Cranio           10.4920     0.4257  24.648  < 2e-16 ***
Tipo.partoNat    29.4080    12.0875   2.433  0.01505 *  
Ospedaleosp2    -10.8939    13.4447  -0.810  0.41786    
Ospedaleosp3     28.7917    13.4969   2.133  0.03301 *  
SessoM           77.4657    11.1842   6.926 5.48e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 274 on 2489 degrees of freedom
Multiple R-squared:  0.7287,    Adjusted R-squared:  0.7278 
F-statistic: 835.7 on 8 and 2489 DF,  p-value: < 2.2e-16

• Multiple R-squared modello completo: 0.7289

• Multiple R-squared modello semplificato: 0.7287

• R-squared modello completo: 0.7278

• R-squared modello semplificato: 0.7278

Modello completo senza senza le variabili Fumatrici, Anni.madre e Ospedale

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Tipo.parto + Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1129.14  -181.97   -16.26   160.95  2638.18 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6708.0171   136.0715 -49.298  < 2e-16 ***
N.gravidanze     12.7356     4.3385   2.935  0.00336 ** 
Gestazione       32.3253     3.7969   8.514  < 2e-16 ***
Lunghezza        10.2833     0.3009  34.177  < 2e-16 ***
Cranio           10.5063     0.4263  24.648  < 2e-16 ***
Tipo.partoNat    30.1601    12.1027   2.492  0.01277 *  
SessoM           77.9171    11.1994   6.957 4.42e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 274.4 on 2491 degrees of freedom
Multiple R-squared:  0.7277,    Adjusted R-squared:  0.727 
F-statistic:  1109 on 6 and 2491 DF,  p-value: < 2.2e-16

• Multiple R-squared modello completo: 0.7289

• Multiple R-squared modello semplificato: 0.7277

• R-squared modello completo: 0.7278

• R-squared modello semplificato: 0.727

Modello completo senza senza le variabili Fumatrici, Anni.madre e Ospedale, Tipo.parto

Si è deciso di eliminare dal modello la variabile Tipo.partoNat in qunto si ritiene non avere un effetto significativo.

Come si evince difatti dal summary sotto riportato, il valore di Adjusted R-squared del nuovo modello è rimasto pressochè identico al valore del modello completo.

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1149.37  -180.98   -15.57   163.69  2639.09 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
Cranio          10.5410     0.4265  24.717  < 2e-16 ***
SessoM          77.9807    11.2111   6.956 4.47e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 274.7 on 2492 degrees of freedom
Multiple R-squared:  0.727, Adjusted R-squared:  0.7265 
F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

• Multiple R-squared modello completo: 0.7289

• Multiple R-squared modello semplificato: 0.727

• R-squared modello completo: 0.7278

• R-squared modello semplificato: 0.7265

Definizione del modello migliore tramite la procedura stepwise utilizzando le funzioni di R

Di seguito lo script in R che effettua le seguenti operazioni:

• genera 3 nuovi modelli tramite la procedura stepwise

• confronta il modello completo, il modello personale ed i 3 modelli di cui sopra, usando i criteri AIC e BIC

# Forward
modello_forward <- MASS::stepAIC(modello_completo,
                                 direction = "forward",
                                 k=log(n))

# Backward
modello_backward <- MASS::stepAIC(modello_completo,
                                  direction = "backward",
                                  k=log(n))

# Mixed
modello_stepwise <- MASS::stepAIC(modello_completo,
                                  direction = "both",
                                  k=log(n))

# Confronto dei modelli
summary(modello_completo)
summary(modello_forward)
summary(modello_backward)
summary(modello_stepwise)
summary(modello_manuale)

# Calcolo AIC e BIC per i modelli
aic_completo <- AIC(modello_completo)
aic_forward <- AIC(modello_forward)
aic_backward <- AIC(modello_backward)
aic_stepwise <- AIC(modello_stepwise)
aic_modello_test <- AIC(modello_manuale)

bic_completo <- BIC(modello_completo)
bic_forward <- BIC(modello_forward)
bic_backward <- BIC(modello_backward)
bic_stepwise <- BIC(modello_stepwise)
bic_modello_test <- BIC(modello_manuale)

# Stampa i risultati di AIC e BIC
cat("AIC:\n")
cat("Completo:", aic_completo, "\nForward:", aic_forward, "\nBackward:", aic_backward, "\nStepwise:", aic_stepwise, "\nModello manuale:", aic_modello_test, "\n")

cat("BIC:\n")
cat("Completo:", bic_completo, "\nForward:", bic_forward, "\nBackward:", bic_backward, "\nStepwise:", bic_stepwise, "\nModello manuale:", bic_modello_test, "\n")

Viene di seguito riportato il risultato dei confronti:

• AIC:

  • Completo: 35145.57

  • Forward: 35145.57

  • Backward: 35152.89

  • Stepwise: 35152.89

  • Modello personale: 35152.89

• BIC:

  • Completo: 35215.45

  • Forward: 35215.45

  • Backward: 35193.65

  • Stepwise: 35193.65

  • Modello personale: 35193.65

Il criterio di selezione BIC personalmente è da preferirsi in quanto tende a penalizzare i modelli più complessi.

Secondo questo criterio, i modelli migliori sono il modello manuale, il backward e lo Stepwise both (sono identici), con un valore di 35193.65.

9) Si potrebbero considerare interazioni o effetti non lineari?

La relazione tra le variabili Peso e Lunghezza e Peso e Cranio suggerisce una relazione non lineare per cui si andrà ad aggiungere al modello di riferimento un termine quadratico, per poi osservarne i risultati.

# Modello di regressione lineare con interazioni e termini quadratici
modello_quadratico <- lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + I(Cranio^2) + Sesso, data = dataset)

# Sommario del modello
summary(modello_quadratico)

aic_modello_quadratico <- AIC(modello_quadratico)
bic_modello_quadratico <- BIC(modello_quadratico)

print(aic_modello_quadratico)
print(bic_modello_quadratico)

Viene mostrato di seguito l’output del modello:

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + 
    I(Cranio^2) + Sesso, data = dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-1157.94  -179.42   -12.32   165.84  2366.68 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -2.653e+03  1.176e+02 -22.570  < 2e-16 ***
N.gravidanze    1.313e+01  4.289e+00   3.062  0.00222 ** 
Gestazione      3.651e+01  3.695e+00   9.882  < 2e-16 ***
I(Lunghezza^2)  1.083e-02  3.062e-04  35.360  < 2e-16 ***
I(Cranio^2)     1.560e-02  6.218e-04  25.081  < 2e-16 ***
SessoM          7.338e+01  1.108e+01   6.620 4.38e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 271.4 on 2492 degrees of freedom
Multiple R-squared:  0.7336,    Adjusted R-squared:  0.7331 
F-statistic:  1372 on 5 and 2492 DF,  p-value: < 2.2e-16

Il valori di Adjusted R-squared: 0.7331 e BIC = 35132.06 sono i migliori di tutti quelli osservati tra i vari modelli.

---

Sono riportati di seguito i grafici Component + Residual Plots del miglior modello con tutti i termini lineari ed il modello con termini quadratici, al fine di visualizzare le relazioni tra la variabile Peso rispetto le altre.

Grafico Component + Residual Plot del miglior modello trovato con tutti i termini lineari

Grafico Component + Residual Plot del modello con termini quadratici per le variabili Lunghezza e Cranio

---

La differenza tra i grafici Lunghezza e Cranio dei due modelli di regressione evidenzia che i termini quadratici spieghino meglio la relazione tra le variabili Peso-Lunghezza e Peso-Cranio.

Verrà pertanto preso come modello di riferimento per la diagnostica dei residui e per le previsioni.

10) Effettua una diagnostica approfondita dei residui del modello e di potenziali valori influenti. Se ne trovi prova a verificare la loro effettiva influenza

Grafici standard per i residui ed il modello

Di seguito i 4 grafici diagnostici standard per i residui e il modello

1. Residuals vs Fitted (Residui vs Valori Predetti)

• Verifica l’ipotesi di omoschedasticità (varianza costante dei residui) e la linearità del modello.

• Si osserva un pattern abbastanza casuale per cui non vi sono violazioni evidenti delle assunzioni.

• Tuttavia si rilevano alcuni punti fuori dalla media (come il punto identificato come 1551), che potrebbe indicare osservazioni influenti o possibili outlier.

2. Q-Q Plot dei Residui

• Verifica l’ipotesi di normalità dei residui.

• I punti per la maggior parte si allineano bene lungo la linea retta diagonale, ma agli estremi si discostano dalla linea, soprattutto a destra. Questo suggerisce che ci sono alcuni outlier o valori estremi che non seguono la distribuzione normale.

3. Scale-Location (o Spread-Location) Plot

• Valuta l’eteroschedasticità (varianza non costante).

• La dispersione appare abbastanza casuale attorno alla linea orizzontale per la maggior parte dei punti, ma ci sono alcuni outlier che potrebbero indicare un aumento della varianza in alcune regioni (come nel caso dell’osservazione 1551 e 155).

4. Residuals vs Leverage (Residui vs Leverage) con la Distanza di Cook

• Identifica osservazioni influenti e outlier. Il leverage misura quanto un’osservazione è distante dai valori medi delle variabili indipendenti. La distanza di Cook (curve grigie) misura l’influenza combinata del leverage e dei residui.

• Le osservazioni che hanno sia un leverage elevato che una grande distanza di Cook possono avere un’influenza significativa sul modello e potrebbero essere considerate influenti. Questi punti sono potenzialmente distorsivi.

• Alcuni punti come 1551, 310, e 1780 sembrano avere leverage elevato o una distanza di Cook significativa. Potrebbero pertanto essere osservazioni influenti

In definitiva, I grafici mostrano un modello che potrebbe essere influenzato da alcuni outlier, in particolare le osservazioni identificate come 1551 e altre nei vari grafici.

Tuttavia, per la maggior parte dei dati, il modello quadrativo sembra soddisfare le assunzioni, con una leggera deviazione dalla normalità e alcuni possibili outlier o punti influenti.

Analisi leverage

Viene mostrato, di seguito, il grafico che mostra i valori di leverage per ogni osservazione, che misura quanto un’ osservazione sia distante dal centro del campione.

Le osservazioni con leverage alto sono quelle che hanno un valore di leverage maggiore di 2 * mean(hat_values), ovvero tutte le osservazioni che si trovano al di sopra della linea rossa.

Distanza di Cook

Il grafico seguente calcola la distanza di Cook per ogni osservazione, che misura l’influenza di ogni osservazione sul modello complessivo.

Il valore soglia è calcolato tramite la formula (4 * <numero osservazioni>) ovvero tutte le osservazioni che si trovano al di sopra della linea rossa.

Considerazioni

E’ stato identificato il set delle osservazioni influenti del modello di regressione tramite le influenze basate sulla distanza di Cook e le osservazioni con lavarage alto.

E’ stata effettuata un’analisi visiva approfondita su questi punti, con particolare riferimento alle osservazioni presenti nei due gruppi di cui sopra.

Non sono state rilevate anomalie sui dati, pertanto tali osservazioni sono da considerarsi corrette.

Viene riportato di seguito lo script in R.

# 1. Analisi dei residui

# Grafico dei residui standardizzati vs valori predetti
par(mfrow = c(2, 2))  # Disposizione 2x2 per più grafici
plot(modello_quadratico)

#test di omoscedasticità
lmtest::bptest(modello_quadratico)

#test di verifica autocorrelazione dei residui
lmtest::dwtest(modello_quadratico)

#test di normalità dei residui
shapiro.test(modello_quadratico$residuals)

#stima densita dei residui
plot(density(residuals(modello_quadratico)))

#leverage
# Identifica punti con alta leva
hatvalues <- hatvalues(modello_quadratico)
plot(hatvalues, main="Valori di leva", ylab="Valori di leva", xlab="Indice")
abline(h = 2*mean(hatvalues), col="red")

leverage_threshold <- 2 * mean(hat_values)
high_leverage_points <- which(hat_values > leverage_threshold)
high_leverage_points


#outliers
plot(rstudent(modello_quadratico))
abline(h=c(-2,2))
car::outlierTest(modello_quadratico)


# Identificazione di valori influenti: Distanza di Cook
cook_distances <- cooks.distance(modello_quadratico)
influential_points <- which(cook_distances > (4 / nrow(dataset)))
influential_points

# Visualizzazione dei valori influenti
plot(cook_distances, type="h", main="Distanza di Cook", ylab="Distanza di Cook")
abline(h = 4 / nrow(dataset), col = "red")


# Visualizzare leverage e distanza di Cook insieme
plot(hat_values, cook_distances, xlab = "Leverage", ylab = "Distanza di Cook", main = "Leverage vs Distanza di Cook")
abline(h = 4 / nrow(dataset), col = "red")
abline(v = leverage_threshold, col = "blue")

# Grafico di influenza
influencePlot(modello_quadratico, id.method="identify", main="Grafico di influenza", sub="La dimensione del cerchio è proporzionale alla distanza di Cook")

8) Fai la tua migliore previsione per il peso di una neonata, considerato che la madre è alla terza gravidanza e partorirà alla 39esima settimana. Niente misure dall’ecografia.

new_data <- data.frame(
  Anni.madre = mean(dataset$Anni.madre, na.rm = TRUE),  # età media della madre
  N.gravidanze = 3,
  Fumatrici = 0,  # Supponiamo che la madre non fumi
  Gestazione = 39,
  Lunghezza = mean(dataset$Lunghezza, na.rm = TRUE),  # lunghezza media del neonato
  Cranio = mean(dataset$Cranio, na.rm = TRUE),  # circonferenza cranica media
  Tipo.parto = factor("Nat", levels = levels(dataset$Tipo.parto)),  # parto naturale
  Ospedale = factor("osp1", levels = levels(dataset$Ospedale)),  # esempio di ospedale
  Sesso = factor("F", levels = levels(dataset$Sesso))  # sesso femminile
)

predicted_values = predict(modello_quadratico, newdata = new_data)

print(predicted_values)

Per le variabili Lunghezza e Cranio viene passata la media dei valori presenti nel campione.

La previsione del peso è 3263.211 grammi.

9) Cerca di creare qualche rappresentazione grafica che aiuti a visualizzare il modello. Se è il caso semplifica quest’ultimo!

Allo scopo verrà generato un 3d Scatter plot utilizzando le variabili Lunghezza - Gestazione e Peso.

Cliccare sull’immagine per visualizzare il grafico 3d interattivo.