Modello statistico per la previsione del peso neonatale

1. Raccolta dei dati e struttura del dataset

# Caricamento del dataset dal file CSV "neonati.csv"
df <- read.csv(
  "neonati.csv",           # Nome del file
  stringsAsFactors = TRUE  # Converte automaticamente le variabili testuali in fattori 
)

# Visualizzazione delle prime 10 righe del dataset in tabella, con formattazione kableExtra
kable(head(df, 10)) %>%  # Mostra le prime 10 osservazioni per un’anteprima del contenuto
  kable_styling(
    full_width = TRUE,                        # Tabella a larghezza piena
    bootstrap_options = c("striped", "hover") # Righe alternate a righe colorate (striped) e effetto hover
  )

Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F
32	0	40	3200	495	340	Nat	osp2	F
26	1	39	3100	480	345	Nat	osp3	F
25	0	40	3580	510	349	Nat	osp1	M
22	1	40	3670	500	335	Ces	osp2	F
23	0	41	3700	510	362	Ces	osp2	F

# Stampa una stringa HTML con le dimensioni del dataset in output "grezzo"
knitr::asis_output(
  paste0(
    "<h4><strong>Dimensioni del Dataset:</strong> ", 
    nrow(df),  # Numero di righe del dataset
    " righe x ", 
    ncol(df),  # Numero di colonne del dataset
    " colonne</h4>"
  )
)

Dimensioni del Dataset: 2500 righe x 10 colonne

# Crea una tabella con le informazioni base su ogni variabile del dataset
glimpse_info <- data.frame(
  # Colonna "Tipo": tipo di dato per ogni variabile
  Tipo = sapply(df, class),  
  
  # Colonna "Esempi": primi 3 valori per ogni variabile, uniti da virgole
  Esempi = sapply(df, function(x) paste(head(x, 3), collapse = ", ")) 
)

# Visualizza la tabella con kable con uno stile per una migliore leggibilità
kable(glimpse_info, caption = "**Struttura del dataset**") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover"))

**Struttura del dataset**
	Tipo	Esempi
Anni.madre	integer	26, 21, 34
N.gravidanze	integer	0, 2, 3
Fumatrici	integer	0, 0, 0
Gestazione	integer	42, 39, 38
Peso	integer	3380, 3150, 3640
Lunghezza	integer	490, 490, 500
Cranio	integer	325, 345, 375
Tipo.parto	factor	Nat, Nat, Nat
Ospedale	factor	osp3, osp1, osp2
Sesso	factor	M, F, M

Anni.madre

Età della madre (in anni)
Tipo: Quantitativa continua (tempo)

N.gravidanze

Quante gravidanze ha avuto la madre
teniamo conto che 0 indica che è alla prima gravidanza ma non ha mai partorito prima
Tipo: Quantitativa discreta in scala di rapporti

Fumatrici ⁽¹⁾

Se la mamma fuma (0=non fumatrice, 1=fumatrice)
Tipo: Qualitativa dicotomica

Gestazione

Durata della gravidanza (in settimane)
Tipo: Quantitativa continua in scala di rapporti

Peso => VARIABILE TARGET

Peso del neonato (in grammi)
Tipo: Quantitativa continua in scala di rapporti

Lunghezza

Lunghezza del neonato (in millimetri)
Tipo: Quantitativa continua in scala di rapporti

Cranio ⁽²⁾

Circonferenza del cranio (in millimetri) e non diametro del cranio
Tipo: Quantitativa continua in scala di rapporti

Tipo.parto ⁽³⁾

Tipo di parto (Nat = Naturale, Ces = Cesareo)
Tipo: Qualitativa nominale

Ospedale ⁽⁴⁾

Ospedale di nascita (osp1, osp2, osp3)
Tipo: Qualitativa nominale

Sesso

Sesso del neonato (M = Maschio, F = Femmina)
Tipo: Qualitativa dicotomica

2. Analisi e modellizzazione

Clicca per mostrare il codice

# formatta il risultato del p-value
format_p_value <- function(p) {
  if (is.na(p)) {
    return("NA")
  } else if (p < 2e-16) {
    return("&lt; 2e-16")  
  } else if (p < 0.001) {
    return(formatC(p, format = "e", digits = 2))
  } else {
    return(format(round(p, 4), nsmall = 4))
  }
}

analyze_quantitative <- function(data, var_name, show_density_line = TRUE) {
  # Estrae la variabile quantitativa dal dataset
  var <- data[[var_name]]
  
  # Rimuove i valori mancanti
  var <- na.omit(var)
  
  # Calcola statistiche descrittive
  min_val <- min(var, na.rm = TRUE)
  mean_val <- mean(var, na.rm = TRUE)
  max_val <- max(var, na.rm = TRUE)
  q1_val <- quantile(var, 0.25, na.rm = TRUE)
  median_val <- median(var, na.rm = TRUE)
  q3_val <- quantile(var, 0.75, na.rm = TRUE)
  iqr_val <- IQR(var, na.rm = TRUE)
  
  # Calcola i limiti IQR per l’identificazione degli outlier
  lower_bound <- q1_val - 1.5 * iqr_val
  upper_bound <- q3_val + 1.5 * iqr_val
  outliers <- var[var < lower_bound | var > upper_bound]
  
  # Calcola skewness e kurtosis (asimmetria e forma)
  skewness_val <- skewness(var, na.rm = TRUE)
  kurtosis_val <- kurtosis(var, na.rm = TRUE) - 3  # Centro su 0 per confronto con la normale

  # Esegue il test di normalità di Shapiro-Wilk
  shapiro_test <- shapiro.test(var)
  
  # Crea un dataframe con tutte le statistiche calcolate
  results <- data.frame(
    Misura = c("Min", "Media", "Max", "Q1", "Mediana", "Q3", "IQR", 
               "Limite Inf. IQR", "Limite Sup. IQR", "Skewness", "Kurtosis", 
               "Shapiro-Wilk W", "Shapiro-Wilk p-value", "Numero Outlier"),
    Valore = c(min_val, mean_val, max_val, q1_val, median_val, q3_val, iqr_val, 
               lower_bound, upper_bound, skewness_val, kurtosis_val, 
               shapiro_test$statistic, shapiro_test$p.value, length(outliers))
  )
  
  # Applica la funzione di formattazione solo al p-value del test di Shapiro
  shapiro_index <- which(results$Misura == "Shapiro-Wilk p-value")
  results$Valore[shapiro_index] <- format_p_value(as.numeric(results$Valore[shapiro_index]))

  # Palette colori blu per i grafici
  blue_palette <- c("#1f78b4", "#a6cee3", "#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef")

  # Grafico della distribuzione: istogramma (+ densità opzionale)
  dist_plot <- ggplot(data, aes(x = .data[[var_name]])) +
    geom_histogram(aes(y = ..density..), bins = 30, fill = blue_palette[1], color = "black", alpha = 0.7) +
    labs(title = paste("Distribuzione di", var_name), x = var_name, y = "Densità") +
    theme_minimal()

  if (show_density_line) {
    dist_plot <- dist_plot + 
      geom_density(aes(color = "Dati Osservati"), linewidth = 1, color = blue_palette[3]) +
      theme(legend.position = "bottom") +
      guides(color = guide_legend(title = NULL))
  }

  # Boxplot per la variabile
  box_plot <- ggplot(data, aes(y = .data[[var_name]])) +
    geom_boxplot(fill = blue_palette[2], alpha = 0.7, color = blue_palette[4], outlier.color = "red", outlier.shape = 16) +
    labs(title = paste("Boxplot di", var_name), y = var_name) +
    theme_minimal()
  
  # Grafico jitter per evidenziare gli outlier
  outliers_plot <- ggplot(data, aes(x = 1, y = .data[[var_name]])) +
    geom_jitter(width = 0.1, color = blue_palette[6]) +
    geom_hline(yintercept = lower_bound, linetype = "dashed", color = "red") +
    geom_hline(yintercept = upper_bound, linetype = "dashed", color = "red") +
    labs(title = paste("Outlier di", var_name), x = "", y = var_name) +
    theme_minimal()
  
  # Crea un sotto-dataset con gli outlier
  outlier_df <- data[data[[var_name]] < lower_bound | data[[var_name]] > upper_bound, ]
  
  # Stampa la tabella delle statistiche
  print(
    kbl(results, format = "html", escape = FALSE, align = "c", caption = paste("Statistiche di", var_name)) %>%
      kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))
  )
  
  # Visualizza i grafici generati
  print(dist_plot)
  print(box_plot)
  print(outliers_plot)
  
  # Se ci sono outlier, stampa i 5 più bassi e i 5 più alti
  if (nrow(outlier_df) > 0) {
    extreme_outliers <- rbind(
      head(outlier_df[order(outlier_df[[var_name]]), ], 5),  # Primi 5 minori
      tail(outlier_df[order(outlier_df[[var_name]]), ], 5)   # Ultimi 5 maggiori
    )
    
    # Trova l'indice della colonna da evidenziare
    col_index <- which(colnames(extreme_outliers) == var_name)
    
    # Stampa tabella evidenziando la variabile studiata
    print(
      kbl(extreme_outliers, format = "html", caption = paste("Outlier estremi di", var_name)) %>%
        kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed")) %>%
        column_spec(col_index + 1, background = "#a6cee3", bold = TRUE)
    )
    
  } else {
    print("Nessun outlier rilevato.")
  }

  # Ritorna la lista con i dati degli outlier e le statistiche
  return(list(
    Outliers = outlier_df,
    Test_Statistici = results
  ))
}

analyze_qualitative <- function(data, var_name) {
  var <- na.omit(data[[var_name]])
  
  # Creazione della tabella di frequenza con percentuali
  freq_table <- as.data.frame(table(var))
  colnames(freq_table) <- c("Categoria", "Frequenza")
  freq_table$Percentuale <- round(100 * freq_table$Frequenza / sum(freq_table$Frequenza), 2)
  
  # Test del Chi-quadro per verificare la distribuzione uniforme
  frequenze_reali <- table(var)
  chi_test <- chisq.test(frequenze_reali)
  
  # Se binaria, test binomiale per equilibrio
  if (nrow(freq_table) == 2) {
    binom_test <- binom.test(freq_table$Frequenza[1], sum(freq_table$Frequenza))
    binom_pvalue <- binom_test$p.value
  } else {
    binom_pvalue <- NA
  }
  
  # Gradazioni di blu per i grafici
  blue_palette <- c("#1f78b4","#1f78b4","#1f78b4","#1f78b4")
  
  # Grafico a barre con percentuali
  bar_plot <- ggplot(freq_table, aes(x = Categoria, y = Frequenza, fill = Categoria)) +
    geom_bar(stat = "identity", alpha = 0.8) +
    geom_text(aes(label = paste0(Percentuale, "%")), vjust = -0.3, size = 4, color = "black") +
    scale_fill_manual(values = blue_palette[1:length(unique(freq_table$Categoria))]) +  # Usa i colori blu
    labs(title = paste("Distribuzione di", var_name), x = var_name, y = "Frequenza") +
    theme_minimal() +
    theme(legend.position = "none")
  
  # Tabella dei risultati statistici
  test_results <- data.frame(
    Test = c("Chi-quadro", "Test binomiale"),
    P_value = c(chi_test$p.value, binom_pvalue)
  )
  
  # Stampa della tabella di frequenza
  print(
    kbl(freq_table, format = "html", caption = paste("Tabella delle Frequenze di", var_name)) %>%
      kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))
  )
  
  # Stampa della tabella dei test statistici
  print(
    kbl(test_results, format = "html", escape = FALSE, caption = "Risultati dei Test Statistici") %>%
      kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))
  )
  
  # Mostra i grafici
  print(bar_plot)
  
  return(list(Frequenze = freq_table, Test_Statistici = test_results))
}

Analisi preliminare

Analisi di “Anni.madre”

anni_madre_df <- analyze_quantitative(df, "Anni.madre")

Statistiche di Anni.madre
Misura	Valore
Min	0
Media	28.164
Max	46
Q1	25
Mediana	28
Q3	32
IQR	7
Limite Inf. IQR	14.5
Limite Sup. IQR	42.5
Skewness	0.0428115004569053
Kurtosis	0.3804165222879
Shapiro-Wilk W	0.993080970998317
Shapiro-Wilk p-value	1.64e-09
Numero Outlier	13

Outlier estremi di Anni.madre
	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1380	0	0	0	39	3060	490	330	Nat	osp3	M
1152	1	1	0	41	3250	490	350	Nat	osp2	F
138	13	0	0	38	2760	470	325	Nat	osp2	F
1075	14	1	0	39	3510	490	365	Nat	osp2	M
1532	14	0	0	39	3550	500	355	Ces	osp1	M
335	44	0	1	38	3150	465	335	Nat	osp3	F
2026	44	0	0	40	3050	505	345	Nat	osp1	M
2098	44	5	0	38	3280	475	346	Ces	osp1	F
205	45	2	0	38	3850	505	384	Nat	osp3	M
1106	46	5	0	36	2710	470	347	Nat	osp2	M

La variabile segue una distribuzione leggermenete assimmetrica con code meno pronucniate rispetto ad una distribuzione normale. Il testo di Shapiro-Wilk indica che la distribuzione non è normale, ma dall’istrogramma della distribuzione la forma ricorda la distribuzione normale percui potremmo per ora considerare la distribuzione come normale e fare ulteriori analisi successivamente per vedere l’influenza sul modello.

Analizzando gli outlier possiamo vedere valori molto bassi di 0, 1, 13, 14 anni e valori più alti di 44, 45, 46 anni.

I valori 0 e 1 sono valori chiaramente anomali dovuto quasi sicuramente a errori di inserimento.

I valori di 13 e 14 anni ⁽⁵⁾, sono valori plausibili anche se gravidanze a queste giovane età posso comportare maggiori rischi a causa di complicanze.

I valori di 44, 45 e 46 anni sono valori plausibili ⁽⁶⁾ e sta diventando sempre più anche grazie ai progressi della medicina riproduttiva.

Sulla base di questo decido di sostituire i valori 0 e 1 con la mediana delle età delle madri per evitare di elinminarli e perdere informazione e di tenere gli altri valori in quando plausibili e necessari per mantenere significativo il dataset.

Analisi di “N.gravidanze”

gravidanze_df <- analyze_quantitative(df, "N.gravidanze", show_density_line = FALSE)

Statistiche di N.gravidanze
Misura	Valore
Min	0
Media	0.9812
Max	12
Q1	0
Mediana	1
Q3	1
IQR	1
Limite Inf. IQR	-1.5
Limite Sup. IQR	2.5
Skewness	2.51425407710076
Kurtosis	10.9894063615823
Shapiro-Wilk W	0.721249287971459
Shapiro-Wilk p-value	< 2e-16
Numero Outlier	246

Outlier estremi di N.gravidanze
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
3	34	3	38	3640	500	375	Nat	osp2	M
15	33	3	34	2400	470	298	Ces	osp3	M
18	32	3	40	3030	490	335	Ces	osp3	M
40	34	3	38	3050	460	336	Nat	osp3	F
47	29	3	39	2680	450	346	Nat	osp3	F
2221	35	10	39	2950	495	335	Nat	osp1	F
2422	33	10	40	3090	485	353	Nat	osp3	M
2471	34	10	38	2880	470	345	Ces	osp2	M
1130	33	11	43	3400	475	360	Nat	osp1	M
1219	38	12	39	3350	490	344	Nat	osp2	M

La distribuzione è molto assimmetrica a desta con una coda lunga e con molti valori outlier, questo significa che meno di 3 gravidanze e un dato comune, invece più di 3 gravidanze sono casi piuttosto rari che rappresntano circa il 10% dei casi totali.

Analisi di “Fumatrici”

fumatrici_df <- analyze_qualitative(df, "Fumatrici")

Tabella delle Frequenze di Fumatrici
Categoria	Frequenza	Percentuale
0	2396	95.84
1	104	4.16

Risultati dei Test Statistici
Test	P_value
Chi-quadro	0
Test binomiale	0

Nel nostro campione le madri Fumatrici rappresentano circa il 4% rendendo questa distribuzione molto sbilanciata e potrebbe non essere correttamente rappresentativa nel modello.

Ma considerando che la media sulla popolazione sarebbe intorno al 10% ⁽¹⁾ questo permetterebbe di ottenere risultati più significativi e affidabili per il modello, indicando che probabilmente è il nostro campione a non essere perfettamente bilanciato.

Analisi di “Gestazione”

gestazione_df <- analyze_quantitative(df, "Gestazione", show_density_line = FALSE)

Statistiche di Gestazione
Misura	Valore
Min	25
Media	38.9804
Max	43
Q1	38
Mediana	39
Q3	40
IQR	2
Limite Inf. IQR	35
Limite Sup. IQR	43
Skewness	-2.06531325911714
Kurtosis	8.2581495549494
Shapiro-Wilk W	0.833276929649533
Shapiro-Wilk p-value	< 2e-16
Numero Outlier	67

Outlier estremi di Gestazione
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1780	25	2	25	900	325	253	Nat	osp3	F
2452	28	0	26	930	345	245	Ces	osp3	F
2120	32	0	27	1140	370	267	Nat	osp3	F
2437	28	1	27	980	320	265	Nat	osp1	M
310	40	3	28	1560	420	379	Nat	osp3	F
1311	40	6	34	1840	430	305	Nat	osp1	F
1400	22	0	34	2590	485	314	Nat	osp2	M
1893	22	1	34	3030	470	312	Nat	osp2	F
1977	39	4	34	2970	480	350	Ces	osp2	F
2257	40	0	34	1580	400	300	Nat	osp2	F

Le gravidanze durano dalle 25 alle 43 settimane, con un media di 39 settimane che si attesta vicino alla durata tipica di una gravidanza che è di 40 settimane. Sono presenti diversi outlier sopratutto sotto le 35 settimane che identificano i parti prematuri.

Analisi di “Peso”

peso_df <- analyze_quantitative(df, "Peso")

Statistiche di Peso
Misura	Valore
Min	830
Media	3284.0808
Max	4930
Q1	2990
Mediana	3300
Q3	3620
IQR	630
Limite Inf. IQR	2045
Limite Sup. IQR	4565
Skewness	-0.647030832752079
Kurtosis	2.03153184511532
Shapiro-Wilk W	0.970657098407693
Shapiro-Wilk p-value	< 2e-16
Numero Outlier	69

Outlier estremi di Peso
	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
928	25	0	0	28	830	310	254	Nat	osp1	F
1780	25	2	0	25	900	325	253	Nat	osp3	F
2175	37	8	0	28	930	355	235	Nat	osp1	F
2452	28	0	0	26	930	345	245	Ces	osp3	F
2437	28	1	0	27	980	320	265	Nat	osp1	M
2392	28	1	0	40	4720	540	355	Nat	osp3	M
1639	39	3	0	40	4760	550	365	Nat	osp2	F
1433	31	4	0	41	4810	530	364	Ces	osp3	M
1306	23	0	0	41	4900	510	352	Nat	osp2	F
1920	26	0	1	39	4930	550	350	Ces	osp2	F

Siccome il peso medio e il peso mediano hanno valori abbastanza simili in generale la maggior parte dei neonati da un peso introno a questi valori ma con alcuni neonati che hanno peso significativamente più bassi che potrebbe dipendere da nascite premature o condiczione che influenzano la crescita del feto.

Analisi di “Lunghezza”

lunghezza_df <- analyze_quantitative(df, "Lunghezza")

Statistiche di Lunghezza
Misura	Valore
Min	310
Media	494.692
Max	565
Q1	480
Mediana	500
Q3	510
IQR	30
Limite Inf. IQR	435
Limite Sup. IQR	555
Skewness	-1.51469914804306
Kurtosis	6.48717388620972
Shapiro-Wilk W	0.90940573211939
Shapiro-Wilk p-value	< 2e-16
Numero Outlier	59

Outlier estremi di Lunghezza
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
928	25	0	28	830	310	254	Nat	osp1	F
1551	35	1	38	4370	315	374	Nat	osp3	F
2437	28	1	27	980	320	265	Nat	osp1	M
1780	25	2	25	900	325	253	Nat	osp3	F
1619	31	0	31	990	340	278	Ces	osp2	F
2359	25	6	33	2230	430	313	Nat	osp3	F
2458	31	0	31	1730	430	300	Nat	osp3	F
471	29	0	40	3680	560	346	Nat	osp2	M
1513	25	0	40	4620	560	367	Nat	osp2	M
2391	36	1	41	4400	565	366	Nat	osp1	F

La distribuzione è assimetrica con prevalenza di valori più bassi. I valori più bassi si possono riferire ai parti prematuri, mentre quelli più alti potrebbero essere errori di muisurazione o casi particolari rari.

Analisi di “Cranio”

cranio_df <- analyze_quantitative(df, "Cranio")

Statistiche di Cranio
Misura	Valore
Min	235
Media	340.0292
Max	390
Q1	330
Mediana	340
Q3	350
IQR	20
Limite Inf. IQR	300
Limite Sup. IQR	380
Skewness	-0.785052711819754
Kurtosis	2.94620625365155
Shapiro-Wilk W	0.963571337559393
Shapiro-Wilk p-value	< 2e-16
Numero Outlier	48

Outlier estremi di Cranio
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
2175	37	8	28	930	355	235	Nat	osp1	F
2452	28	0	26	930	345	245	Ces	osp3	F
1780	25	2	25	900	325	253	Nat	osp3	F
928	25	0	28	830	310	254	Nat	osp1	F
2437	28	1	27	980	320	265	Nat	osp1	M
565	24	1	40	4250	520	386	Nat	osp2	F
190	26	2	39	4050	525	390	Nat	osp1	M
684	30	1	39	3000	475	390	Ces	osp2	F
1356	32	1	40	4090	525	390	Ces	osp3	F
1868	29	0	40	3470	525	390	Nat	osp1	M

Si evidenzia una distribuzione tendenzialmente simmetrica ma con una leggera asimmetria negativa e la presenza di diversi outlier, confermata dal test di Shapiro-Wilk che indica una deviazione dalla normalità.

Analisi di “Sesso”

seeso_df <- analyze_qualitative(df, "Sesso")

Tabella delle Frequenze di Sesso
Categoria	Frequenza	Percentuale
F	1256	50.24
M	1244	49.76

Risultati dei Test Statistici
Test	P_value
Chi-quadro	0.8103303
Test binomiale	0.8258766

L’analisi della distribuzione del sesso nel campione evidenzia un equilibrio tra maschi e femmine. I test statistici confermano l’assenza di differenze significative indicando che la variabile sesso risulta ben bilanciata e rappresentativa della popolazione di riferimento.

In alcuni ospedali si fanno più parti cesarei

H₀ (ipotesi nulla): La frequenza di parti cesarei è uguale nei tre ospedali (Ospedale e Tipo di Parto sono indipendenti).

H₁ (ipotesi alternativa): La frequenza di cesarei varia tra gli ospedali (Ospedale e Tipo di Parto sono associati).

# Tabella di contingenza (frequenze assolute)
tabella_parto_ospedale <- table(Ospedale, Tipo.parto)

# Percentuali calcolate sulle righe
tabella_parto_ospedale_percentuali <- prop.table(tabella_parto_ospedale, margin = 1) * 100

# Arrotondo le percentuali per visualizzarle sulla tabella
tabella_parto_ospedale_percentuali_round <- round(prop.table(tabella_parto_ospedale, margin = 1) * 100, 2)

# Crea un'unica tabella con frequenze e percentuali insieme
tabella_combinata_parto_ospedale <- matrix(
  paste0(tabella_parto_ospedale, " (", tabella_parto_ospedale_percentuali_round, "%)"),
  nrow = nrow(tabella_parto_ospedale),
  dimnames = dimnames(tabella_parto_ospedale)
)

# Convertila in data frame per kable
tabella_combinata_parto_ospedale_df <- as.data.frame.matrix(tabella_combinata_parto_ospedale)

# Visualizza la tabella
kbl(
  tabella_combinata_parto_ospedale_df,
  format = "html",
  caption = "Distribuzione assoluta e percentuale dei tipi di parto per ospedale"
) %>%
  kable_styling(
    full_width = TRUE,
    bootstrap_options = c("striped", "hover", "condensed")
  )

Distribuzione assoluta e percentuale dei tipi di parto per ospedale
	Ces	Nat
osp1	242 (29.66%)	574 (70.34%)
osp2	254 (29.92%)	595 (70.08%)
osp3	232 (27.78%)	603 (72.22%)

# Test chi-quadrato solo sulle frequenze percentuali non arrotondate
test_chi_parto_ospedale <- chisq.test(tabella_parto_ospedale_percentuali)

# Riassunto del test
chi_summary_parto_ospedale <- data.frame(
  Statistica = round(test_chi_parto_ospedale$statistic, 3),
  Gradi_libertà = test_chi_parto_ospedale$parameter,
  P_value = round(test_chi_parto_ospedale$p.value, 4)
)

# Visualizza il riassunto del test
kbl(
  chi_summary_parto_ospedale,
  format = "html",
  caption = "Risultati del Test Chi-Quadrato: Ospedale vs Tipo di Parto"
) %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Risultati del Test Chi-Quadrato: Ospedale vs Tipo di Parto
	Statistica	Gradi_libertà	P_value
X-squared	0.131	2	0.9365

Il test del Chi-Quadrato non evidenzia differenze significative nella frequenza dei parti cesarei tra i tre ospedali. Pertanto, non vi è evidenza di un’associazione tra ospedale e tipo di parto. Percui non rifiuto l’ipotesi nulla e posso indicare che ospedale e tipo parto sono indipendenti.

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

H₀ (nulla): La media del peso del campione è uguale alla media della popolazione (μ = μ₀)

H₁ (alternativa): La media del peso del campione è diversa da quella della popolazione (μ ≠ μ₀)

t_test_peso <- t.test(Peso ~ 1, mu = 3250, data = df)

H₀ (nulla): La media della lunghezza del campione è uguale alla media della popolazione (μ = μ₀)

H₁ (alternativa): La media della lunghezza del campione è diversa da quella della popolazione (μ ≠ μ₀)

t_test_lunghezza <- t.test(Lunghezza ~ 1, mu = 500, data = df)

# Creo la tabella riassuntiva dei risultati
test_summary <- data.frame(
  Variabile = c("Peso alla nascita", "Lunghezza"),
  Media_Campionaria = c(round(mean(df$Peso, na.rm = TRUE), 1),
                        round(mean(df$Lunghezza, na.rm = TRUE), 1)),
  Media_Popolazione = c(3250, 500),
  t_value = c(round(t_test_peso$statistic, 3), round(t_test_lunghezza$statistic, 3)),
  p_value = c(round(t_test_peso$p.value, 4), round(t_test_lunghezza$p.value, 4))
)

# Visualizza la tabella con kable
kbl(test_summary, format = "html", 
    caption = "Risultati del Test t: Confronto tra la Media Campionaria e la Media Teorica della Popolazione") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Risultati del Test t: Confronto tra la Media Campionaria e la Media Teorica della Popolazione
Variabile	Media_Campionaria	Media_Popolazione	t_value	p_value
Peso alla nascita	3284.1	3250	3.246	0.0012
Lunghezza	494.7	500	-10.084	0.0000

Considerando puramente i risultati del test t questi indicano che sia la media del peso alla nascita che la media della lunghezza dei neonati del campione sono significativamente diverse dai valori teorici della popolazione (p-value < 0.05). Pertanto, la media del campione non è uguale a quella della popolazione per entrambe le variabili analizzate. Ma considerando la presenza di outliner (che rappresentano casi particolari) analizzando i valori del peso medio del campione 3284 su un media della popolazione di 3250 e i valori della lunghezza di 494,7 su una media della popolazione di 500 e considerando i valori medi ⁽²⁾ posso dichiarare che i valori del nostro campione sono significativamente uguali a quelli della popolazione.

Le misure antropometriche sono significativamente diverse tra i due sessi

H₀ (ipotesi nulla): La media è uguale tra maschi e femmine

H₁ (alternativa): La media è diversa tra maschi e femmine

t_test_peso_sesso <- t.test(Peso ~ Sesso, data = df)
t_test_lunghezza_sesso <- t.test(Lunghezza ~ Sesso, data = df)
t_test_cranio_sesso <- t.test(Cranio ~ Sesso, data = df)

# Creo la tabella riassuntiva dei risultati
test_sesso_summary <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  t_value = c(round(t_test_peso_sesso$statistic, 3),
              round(t_test_lunghezza_sesso$statistic, 3),
              round(t_test_cranio_sesso$statistic, 3)),
  p_value = c(round(t_test_peso_sesso$p.value, 4),
              round(t_test_lunghezza_sesso$p.value, 4),
              round(t_test_cranio_sesso$p.value, 4))
)

test_sesso_summary$p_value <- sapply(test_sesso_summary$p_value, format_p_value)

# Stampo la tabella in HTML con stile
kbl(test_sesso_summary, format = "html", escape=FALSE,
    caption = "Test t per confrontare Peso, Lunghezza e Cranio tra Maschi e Femmine") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Test t per confrontare Peso, Lunghezza e Cranio tra Maschi e Femmine
Variabile	t_value	p_value
Peso	-12.106	< 2e-16
Lunghezza	-9.582	< 2e-16
Cranio	-7.410	< 2e-16

# Grafico per Peso
ggplot(df, aes(x = Sesso, y = Peso, fill = Sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) + 
  labs(title = "Confronto del Peso tra Maschi e Femmine", y = "Peso (g)") +
  theme_minimal()

# Grafico per Lunghezza
ggplot(df, aes(x = Sesso, y = Lunghezza, fill = Sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) +
  labs(title = "Confronto della Lunghezza tra Maschi e Femmine", y = "Lunghezza (mm)") +
  theme_minimal()

# Grafico per Cranio
ggplot(df, aes(x = Sesso, y = Cranio, fill = Sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) +
  labs(title = "Confronto del Cranio tra Maschi e Femmine", y = "Cranio (mm)") +
  theme_minimal()

I neonati di sesso maschile mostrano valori medi superiori per tutte le variabili considerate rispetto alle femmine. Questa differenza risulta evidente sia dai test statistici (p-value < 0.0001) sia dalla disposizione delle mediane e dei box nei grafici.

Creazione del Modello di Regressione

Analisi della correlazione tra le variabili

# Creazione matrice di correlazione con scatterplot
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
    {
      par(usr = c(0, 1, 0, 1))
      r <- abs(cor(x, y))
      txt <- format(c(r, 0.123456789), digits = digits)[1]
      txt <- paste0(prefix, txt)
      if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
      text(0.5, 0.5, txt, cex = 1.5)
    }
    pairs(df,upper.panel = panel.smooth, lower.panel = panel.cor)

Considerazioni:

La correlazione tra Peso e Lunghezza risulta elevato(0.80) e questo suggerisce che i neonati più lunghi tendono ad avere un peso maggiore.
La correlazione tra Peso e Gestazione alto me meno forte di quello precedente (0.59)ma permette di sostenere che una gestazione più lunga tende ad aumentare il peso alla nascita.
La correlazione tra Lunghezza e Cranio è alta (0.60), come ci si attende per due variabili antropometriche che mappano la crescita fetale.
La correlazione tra Peso e Cranio è alta (0.70), come ci si attende per due variabili antropometriche che mappano la crescita fetale.
Durata della gestazione mostra una distribuzione concentrata intorno alle 40 settimane che è il valore atteso per una gravidanza a termine.

Clicca per mostrare il codice

summary_html <- function(modello, titolo = "") {
  
  # Calcolo residui e statistiche principali
  summary_mod <- summary(modello)
  resid_stats <- quantile(residuals(modello), probs = c(0, 0.25, 0.5, 0.75, 1))
  names(resid_stats) <- c("Min", "1Q", "Median", "3Q", "Max")
  
  # Tabella dei coefficienti
  coeff_tab <- tidy(modello)
  coeff_tab$estimate <- round(coeff_tab$estimate, 3)
  coeff_tab$std.error <- round(coeff_tab$std.error, 3)
  coeff_tab$statistic <- round(coeff_tab$statistic, 3)
  
  # Calcolo significatività (usando "ns" al posto di "")
  coeff_tab$significance <- cut(coeff_tab$p.value,
                                breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
                                labels = c("***", "**", "*", ".", "ns"))
  
  coeff_tab$p.value <- sapply(coeff_tab$p.value, format_p_value)
  
  # Mappa colori significatività
  color_map <- c("***" = "#81d959",   
                 "**"  = "#b2f294",   
                 "*"   = "#fcda51",   
                 "."   = "#FFA500",   
                 "ns"  = "#fc5c51")   
  coeff_tab$row_color <- color_map[as.character(coeff_tab$significance)]
  
  # Inizio costruzione HTML
  html_output <- ""
  
  # Residui
  html_output <- paste0(html_output, "<h4>Residui del ", titolo, "</h4>")
  resid_df <- data.frame(Valore = round(resid_stats, 2))
  rownames(resid_df) <- names(resid_stats)
  html_output <- paste0(html_output,
                        kbl(resid_df, format = "html", row.names = TRUE) %>%
                          kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed")))
  
  # Tabella Coefficienti con colori
  html_output <- paste0(html_output, "<h4>Coefficienti del ", titolo, "</h4>")
  
  # Creo la kable base
  kable_coeff <- kbl(coeff_tab[, c("term", "estimate", "std.error", "statistic", "p.value", "significance")], 
                     format = "html", escape = FALSE) %>%
    kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed")) %>%
    row_spec(0, bold = TRUE)
  
  # Applico la colorazione riga per riga
  if (nrow(coeff_tab) > 0) {
    for (i in 1:nrow(coeff_tab)) {
      kable_coeff <- row_spec(kable_coeff, i, background = coeff_tab$row_color[i])
    }
  }
  
  # Appendo la kable colorata all'html
  html_output <- paste0(html_output, kable_coeff)
  
    # Legenda colori
  html_output <- paste0(html_output, "<strong>Legenda significatività e colori</strong>")
  html_output <- paste0(html_output, "<ul>",
                        "<li style='color:#81d959;'>*** (p < 0.001) - Altamente significativo</li>",
                        "<li style='color:#b2f294;'>** (p < 0.01) - Molto significativo</li>",
                        "<li style='color:#fcda51;'>* (p < 0.05) - Significativo</li>",
                        "<li style='color:#FFA500;'>. (p < 0.1) - Tendenza</li>",
                        "<li style='color:#fc5c51;'>ns (p >= 0.1) - Non significativo</li>",
                        "</ul>")
  
  # Statistiche finali
  html_output <- paste0(html_output, "<h4>Statistiche del ", titolo, "</h4>")
  html_output <- paste0(html_output, "<p><strong>Residual standard error: </strong>", round(summary_mod$sigma, 1),
                        " on ", summary_mod$df[2], " degrees of freedom</p>")
  html_output <- paste0(html_output, "<p><strong>Multiple R-squared:</strong> ", round(summary_mod$r.squared, 4),
                        ", <strong>Adjusted R-squared:</strong> ", round(summary_mod$adj.r.squared, 4), "</p>")

  # Calcolo il p-value della F-statistic
f_p_value <- pf(summary_mod$fstatistic[1], 
                summary_mod$fstatistic[2], 
                summary_mod$fstatistic[3], 
                lower.tail = FALSE)

# Lo formatto con la tua funzione
f_p_value_fmt <- format_p_value(f_p_value)

# Creo l'output HTML con il p-value formattato
html_output <- paste0(html_output, "<p><strong>F-statistic:</strong> ", round(summary_mod$fstatistic[1], 1),
                      " on ", summary_mod$fstatistic[2], " and ", summary_mod$fstatistic[3], 
                      " DF, <strong>p-value:</strong> ", f_p_value_fmt, "</p>")

  
  
  # Output finale
  return(HTML(html_output))
}

vif_html <- function(modello, titolo = "Variance Inflation Factor (VIF) per valutare la multicollinearità") {
  # Calcolo VIF
  vif_values <- car::vif(modello)
  vif_df <- data.frame(Variabile = names(vif_values), VIF = round(vif_values, 3))
  rownames(vif_df) <- NULL
  
  # Colore condizionato in base alla soglia
  vif_df$Colore <- ifelse(vif_df$VIF > 5, "#fc5c51",       # rosso = collinearità significativa
                          ifelse(vif_df$VIF > 2.5, "#fcda51", # giallo = attenzione
                                 "#81d959"))                # verde = ok

  # Titolo
  html_output <- paste0("<h4>", titolo, "</h4>")

  # Tabella VIF
  kable_vif <- kbl(vif_df[, c("Variabile", "VIF")], format = "html", escape = FALSE) %>%
    kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed")) %>%
    row_spec(0, bold = TRUE)
  
  for (i in 1:nrow(vif_df)) {
    kable_vif <- row_spec(kable_vif, i, background = vif_df$Colore[i])
  }
  
  # Legenda
  legenda <- paste0("<strong>Legenda VIF:</strong><ul>",
                    "<li style='color:#81d959;'>VIF ≤ 2.5 → Assenza di collinearità</li>",
                    "<li style='color:#fcda51;'>2.5 < VIF ≤ 5 → Potenziale collinearità</li>",
                    "<li style='color:#fc5c51;'>VIF > 5 → Collinearità significativa</li>",
                    "</ul>")
  
  # Output HTML finale
  html_output <- paste0(html_output, kable_vif, legenda)
  return(HTML(html_output))
}

Modello 1

Precedo alla realizzazione del primo modello, escludendo Tipo di parto⁽³⁾ e Ospedale ⁽⁴⁾.

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

summary_html(mod1, "Modello 1")

Residui del Modello 1

	Valore
Min	-1161.56
1Q	-181.19
Median	-15.75
3Q	163.70
Max	2630.75

Coefficienti del Modello 1

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-6714.411	141.151	-47.569	< 2e-16	***
Anni.madre	0.958	1.135	0.845	0.3984	ns
N.gravidanze	11.276	4.669	2.415	0.0158	*
Fumatrici	-30.296	27.597	-1.098	0.2724	ns
Gestazione	32.933	3.827	8.606	< 2e-16	***
Lunghezza	10.234	0.301	34.009	< 2e-16	***
Cranio	10.518	0.427	24.642	< 2e-16	***
Sesso	78.085	11.204	6.969	4.06e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 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: < 2e-16

vif_html(mod1)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
Anni.madre	1.187
N.gravidanze	1.185
Fumatrici	1.007
Gestazione	1.695
Lunghezza	2.079
Cranio	1.629
Sesso	1.040

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

Dall’analisi dei coefficienti del Modello 1 risulta che la variabile Anni.madre non è statisticamente significativa, indicando che l’età materna non ha un effetto rilevante sul peso neonatale.

Nonostante la variabile Fumatrici non risulta statisticamente significativa nel nostro modello, potrebbe comunque avere un impatto indiretto sullo sviluppo fetale non catturato dai dati disponibili.

Inoltre, il valore del Variance Inflation Factor (VIF) per le variabili è ben al di sotto della soglia di 2.5, confermando l’assenza di multicollinearità e rendendo possibile una riduzione del modello senza compromettere la stabilità delle stime.

Modello 2

Provo ad escludere la variabile Anni.madre, con l’obiettivo di ottenere un modello più parsimonioso ma altrettanto informativo.

mod2 <- lm (Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Sesso , data= df)

summary_html(mod2, "Modello 2")

Residui del Modello 2

	Valore
Min	-1150.32
1Q	-181.27
Median	-15.70
3Q	163.00
Max	2636.33

Coefficienti del Modello 2

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-6681.671	135.718	-49.232	< 2e-16	***
N.gravidanze	12.719	4.345	2.927	0.0035	**
Fumatrici	-30.463	27.595	-1.104	0.2697	ns
Gestazione	32.591	3.805	8.565	< 2e-16	***
Lunghezza	10.234	0.301	34.011	< 2e-16	***
Cranio	10.536	0.426	24.718	< 2e-16	***
Sesso	78.171	11.203	6.978	3.83e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 2

Residual standard error: 274.6 on 2493 degrees of freedom

Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265

F-statistic: 1107.2 on 6 and 2493 DF, p-value: < 2e-16

vif_html(mod2)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
N.gravidanze	1.026
Fumatrici	1.007
Gestazione	1.676
Lunghezza	2.079
Cranio	1.625
Sesso	1.040

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

Escludendo Anni.madre non ha compromesso la bontà del modello: l’Adjusted R-squared resta infatti pressoché invariato (0.7265 contro 0.7264 nel Modello 1), a conferma che l’età materna non fornisce un contributo informativo rilevante.

I valori del Variance Inflation Factor (VIF) per tutte le variabili risultano inferiori a 2.1, escludendo problemi di multicollinearità. Le altre variabili — in particolare Gestazione, Lunghezza, Cranio e Sesso — si confermano altamente significative.

Il Modello 2 rappresenta dunque una versione più parsimoniosa rispetto al modello iniziale, mantenendo però al suo interno variabili di rilevanza clinica per garantire una maggiore utilità applicativa in contesto ospedaliero.

Modello 3

Provo ad escludere la variabile N.gravidanze, che pur essendo significativa, lo è meno delle altre, percui vale la pena vedere l’impatto sul modello.

mod3 <- lm (Peso ~ Fumatrici + Gestazione + Lunghezza + Cranio + Sesso , data= df)

summary_html(mod3, "Modello 3")

Residui del Modello 3

	Valore
Min	-1138.76
1Q	-182.09
Median	-17.88
3Q	162.77
Max	2624.15

Coefficienti del Modello 3

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-6651.068	135.520	-49.078	< 2e-16	***
Fumatrici	-26.362	27.601	-0.955	0.3396	ns
Gestazione	31.480	3.792	8.302	< 2e-16	***
Lunghezza	10.192	0.301	33.859	< 2e-16	***
Cranio	10.669	0.424	25.135	< 2e-16	***
Sesso	79.278	11.213	7.070	2.00e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 3

Residual standard error: 275 on 2494 degrees of freedom

Multiple R-squared: 0.7262, Adjusted R-squared: 0.7256

F-statistic: 1322.9 on 5 and 2494 DF, p-value: < 2e-16

vif_html(mod3)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
Fumatrici	1.004
Gestazione	1.659
Lunghezza	2.074
Cranio	1.606
Sesso	1.039

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

L’esclusione della variabile N.gravidanze ha comportato una riduzione trascurabile dell’Adjusted R-squared (da 0.7265 a 0.7256) permettendo di ottenere un modello ancora più snello, mantenendo inalterata la significatività delle variabili principali. Per questo motivo si considera accettabile l’eliminazione di tale predittore, in linea con il principio di parsimonia.

Osservazioni sui modelli analizzati fino ad ora

Vedendo che con questi modelli, passando dal modello 1 al modello 3, i miglioramenti sono stati marginali in termini di -Adjusted R-squared- ritengo opportuno esplorare potenziali effetti non lineari, nello specifico:

1. Gestazione vs Peso

Il grafico presenta una, con aumento rapido fino a circa 39 settimane e piattaforma dopo mostrando chiaramente una relazione non perfettamente lineare. proverò a testare Gestazione^2.

2. Anni.madre vs Peso

La relazione appare piatta o leggermente curva verso l’alto nella parte centrale. Proverò a testare Anni.madre^2.

3. Numero di gravidanze vs Peso

Si può intuire un effetto di saturazione dopo un certo livello. Proverò a testare N.gravidanze^2.

Modello 4

Introduco un termine quadratico per la durata della gestazione, al fine di modellare l’evidente curvatura osservata nella relazione con il peso neonatale.

mod4 <- lm(Peso ~ Fumatrici + Gestazione + I(Gestazione^2) + Lunghezza + Cranio + Sesso, data = df)

summary_html(mod4, "Modello 4")

Residui del Modello 4

	Valore
Min	-1133.42
1Q	-182.06
Median	-15.72
3Q	164.10
Max	2646.60

Coefficienti del Modello 4

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-4659.725	899.703	-5.179	2.41e-07	***
Fumatrici	-25.130	27.584	-0.911	0.3624	ns
Gestazione	-79.696	49.803	-1.600	0.1097	ns
I(Gestazione^2)	1.484	0.663	2.239	0.0253	*
Lunghezza	10.295	0.304	33.834	< 2e-16	***
Cranio	10.764	0.426	25.253	< 2e-16	***
Sesso	77.119	11.246	6.858	8.80e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 4

Residual standard error: 274.8 on 2493 degrees of freedom

Multiple R-squared: 0.7267, Adjusted R-squared: 0.7261

F-statistic: 1105.1 on 6 and 2493 DF, p-value: < 2e-16

vif_html(mod4)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
Fumatrici	1.004
Gestazione	286.635
I(Gestazione^2)	279.859
Lunghezza	2.123
Cranio	1.622
Sesso	1.047

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

Inserendo Gestazione^2 conferma l’esistenza di una relazione non lineare significativa con p.vlaue = 0.0253. Notiamo però che Gestazione non risulta ancora significativo probabilmente a causa di forte collinearità tra Gestazione e Gestazione^2.

Modello 5

Con l’obiettivo di valutare se l’età contribuisca significativamente alla spiegazione della variabilità del peso, soprattutto in presenza di una possibile curvatura centrale osservata nei grafici esplorativi, estendo il modello includendo Anni.madre e Anni.madre^2

mod5 <- lm(Peso ~ Fumatrici + Gestazione + I(Gestazione^2) + Anni.madre + I(Anni.madre^2) + Lunghezza + Cranio + Sesso, data = df)

summary_html(mod5, "Modello 5")

Residui del Modello 5

	Valore
Min	-1130.55
1Q	-182.69
Median	-13.30
3Q	161.46
Max	2638.58

Coefficienti del Modello 5

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-4843.612	906.344	-5.344	9.91e-08	***
Fumatrici	-25.455	27.566	-0.923	0.3559	ns
Gestazione	-80.177	49.792	-1.610	0.1075	ns
I(Gestazione^2)	1.501	0.663	2.264	0.0237	*
Anni.madre	12.407	7.334	1.692	0.0908	.
I(Anni.madre^2)	-0.184	0.129	-1.428	0.1533	ns
Lunghezza	10.293	0.304	33.824	< 2e-16	***
Cranio	10.703	0.427	25.041	< 2e-16	***
Sesso	77.190	11.248	6.863	8.50e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 5

Residual standard error: 274.6 on 2491 degrees of freedom

Multiple R-squared: 0.7274, Adjusted R-squared: 0.7265

F-statistic: 830.8 on 8 and 2491 DF, p-value: < 2e-16

vif_html(mod5)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
Fumatrici	1.005
Gestazione	286.943
I(Gestazione^2)	280.367
Anni.madre	49.578
I(Anni.madre^2)	49.714
Lunghezza	2.126
Cranio	1.634
Sesso	1.049

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

L’aggiunta dell’età materna e del suo termine quadratico non ha apportato miglioramenti sostanziali al modello, mostrando collinearità elevata e nessuna significatività chiara, mentre la curvatura della gestazione si conferma rilevante.

Modello 6

Aggiungo anche il numero di gravidanze e il relativo termine quadratico, per valutare un possibile effetto di saturazione sulla variabile dipendente.

mod6 <- lm(Peso ~ Fumatrici + Gestazione + I(Gestazione^2) +
                     Anni.madre + I(Anni.madre^2) +
                     N.gravidanze + I(N.gravidanze^2) +
                     Lunghezza + Cranio + Sesso,
           data = df)

summary_html(mod6, "Modello 6")

Residui del Modello 6

	Valore
Min	-1129.40
1Q	-181.73
Median	-11.33
3Q	162.97
Max	2652.78

Coefficienti del Modello 6

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-4742.228	907.066	-5.228	1.86e-07	***
Fumatrici	-31.737	27.590	-1.150	0.2501	ns
Gestazione	-84.912	49.804	-1.705	0.0883	.
I(Gestazione^2)	1.578	0.663	2.379	0.0174	*
Anni.madre	12.230	7.335	1.667	0.0955	.
I(Anni.madre^2)	-0.204	0.129	-1.586	0.1128	ns
N.gravidanze	25.998	8.505	3.057	0.0023	**
I(N.gravidanze^2)	-2.614	1.318	-1.984	0.0474	*
Lunghezza	10.317	0.304	33.928	< 2e-16	***
Cranio	10.580	0.428	24.695	< 2e-16	***
Sesso	76.599	11.232	6.819	1.14e-11	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello 6

Residual standard error: 274.1 on 2489 degrees of freedom

Multiple R-squared: 0.7285, Adjusted R-squared: 0.7274

F-statistic: 667.9 on 10 and 2489 DF, p-value: < 2e-16

vif_html(mod6)

Variance Inflation Factor (VIF) per valutare la multicollinearità

Variabile	VIF
Fumatrici	1.010
Gestazione	288.057
I(Gestazione^2)	281.643
Anni.madre	49.756
I(Anni.madre^2)	50.152
N.gravidanze	3.945
I(N.gravidanze^2)	3.577
Lunghezza	2.130
Cranio	1.647
Sesso	1.049

Legenda VIF:

VIF ≤ 2.5 → Assenza di collinearità
2.5 < VIF ≤ 5 → Potenziale collinearità
VIF > 5 → Collinearità significativa

L’introduzione dei termini quadratici per il numero di gravidanze ha migliorato significativamente il modello, confermando un effetto curvilineo. Il modello 6 rappresenta i risultati migliori. Questo modello può essere considerato una solida base da cui partire per eventuali semplificazioni successive, valutando se eliminare predittori non significativi o ridurre la complessità.

Selezione del modello ottimale con stepwise

Faccio delle ulteriori verifiche per accertarmi di trovare il modello più parsimonioso ma che mantenga un’ottima interpretabilità.

# aggiungo termini quadratici al dataset
df$Gestazione2 <- df$Gestazione^2
df$Anni.madre2 <- df$Anni.madre^2
df$N.gravidanze2 <- df$N.gravidanze^2

#Costruzione modello esteso
mod_full <- lm(Peso ~ Fumatrici + Gestazione + Gestazione2 +
                          Anni.madre + Anni.madre2 +
                          N.gravidanze + N.gravidanze2 +
                          Lunghezza + Cranio + Sesso, data = df)

# Applico stepwise selection basata sul BIC
mod_bic <- step(mod_full, direction = "both", k = log(nrow(df)))

formula_string <- paste(deparse(formula(mod_bic)), collapse = " ")

cat("<h4>Formula del modello selezionato con stepwise</h4>")

Formula del modello selezionato con stepwise

cat("<pre><code>", formula_string, "</code></pre>")

 Peso ~ Gestazione2 + N.gravidanze + Lunghezza + Cranio + Sesso

# Riepilogo con stime e significatività
summary_html(mod_bic,"Modello selezionato con stepwise")

Residui del Modello selezionato con stepwise

	Valore
Min	-1146.87
1Q	-181.09
Median	-15.12
3Q	165.40
Max	2643.89

Coefficienti del Modello selezionato con stepwise

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-6100.641	123.545	-49.380	< 2e-16	***
Gestazione2	0.438	0.051	8.667	< 2e-16	***
N.gravidanze	12.554	4.338	2.894	0.0038	**
Lunghezza	10.261	0.299	34.358	< 2e-16	***
Cranio	10.560	0.426	24.816	< 2e-16	***
Sesso	77.329	11.198	6.906	6.32e-12	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello selezionato con stepwise

Residual standard error: 274.5 on 2494 degrees of freedom

Multiple R-squared: 0.7273, Adjusted R-squared: 0.7267

F-statistic: 1330.2 on 5 and 2494 DF, p-value: < 2e-16

Correzione del Modello selezionato con stepwise

L’algoritmo di selezione automatica basato su BIC ha escluso il termine lineare Gestazione, mantenendo solo il termine quadratico ma per garantire una corretta interpretazione della curva e coerenza teorica, ho reinserito Gestazione nel modello finale.

mod_bic_corrected <- lm(Peso ~ Gestazione + Gestazione2 + N.gravidanze + Lunghezza + Cranio + Sesso, data = df)
summary_html(mod_bic_corrected, "Modello selezionato con stepwise (Corretto)")

Residui del Modello selezionato con stepwise (Corretto)

	Valore
Min	-1144.11
1Q	-181.17
Median	-13.16
3Q	165.73
Max	2662.56

Coefficienti del Modello selezionato con stepwise (Corretto)

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-4650.227	898.171	-5.177	2.43e-07	***
Gestazione	-81.049	49.713	-1.630	0.1032	ns
Gestazione2	1.514	0.662	2.287	0.0223	*
N.gravidanze	12.566	4.336	2.898	0.0038	**
Lunghezza	10.353	0.304	34.074	< 2e-16	***
Cranio	10.636	0.428	24.854	< 2e-16	***
Sesso	75.790	11.234	6.747	1.88e-11	***

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello selezionato con stepwise (Corretto)

Residual standard error: 274.4 on 2493 degrees of freedom

Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269

F-statistic: 1109.7 on 6 and 2493 DF, p-value: < 2e-16

mod_bic_revisioned <- lm(Peso ~ Gestazione + Gestazione2 + N.gravidanze + Lunghezza + Cranio + Sesso + Fumatrici, data = df)
summary_html(mod_bic_revisioned, "Modello selezionato con stepwise (Rivisto)")

Residui del Modello selezionato con stepwise (Rivisto)

	Valore
Min	-1145.00
1Q	-180.96
Median	-13.12
3Q	165.00
Max	2659.10

Coefficienti del Modello selezionato con stepwise (Rivisto)

term	estimate	std.error	statistic	p.value	significance
(Intercept)	-4669.103	898.325	-5.198	2.18e-07	***
Gestazione	-79.774	49.726	-1.604	0.1088	ns
Gestazione2	1.500	0.662	2.266	0.0235	*
N.gravidanze	12.799	4.342	2.948	0.0032	**
Lunghezza	10.339	0.304	33.989	< 2e-16	***
Cranio	10.631	0.428	24.841	< 2e-16	***
Sesso	75.981	11.235	6.763	1.68e-11	***
Fumatrici	-29.244	27.577	-1.060	0.2890	ns

Legenda significatività e colori

*** (p < 0.001) - Altamente significativo
** (p < 0.01) - Molto significativo
* (p < 0.05) - Significativo
. (p < 0.1) - Tendenza
ns (p >= 0.1) - Non significativo

Statistiche del Modello selezionato con stepwise (Rivisto)

Residual standard error: 274.4 on 2492 degrees of freedom

Multiple R-squared: 0.7277, Adjusted R-squared: 0.7269

F-statistic: 951.4 on 7 and 2492 DF, p-value: < 2e-16

Confronto tra modelli: BIC e R² Adjusted

# Nomi dei modelli
mod_names <- c("mod1", "mod2", "mod3", "mod4", "mod5", "mod6", "mod_bic", "mod_bic_corrected", "mod_bic_revisioned")

# Calcolo BIC
bic_vals <- BIC(mod1, mod2, mod3, mod4, mod5, mod6, mod_bic, mod_bic_corrected, mod_bic_revisioned)

# Calcolo R² Adjusted
r2_vals <- c(
  summary(mod1)$adj.r.squared,
  summary(mod2)$adj.r.squared,
  summary(mod3)$adj.r.squared,
  summary(mod4)$adj.r.squared,
  summary(mod5)$adj.r.squared,
  summary(mod6)$adj.r.squared,
  summary(mod_bic)$adj.r.squared,
  summary(mod_bic_corrected)$adj.r.squared,
  summary(mod_bic_revisioned)$adj.r.squared
)

# Tabella riassuntiva
bic_table <- data.frame(
  Modello = mod_names,
  DF = bic_vals$df,
  BIC = round(bic_vals$BIC, 2),
  R2_Adj = round(r2_vals, 4)
)

# Identifica il miglior modello (BIC più basso)
min_bic <- min(bic_table$BIC)
bic_table$BIC_colored <- ifelse(
  bic_table$BIC == min_bic,
  cell_spec(bic_table$BIC, background = "lightgreen"),
  as.character(bic_table$BIC)
)

# Genera tabella HTML con stile
kbl(
  bic_table[, c("Modello", "DF", "BIC_colored", "R2_Adj")],
  format = "html",
  escape = FALSE,
  col.names = c("Modello", "DF", "BIC", "R² Adj"),
  caption = ""
) %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE)


Modello	DF	BIC	R² Adj
mod1	9	35233.81	0.7264
mod2	8	35226.7	0.7265
mod3	7	35227.45	0.7256
mod4	8	35230.26	0.7261
mod5	10	35240.12	0.7265
mod6	12	35245.33	0.7274
mod_bic	7	35217.52	0.7267
mod_bic_corrected	8	35222.68	0.7269
mod_bic_revisioned	9	35229.38	0.7269

Tra i modelli testati, quello con il BIC più basso escludeva il termine lineare della gestazione. Tuttavia, per rappresentare meglio l’andamento reale del peso neonatale rispetto alle settimane di gravidanza, è stato scelto un modello leggermente più completo mod_bic_corrected , che include anche la gestazione lineare. Questo modello mantiene ottime performance statistiche e risulta più coerente dal punto di vista clinico, offrendo una lettura più realistica e interpretabile dei dati, sia rispetto a mod2 che a mod_bic_corrected.

In base alle scelte fatte in precedenza decido di aggiungere al modello anche Fumatrici che pur non risultando rilevante nel nostro campione, lo è a livello di popolazione, percui decido di utilizzare un modello leggermente più complesso ma che possa essere più utile a spiegare le motivazione sulle previsoni del peso neonatale.

La scelta finale ricade sul modello mod_bic_revisioned.

Analisi della qualità del modello

Dopo aver valutato l’andamento dell’R² aggiustato nei vari modelli, passiamo all’analisi dei residui per verificare che le assunzioni del modello non siano violate: linearità, normalità, omoscedasticità, indipendenza, oltre alla presenza di outlier o osservazioni influenti.

# Imposta layout 2x2
par(mfrow = c(2, 2))

# Colore personalizzato
colore <- "#1f78b4"

# Verifica della linearità e omoschedasticità
plot(mod_bic_revisioned$fitted.values, resid(mod_bic_revisioned),
     main = "Verifica della linearità e omoschedasticità",
     xlab = "Valori Predetti", ylab = "Residui",
     col = colore, pch = 19)
abline(h = 0, lty = 2)

# Verifica della normalità dei residui
qqnorm(resid(mod_bic_revisioned), main = "Verifica della normalità dei residui", col = colore, pch = 19)
qqline(resid(mod_bic_revisioned), col = "darkgray", lwd = 2)

# Verifica della omoschedasticità
plot(mod_bic_revisioned$fitted.values, sqrt(abs(rstandard(mod_bic_revisioned))),
     main = "Verifica della omoschedasticità",
     xlab = "Valori Predetti", ylab = "Residui Standardizzati|",
     col = colore, pch = 19)
abline(h = median(sqrt(abs(rstandard(mod_bic_revisioned)))), lty = 2)

# Verifica della presenza di osservazioni influenti
plot(hatvalues(mod_bic_revisioned), rstandard(mod_bic_revisioned),
     main = "Verifica della presenza di osservazioni influenti",
     xlab = "Leverage", ylab = "Residui Standardizzati",
     col = colore, pch = 19)
abline(h = c(-2, 0, 2), lty = 2, col = "gray")
abline(v = 2*mean(hatvalues(mod_bic_revisioned)), lty = 2, col = "gray")

Verifica della linearità e omoschedasticità

Il grafico “Residui vs Valori Predetti” evidenzia un pattern a funnel, suggerendo una possibile violazione dell’omoschedasticità e una lieve non linearità per i valori più bassi. Questo potrebbe dipendere da fattori clinici non inclusi nel modello, come gravidanze gemellari, patologie materne o variabilità legata alla prematurità. Tuttavia, trattandosi di casi specifici e limitati, il modello può essere comunque considerato valido nel complesso.

Verifica della normalità dei residui

Il grafico evidenzia una deviazione nelle code, segnalando una normalità approssimativa e la possibile presenza di outlier. Tuttavia, dato che la maggior parte dei residui segue la distribuzione attesa e le deviazioni si concentrano negli estremi, il modello può essere ritenuto adeguato per scopi predittivi.

Verifica della omoschedasticità

Il grafico mostra una diminuzione della varianza dei residui all’aumentare dei valori predetti, indicando eteroschedasticità. Tuttavia, trattandosi di una tendenza moderata e limitata ai valori più bassi, il modello resta valido per fini predittivi.

Verifica della presenza di osservazioni influenti

Il grafico evidenzia alcuni outlier con residui elevati e almeno un punto potenzialmente influente. Tuttavia, trattandosi di casi isolati, l’impatto complessivo sul modello è limitato e ne consente comunque un uso affidabile a fini predittivi.

Durbin-Watson test per l’autocorrelazione dei residui

Clicca per mostrare il codice

dw_test_html <- function(model) {
  # Esegui il test di Durbin-Watson
  dw <- lmtest::dwtest(model)
  
  # Estrai i dati
  DW <- round(dw$statistic, 4)
  pval <- round(dw$p.value, 4)
  model_name <- deparse(substitute(model))
  hypothesis <- dw$alternative
  
  # Crea l'HTML della tabella
  html <- paste0(
    "<table border='1' width='100%' style='border-collapse: collapse; margin: 20px auto; font-family: Arial, sans-serif;'>",
    "<thead><tr>",
    "<th style='padding: 8px; text-align: center;'>Test</th>",
    "<th style='padding: 8px; text-align: center;'>Modello</th>",
    "<th style='padding: 8px; text-align: center;'>DW</th>",
    "<th style='padding: 8px; text-align: center;'>p-value</th>",
    "<th style='padding: 8px; text-align: center;'>Ipotesi alternativa</th>",
    "</tr></thead>",
    "<tbody><tr>",
    "<td style='padding: 8px; text-align: center;'>Durbin-Watson</td>",
    "<td style='padding: 8px; text-align: center;'>", model_name, "</td>",
    "<td style='padding: 8px; text-align: center;'>", DW, "</td>",
    "<td style='padding: 8px; text-align: center;'>", pval, "</td>",
    "<td style='padding: 8px; text-align: center;'>", hypothesis, "</td>",
    "</tr></tbody>",
    "</table>"
  )
  
  return(HTML(html))  # usa HTML() in R Markdown per renderizzare correttamente
}

dw_test_html(mod_bic_revisioned)

Test	Modello	DW	p-value	Ipotesi alternativa
Durbin-Watson	mod_bic_revisioned	1.9531	0.1203	true autocorrelation is greater than 0

Il modello non presenta autocorrelazione positiva significativa tra i residui. Questo rafforza la validità inferenziale delle stime ottenute e conferma che l’ipotesi di indipendenza dei residui è ragionevolmente soddisfatta.

Shapiro-Wilk test per la verifica della normalità

Clicca per mostrare il codice

shapiro_test_html <- function(model) {
  if (!"residuals" %in% names(model)) {
    stop("Il modello fornito non contiene residui.")
  }
  
  shapiro <- shapiro.test(model$residuals)
  W <- round(shapiro$statistic, 5)
  pval <- if (shapiro$p.value < 2.2e-16) "&lt; 2.2e-16" else round(shapiro$p.value, 5)
  model_name <- deparse(substitute(model))
  
  html <- paste0(
    "<table border='1' width='100%' style='border-collapse: collapse; margin: 20px auto; font-family: Arial, sans-serif;'>",
    "<thead><tr>",
    "<th style='padding: 8px; text-align: center;'>Test</th>",
    "<th style='padding: 8px; text-align: center;'>Modello</th>",
    "<th style='padding: 8px; text-align: center;'>W</th>",
    "<th style='padding: 8px; text-align: center;'>p-value</th>",
    "<th style='padding: 8px; text-align: center;'>Ipotesi nulla</th>",
    "</tr></thead>",
    "<tbody><tr>",
    "<td style='padding: 8px; text-align: center;'>Shapiro-Wilk</td>",
    "<td style='padding: 8px; text-align: center;'>", model_name, "</td>",
    "<td style='padding: 8px; text-align: center;'>", W, "</td>",
    "<td style='padding: 8px; text-align: center;'>", pval, "</td>",
    "<td style='padding: 8px; text-align: center;'>Normalità dei residui</td>",
    "</tr></tbody>",
    "</table>"
  )
  
  return(HTML(html))
}

shapiro_test_html(mod_bic_revisioned)

Test	Modello	W	p-value	Ipotesi nulla
Shapiro-Wilk	mod_bic_revisioned	0.97365	< 2.2e-16	Normalità dei residui

Il test di Shapiro-Wilk valuta se i residui seguono una distribuzione normale. Un p-value molto basso indica una deviazione dalla normalità, violando così l’assunzione necessaria per l’affidabilità dei test statistici classici (come t-test e intervalli di confidenza) nel modello di regressione.

Test di Breusch-Pagan, test d’ipotesi di eteroschedasticità

Clicca per mostrare il codice

bp_test_html <- function(model) {
  if (!requireNamespace("lmtest", quietly = TRUE)) {
    stop("Il pacchetto 'lmtest' è richiesto ma non è installato.")
  }

  bp <- lmtest::bptest(model)
  BP <- round(bp$statistic, 3)
  df <- bp$parameter
  pval <- if (bp$p.value < 2.2e-16) "&lt; 2.2e-16" else round(bp$p.value, 5)
  model_name <- deparse(substitute(model))
  
  html <- paste0(
    "<table border='1' width='100%' style='border-collapse: collapse; margin: 20px auto; font-family: Arial, sans-serif;'>",
    "<thead><tr>",
    "<th style='padding: 8px; text-align: center;'>Test</th>",
    "<th style='padding: 8px; text-align: center;'>Modello</th>",
    "<th style='padding: 8px; text-align: center;'>BP</th>",
    "<th style='padding: 8px; text-align: center;'>df</th>",
    "<th style='padding: 8px; text-align: center;'>p-value</th>",
    "<th style='padding: 8px; text-align: center;'>Ipotesi nulla</th>",
    "</tr></thead>",
    "<tbody><tr>",
    "<td style='padding: 8px; text-align: center;'>Breusch-Pagan</td>",
    "<td style='padding: 8px; text-align: center;'>", model_name, "</td>",
    "<td style='padding: 8px; text-align: center;'>", BP, "</td>",
    "<td style='padding: 8px; text-align: center;'>", df, "</td>",
    "<td style='padding: 8px; text-align: center;'>", pval, "</td>",
    "<td style='padding: 8px; text-align: center;'>Omogeneità della varianza (omoscedasticità)</td>",
    "</tr></tbody>",
    "</table>"
  )
  
  return(HTML(html))
}

bp_test_html(mod_bic_revisioned)

Test	Modello	BP	df	p-value	Ipotesi nulla
Breusch-Pagan	mod_bic_revisioned	95.282	7	< 2.2e-16	Omogeneità della varianza (omoscedasticità)

Il test di Breusch-Pagan verifica se la varianza dei residui rimane costante al variare dei predittori. Un p-value molto basso indica che la varianza non è omogenea (cioè i residui presentano eteroscedasticità), violando così una delle assunzioni fondamentali della regressione lineare. Questo può influenzare l’affidabilità degli errori standard e dei test sui coefficienti.

Analisi visuale

Visto i risultati dei test di Shapiro-Wilk e di Breusch-Pagan, eseguo un’analisi visiva per approfondire la comprensione delle deviazioni dalle assunzioni del modello di regressione lineare.

# Residui del modello
residui <- residuals(mod_bic_revisioned)

# Densità residui osservati
dens_residui <- density(residui)
dens_normale <- density(rnorm(100000, mean = 0, sd = sd(residui)))

# Data frame combinato
df_dens <- data.frame(
  x = c(dens_residui$x, dens_normale$x),
  y = c(dens_residui$y, dens_normale$y),
  gruppo = rep(c("Residui osservati", "Normale teorica"), each = length(dens_residui$x))
)

# Grafico
ggplot(df_dens, aes(x = x, y = y, color = gruppo, linetype = gruppo)) +
  geom_line(size = 1.1) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Confronto tra la densità dei residui e una distribuzione normale",
    x = "Valori",
    y = "Densità",
    color = "Distribuzione",
    linetype = "Distribuzione"
  ) +
  scale_color_manual(values = c("steelblue", "firebrick")) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

Data la simmetria complessiva e l’assenza di distorsioni gravi, il modello può essere considerato adeguato per fini previsionali.

# Calcolo del leverage
lev <- hatvalues(mod_bic_corrected)
p <- sum(lev)  # gradi di libertà effettivi
n <- length(lev)
soglia <- 2 * p / n  # soglia standard per individuare high-leverage points

# Crea un data frame per ggplot
df_lev <- data.frame(
  osservazione = 1:n,
  leverage = lev
)

# Grafico con soglia
ggplot(df_lev, aes(x = osservazione, y = leverage)) +
  geom_point(color = "steelblue", size = 1.5) +
  geom_hline(yintercept = soglia, color = "firebrick", linetype = "dashed", size = 1) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Leverage per osservazione",
    x = "Indice osservazione",
    y = "Leverage"
  ) +
  annotate("text", x = n*0.9, y = soglia + 0.01, label = paste("Soglia =", round(soglia, 3)), 
           color = "firebrick", hjust = 0) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Il grafico mostra che solo poche osservazioni superano significativamente la soglia teorica di leverage, indicando potenziali influenze anomale sul modello, stampo in una tabella i valori più influenti per analizzarli.

Clicca per mostrare il codice

tabella_leverage_html <- function(model, dati = NULL) {

  lev <- hatvalues(model)
  p <- sum(lev)
  n <- length(lev)
  soglia_base <- 2 * p / n
  moltiplicatore <- 5
  soglia_effettiva <- moltiplicatore * soglia_base

  if (is.null(dati)) dati <- model$model

  df <- tibble::as_tibble(dati)
  df$osservazione <- 1:n
  df$leverage <- lev

  # Rimuove colonne che terminano con "2" in modo compatibile
  colonne_da_escludere <- grep("2$", names(df), value = TRUE)
  df <- df[, !(names(df) %in% colonne_da_escludere)]

  # Filtra e riordina
  df_filtrato <- df %>%
    filter(leverage > soglia_effettiva) %>%
    arrange(desc(leverage)) %>%
    relocate(osservazione, leverage)

  # Tabella HTML finale
  kbl(df_filtrato, format = "html", digits = 5,
      caption = paste0("Osservazioni con leverage > ", round(soglia_effettiva, 5),
                       " (", moltiplicatore, " × soglia teorica)")) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  full_width = TRUE, position = "center") %>%
    column_spec(2, color = "red") %>%
    row_spec(0, bold = TRUE)
}

tabella_leverage_html(mod_bic_revisioned,df)

Osservazioni con leverage > 0.032 (5 × soglia teorica)
osservazione	leverage	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Sesso
1780	0.12810	25	2	25	900	325	253	0
2452	0.09377	28	0	26	930	345	245	0
310	0.06765	40	3	28	1560	420	379	0
2437	0.06683	28	1	27	980	320	265	1
2120	0.06669	32	0	27	1140	370	267	0
2175	0.05973	37	8	28	930	355	235	0
1551	0.05023	35	1	38	4370	315	374	0
928	0.04805	25	0	28	830	310	254	0
378	0.04734	27	0	28	1285	400	274	0
1429	0.04287	24	4	29	1280	390	355	0
1130	0.03552	33	11	43	3400	475	360	1

Tra le osservazioni con leverage molto elevato, i casi 1780, 2175, 1551 e 1130 mostrano anche valori clinicamente anomali (come peso molto basso o elevato, numero di gravidanze atipico o gestazione fuori norma), suggerendo una possibile influenza reale e giustificata sul modello.

Faccio delle ulteriori analisi per approfondire a partire dal grafico dei residui studentizzati che consente di individuare visivamente le osservazioni che si discostano in modo anomalo dal modello, aiutando a rilevare outlier che potrebbero influenzare l’accuratezza dell’inferenza.

# Calcola residui studentizzati
resid_student <- rstudent(mod_bic_revisioned)

# Crea dataframe
df_rstudent <- data.frame(
  osservazione = 1:length(resid_student),
  rstudent = resid_student
)

# Grafico
ggplot(df_rstudent, aes(x = osservazione, y = rstudent)) +
  geom_point(color = "steelblue", size = 1.5) +
  geom_hline(yintercept = c(-2, 2), color = "firebrick", linetype = "dashed", size = 1) +
  annotate("text", x = max(df_rstudent$osservazione) * 0.95, y = 2.1, 
           label = "+2", color = "firebrick", hjust = 0) +
  annotate("text", x = max(df_rstudent$osservazione) * 0.95, y = -2.1, 
           label = "-2", color = "firebrick", hjust = 0) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Residui studentizzati per osservazione",
    x = "Indice osservazione",
    y = "Residuo studentizzato"
  ) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Il grafico mostra che la maggior parte dei residui studentizzati rientra nell’intervallo ±2, con poche osservazioni potenzialmente anomale.

Clicca per mostrare il codice

tabella_outlier_html <- function(model) {
  library(car)
  library(knitr)
  library(kableExtra)

  # Applica il test
  out <- car::outlierTest(model)

  # Se nessun outlier rilevato
  if (is.null(out)) {
    return(HTML("<p><strong>Nessun outlier significativo rilevato dal test di Bonferroni.</strong></p>"))
  }

  # Estrai valori (è una named vector)
  rstudent <- as.numeric(out$rstudent)
  pval <- as.numeric(out$bonf.p)
  nomi <- names(out$rstudent)

  # Costruisci data frame manualmente
  df <- data.frame(
    osservazione = as.numeric(nomi),
    Residuo_studentizzato = round(rstudent, 4),
    `p-value (Bonferroni)` = ifelse(pval < 0.05,
                                    paste0("<span style='color:black;'>", format(pval, digits = 4), "</span>"),
                                    format(pval, digits = 4)),
    check.names = FALSE
  )

  # Tabella formattata
  knitr::kable(df, format = "html", escape = FALSE,
               caption = "Osservazioni identificate come outlier dal test di Bonferroni") %>%
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                              full_width = TRUE, position = "center") %>%
    kableExtra::row_spec(0, bold = TRUE)
}

tabella_outlier_html(mod_bic_revisioned)

Osservazioni identificate come outlier dal test di Bonferroni
osservazione	Residuo_studentizzato	p-value (Bonferroni)
1551	10.1463	2.473e-20
155	5.0897	9.636e-04
1306	4.7773	4.700e-03

Le osservazioni 1551, 155, e 1306 presentano residui studentizzati molto elevati.

I rispettivi p-value corretti con Bonferroni sono estremamente piccoli, il che significa che la probabilità che tali scostamenti siano casuali è molto bassa.

# Outlier ID
outlier_id <- c(1551, 155, 1306)

# Tabella pulita senza colonne duplicate
df_outlier <- df[outlier_id, ] %>%
  dplyr::mutate(osservazione = outlier_id) %>%
  dplyr::relocate(osservazione) %>%
  dplyr::select(-Gestazione2, -Anni.madre2, -N.gravidanze2)

# Tabella HTML senza rownames
knitr::kable(df_outlier, format = "html", row.names = FALSE,
             caption = "Osservazioni outlier con residui anomali (test di Bonferroni)") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                            full_width = TRUE, position = "center") %>%
  kableExtra::row_spec(0, bold = TRUE)

Osservazioni outlier con residui anomali (test di Bonferroni)
osservazione	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Sesso
1551	35	1	38	4370	315	374	0
155	30	0	36	3610	410	330	1
1306	23	0	41	4900	510	352	0

Calcoliamo ora la distanza di Cook dopo aver individuato gli outlier in modo da distinguere outlier influenti da outlier innocui, e poter decidere se tenerli, trattarli o rimuoverli.

# Calcolo della distanza di Cook
cook <- cooks.distance(mod_bic_revisioned)

# Soglia di riferimento (4 / n è un criterio comunemente usato)
soglia <- 4 / length(cook)

# Crea il dataframe per ggplot
df_cook <- data.frame(
  osservazione = 1:length(cook),
  cook = cook
)

# Grafico
ggplot(df_cook, aes(x = osservazione, y = cook)) +
  geom_point(color = "steelblue", size = 1.5) +
  geom_hline(yintercept = 0.5, color = "darkred", linetype = "dashed", size = 1) +
  annotate("text", x = max(df_cook$osservazione) * 0.95, y = 0.52,
           label = "Soglia = 0.5", color = "darkred", hjust = 0) +
  coord_cartesian(ylim = c(0, max(cook) * 1.1)) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Distanza di Cook per osservazione",
    x = "Indice osservazione",
    y = "Cook’s Distance"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Il grafico mostra una sola osservazione oltre la soglia di 0.5, su cui porrò la mia attenzione in quanto l’unica realmente anomala nel dataset, vediamo di quale osservazione si tratta.

# Crea tabella pulita senza rownames
cook_df <- data.frame(
  osservazione = which(cook > 0.5),
  `Cook’s Distance` = round(as.numeric(cook[cook > 0.5]), 5),
  check.names = FALSE
)

# Tabella HTML finale
kbl(cook_df, format = "html", row.names = FALSE,
    caption = "Osservazioni con Cook’s Distance > 0.5") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, position = "center") %>%
  row_spec(0, bold = TRUE)

Osservazioni con Cook’s Distance > 0.5
osservazione	Cook’s Distance
1551	0.65388

L’osservazione 1551 presenta un peso e un diametro cranico superiori alla media, ed è influente sul modello secondo la distanza di Cook e il test di Bonferroni. Tuttavia, trattandosi di un caso clinicamente plausibile di macrosomia neonatale ⁽⁷⁾, decido di mantenerla nell’analisi.

In conclusione

Nonostante le violazioni delle assunzioni classiche (normalità ed eteroscedasticità dei residui), il modello mod_bic_revisioned mostra un comportamento stabile in termini di indipendenza dei residui e conserva una buona capacità predittiva. In contesti applicativi come quello clinico, la previsione accurata è spesso più rilevante della perfetta aderenza ai presupposti teorici. Pertanto, il modello può essere considerato valido per fini predittivi.

3. Previsioni e risultati

Ora userò il modello mod_bic_revisioned per predire il peso di una neonata, da una madre alla terza gravidanza e 39 settimane di gestazione.

# Calcolo valori mancanti stimati
media_lunghezza <- round(mean(df$Lunghezza, na.rm = TRUE), 0)
media_cranio <- round(mean(df$Cranio, na.rm = TRUE), 0)

# Crea tabella descrittiva
tabella_variabili <- data.frame(
  Variabile = c("Gestazione", "Gestazione2", "N.gravidanze", "Lunghezza", "Cranio", "Sesso", "Fumatrici"),
  Tipo = c("Numerica (intera)", "Derivata", "Intera", "Numerica", "Numerica", "Binaria (0 = F)", "Binaria (0/1)"),
  Obbligatoria = rep("✅ sì", 7),
  `Come ricavarla se mancante` = c(
    "Fornita nel caso d’esempio (39)",
    "Calcolata come Gestazione^2 = 1521",
    "Fornita (es. 3)",
    paste0("Usata media del dataset (", media_lunghezza, " mm)"),
    paste0("Usata media del dataset (", media_cranio, " mm)"),
    "Indicata nel caso: femmina = 0",
    "Assunta come 0 (non fumatrice)"
  ),
  check.names = FALSE
)

# Tabella HTML
kbl(tabella_variabili, format = "html", escape = FALSE,
    caption = "Variabili richieste per la previsione del peso neonatale") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, position = "center") %>%
  row_spec(0, bold = TRUE)

Variabili richieste per la previsione del peso neonatale
Variabile	Tipo	Obbligatoria	Come ricavarla se mancante
Gestazione	Numerica (intera)	✅ sì	Fornita nel caso d’esempio (39)
Gestazione2	Derivata	✅ sì	Calcolata come Gestazione^2 = 1521
N.gravidanze	Intera	✅ sì	Fornita (es. 3)
Lunghezza	Numerica	✅ sì	Usata media del dataset (495 mm)
Cranio	Numerica	✅ sì	Usata media del dataset (340 mm)
Sesso	Binaria (0 = F)	✅ sì	Indicata nel caso: femmina = 0
Fumatrici	Binaria (0/1)	✅ sì	Assunta come 0 (non fumatrice)

Eseguo la previsione con i valori definiti nella tabella sovrastante.

nuovo_caso <- data.frame(
  Gestazione = 39,
  Gestazione2 = 39^2,
  N.gravidanze = 3,
  Lunghezza = media_lunghezza,
  Cranio = media_cranio,
  Sesso = 0,        # femmina
  Fumatrici = 0     # se non indicato, assumiamo non fumatrice
)

# Risultato della previsione
risultato <- predict(mod_bic_revisioned, newdata = nuovo_caso, interval = "confidence")

# Convertilo in data frame formattato
df_prev <- data.frame(
  `Stima del peso (g)` = round(risultato[1, "fit"], 0),
  `Limite inferiore (g)` = round(risultato[1, "lwr"], 0),
  `Limite superiore (g)` = round(risultato[1, "upr"], 0)
)

# Tabella HTML
kbl(df_prev, format = "html", caption = "Stima del peso neonatale con intervallo di confidenza al 95%") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = TRUE, position = "center") %>%
  row_spec(0, bold = TRUE)

Stima del peso neonatale con intervallo di confidenza al 95%
Stima.del.peso..g.	Limite.inferiore..g.	Limite.superiore..g.
3272	3248	3295

La stima del peso alla nascita per una neonata, figlia di una madre alla terza gravidanza che partorirà alla 39ª settimana, è pari a 3272 grammi, con un intervallo di confidenza al 95% compreso tra 3248 e 3295 grammi.

4. Visualizzazioni

# Converte sesso numerico in fattore con etichette M/F
df$SessoFatt <- factor(df$Sesso, levels = c(1, 0), labels = c("M", "F"))

# Grafico
ggplot(data = df, aes(x = Gestazione, y = Peso)) +
  geom_jitter(aes(color = SessoFatt), width = 0.3, height = 0, size = 1.5, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE, color = "black", size = 1) +
  labs(
    title = "Andamento del peso al variare delle settimane di gestazione per genere",
    x = "Settimane di gestazione",
    y = "Peso alla nascita (g)",
    color = "Sesso"
  ) +
  scale_color_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "top"
  )

Il grafico mostra una chiara relazione crescente tra le settimane di gestazione e il peso alla nascita. Tale andamento è coerente tra maschi e femmine, suggerendo che la durata della gestazione incide positivamente sul peso indipendentemente dal genere.

# Modello lineare semplice: Peso ~ Gestazione
mod_gest <- lm(Peso ~ Gestazione, data = df)

# Estrai i coefficienti e l'intervallo di confidenza
stima <- summary(mod_gest)$coefficients["Gestazione", ]
ci <- confint(mod_gest, level = 0.95)["Gestazione", ]

# Crea data frame per la tabella
tabella <- data.frame(
  `Stima incremento medio (g/settimana)` = round(stima["Estimate"], 2),
  `Errore standard` = round(stima["Std. Error"], 2),
  `p-value` = format.pval(stima["Pr(>|t|)"], digits = 3, eps = 1e-4),
  `CI 2.5%` = round(ci[1], 2),
  `CI 97.5%` = round(ci[2], 2)
)

# Visualizza tabella formattata
knitr::kable(tabella, format = "html", digits = 2,
             caption = "Stima dell'incremento del peso neonatale per ogni settimana di gestazione") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                            full_width = TRUE, position = "center") %>%
  kableExtra::row_spec(0, bold = TRUE)

Stima dell’incremento del peso neonatale per ogni settimana di gestazione
	Stima.incremento.medio..g.settimana.	Errore.standard	p.value	CI.2.5.	CI.97.5.
Estimate	166.27	4.53	<1e-04	157.39	175.16

Abbiamo un incremento medio del peso di circa 166g per ogni settima di gestione.

# Converte in fattore con etichette leggibili
df$FumatriciFatt <- factor(df$Fumatrici, levels = c(0, 1),
                           labels = c("Non fumatrici", "Fumatrici"))

# Grafico
ggplot(data = df, aes(x = Gestazione, y = Peso)) +
  geom_jitter(aes(color = FumatriciFatt), width = 0.3, height = 0, size = 1.5, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE, color = "black", size = 1) +
  labs(
    title = "Andamento del peso neonatale in funzione della durata della gestazione tra fumatrici e non fumatrici",
    x = "Settimane di gestazione",
    y = "Peso alla nascita (g)",
    color = "Abitudine al fumo"
  ) +
  scale_color_manual(values = c("Fumatrici" = "black", "Non fumatrici" = "#1E90FF")) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "top"
  )

Il grafico suggerisce che i neonati di madri fumatrici tendono ad avere un peso alla nascita inferiore rispetto a quelli di non fumatrici a parità di settimane di gestazione, supportando l’inclusione di questa variabile nel modello predittivo.

# Calcolo valori medi
media_lunghezza <- mean(df$Lunghezza, na.rm = TRUE)
media_cranio <- mean(df$Cranio, na.rm = TRUE)

# Nuovi casi: fumatrici vs non fumatrici
nuovi_casi <- data.frame(
  Gestazione = 39,
  Gestazione2 = 39^2,
  N.gravidanze = 3,
  Lunghezza = media_lunghezza,
  Cranio = media_cranio,
  Sesso = 0,  # femmina
  Fumatrici = c(0, 1)  # non fumatrice, fumatrice
)

# Previsione del peso con intervallo di confidenza
predizioni <- predict(mod_bic_revisioned, newdata = nuovi_casi, interval = "confidence")

# Tabella dei risultati
tabella <- data.frame(
  `Abitudine al fumo` = c("Non fumatrice", "Fumatrice"),
  `Stima del peso (g)` = round(predizioni[, "fit"]),
  `Limite inferiore (g)` = round(predizioni[, "lwr"]),
  `Limite superiore (g)` = round(predizioni[, "upr"])
)

# Visualizza in tabella HTML formattata
kbl(tabella, format = "html", escape = FALSE,
    caption = "Stima del peso neonatale per fumatrici e non fumatrici a 39 settimane") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, position = "center") %>%
  row_spec(0, bold = TRUE)

Stima del peso neonatale per fumatrici e non fumatrici a 39 settimane
Abitudine.al.fumo	Stima.del.peso..g.	Limite.inferiore..g.	Limite.superiore..g.
Non fumatrice	3269	3245	3293
Fumatrice	3240	3183	3296

ggplot(df, aes(x = Lunghezza, y = Peso, color = SessoFatt)) +
  geom_jitter(width = 0.3, alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, size = 1) +
  scale_color_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) +
  labs(
    title = "Relazione tra lunghezza e peso alla nascita per genere",
    x = "Lunghezza (cm)",
    y = "Peso alla nascita (g)",
    color = "Sesso"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Il grafico evidenzia una relazione fortemente positiva tra la lunghezza del neonato e il peso alla nascita, confermando l’elevata rilevanza predittiva di questa variabile all’interno del modello.

ggplot(df, aes(x = Peso, fill = SessoFatt)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("M" = "#1E90FF", "F" = "#FF69B4")) +
  labs(
    title = "Distribuzione del peso neonatale per genere",
    x = "Peso alla nascita (g)",
    y = "Densità",
    fill = "Sesso"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

Il grafico mostra che i neonati maschi tendono ad avere un peso alla nascita maggiore rispetto alle femmine, evidenziando un impatto significativo del sesso nella distribuzione del peso neonatale.

NOTE

⁽¹⁾ Impatto del fumo durante la gravidanza sul peso neonatale:

Ministero della salute: rischio di basso peso alla nascita come conseguenza della riduzione dell’apporto di ossigeno necessario per la corretta crescita del bambino.
[FONTE: https://www.salute.gov.it/portale/fumo/dettaglioContenutiFumo.jsp?lingua=italiano&id=5583&area=fumo&menu=vuoto].
Fondazione Alessandra Graziottin: il fumo materno in gravidanza correla significativamente con il rischio di basso peso alla nascita (OR = 1.89, 95% CI = 1.80-1.98):
- OR = 1.89: il fumo aumenta il rischio di basso peso alla nascita di circa 89% rispetto alle non fumatrici.
- CI 95% = 1.80 - 1.98: conferma che l’associazione è forte e statisticamente significativa.
- Poiché il CI non include 1, l’effetto del fumo sul peso del neonato è quasi sicuramente reale e non dovuto al caso.
- [FONTE: https://www.fondazionegraziottin.org/it/articolo.php/Gravidanza-il-fumo-aumenta-il-rischio-di-basso-peso-alla-nascita?EW_CHILD=29894]
Fondazione Umberto Veronesi: fumare in gravidanza non pone a rischio soltanto la salute materna, ma pure quella del nascituro che, nel lungo periodo, rischia di sviluppare una capacità aerobica inferiore a quella dei coetanei nati da donne non fumatrici.
[FONTE: https://www.fondazioneveronesi.it/magazine/articoli/ginecologia/cosi-il-fumo-gravidanza-danneggia-il-bambino]
Istituto Superiore di Sanità riporta che circa il 10% delle donne in gravidanza in Italia ha dichiarato di fumare. [FONTE: https://www.epicentro.iss.it/sorveglianza02anni/indagine-2022-risultati]

⁽²⁾ Valori medi di un neonato alla nascita:

ho trovato un documento del reparto pediatrico del Polincino Gemelli di Roma che si dimostra una fonte autorevole [FONTE: https://pspediatrico.policlinicogemelli.it/wp-content/uploads/pdfScaricabili-12-2017/2-Esame_obiettivo_neonatale.pdf] il quale riporta che in media un neonato normale ha:

un peso di 3100-3400 grammi
una lunghezza di 48-52 cm
una circonferenza cranica di 33-36 cm

il diametro del cranio dovrebbe essere tra 90 mm e 110 mm (9-11 cm), percui le misure presenti nel dataset si riferiscono alla circonferenza del cranio e non al diametro.

⁽³⁾ Impatto del tipo di parto sul peso neonatale:

Studi suggeriscono che il tipo di parto (naturale o cesareo) non è un determinante diretto del peso del neonato. Il cesareo viene spesso eseguito per motivi indipendenti dal peso (ad esempio, posizione del feto, sofferenza fetale, condizioni materne), il che rende questa variabile non causalmente legata al peso alla nascita.

Fonti:

Kuhle et al., Cesarean Delivery Is Not Associated with Obesity in Childhood, The Journal of Pediatrics (2015).
H. Blustein et al., Association of caesarean delivery with offspring overweight and obesity: a systematic review and meta-analysis, International Journal of Obesity (2013).

⁽⁴⁾ Impatto dell’ospedale sul peso neonatale:

Il peso neonatale è più influenzato da fattori materni e fetali (età della madre, durata della gestazione, fumo, condizioni mediche, ecc.) piuttosto che dal luogo di nascita. Le differenze tra ospedali sono spesso legate a differenze nella popolazione trattata o nei protocolli di assistenza, ma non rappresentano un fattore determinante diretto.

Fonti:

Kogan et al., The Role of Clinical Practices in Variation in Birth Weight Among US Hospitals, Maternal and Child Health Journal (2010).
Kramer et al., Determinants of low birth weight: Methodological assessment and meta-analysis, Bulletin of the World Health Organization (1987).

⁽⁵⁾ Gravidanza in minorenni:

E’ possibile che una ragazza di 13 o 14 anni rimanga incinta. La fertilità femminile inizia con la prima mestruazione, che in media si verifica intorno ai 12 anni. [FONTE: https://www.ospedalebambinogesu.it/gravidanza-in-minorenni-96849]

⁽⁶⁾ Gravidanza oltre i 40 anni:

Diventare mamma a 40 anni non è affatto strano, ma sta quasi diventando la normalità. [FONTE:https://imamma.it/gravidanza-a-40-anni-tutte-le-cose-da-sapere/?utm_source=chatgpt.com]

⁽⁷⁾ Macrosomia neonatale:

La macrosomia neonatale è definita come un peso alla nascita superiore a 4000 g. [FONTE:https://www.msdmanuals.com/it/professionale/pediatria/problemi-perinatali/neonato-grande-per-l-et%C3%A0-gestazionale?utm_source=chatgpt.com]
In alcuni casi, si parla di macrosomia severa quando il peso supera i 4500 g. [FONTE: https://it.wikipedia.org/wiki/Macrosomia_fetale?utm_source=chatgpt.com]
Un diametro cranico di 374 mm (equivalente a una circonferenza di circa 37,4 cm) è superiore alla media, ma può essere giustificato in un neonato macrosomico essendo di poco superiore alla media ⁽²⁾.

Nella redazione del presente documento si è fatto uso di strumenti di assistenza alla scrittura basati su intelligenza artificiale, come ChatGPT, con l’obiettivo di migliorare la chiarezza espositiva, la coerenza linguistica e la resa visiva. Le decisioni sui contenuti e l’analisi scientifica sono state curate in modo autonomo e consapevole dall’autore.

Modello statistico per la previsione del peso neonatale

Dario Colombo

1. Raccolta dei dati e struttura del dataset

Dimensioni del Dataset: 2500 righe x 10 colonne

Anni.madre

N.gravidanze

Fumatrici (1)

Gestazione

Peso => VARIABILE TARGET

Lunghezza

Cranio (2)

Tipo.parto (3)

Ospedale (4)

Sesso

2. Analisi e modellizzazione

Analisi preliminare

Analisi di “Anni.madre”

Analisi di “N.gravidanze”

Analisi di “Fumatrici”

Analisi di “Gestazione”

Analisi di “Peso”

Analisi di “Lunghezza”

Analisi di “Cranio”

Analisi di “Sesso”

In alcuni ospedali si fanno più parti cesarei

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

Le misure antropometriche sono significativamente diverse tra i due sessi

Creazione del Modello di Regressione

Analisi della correlazione tra le variabili

Modello 1

Residui del Modello 1

Coefficienti del Modello 1

Statistiche del Modello 1

Variance Inflation Factor (VIF) per valutare la multicollinearità

Modello 2

Residui del Modello 2

Coefficienti del Modello 2

Statistiche del Modello 2

Variance Inflation Factor (VIF) per valutare la multicollinearità

Modello 3

Residui del Modello 3

Coefficienti del Modello 3

Statistiche del Modello 3

Variance Inflation Factor (VIF) per valutare la multicollinearità

Osservazioni sui modelli analizzati fino ad ora

Modello 4

Residui del Modello 4

Coefficienti del Modello 4

Statistiche del Modello 4

Variance Inflation Factor (VIF) per valutare la multicollinearità

Modello 5

Residui del Modello 5

Coefficienti del Modello 5

Statistiche del Modello 5

Variance Inflation Factor (VIF) per valutare la multicollinearità

Modello 6

Residui del Modello 6

Coefficienti del Modello 6

Statistiche del Modello 6

Variance Inflation Factor (VIF) per valutare la multicollinearità

Selezione del modello ottimale con stepwise

Formula del modello selezionato con stepwise

Residui del Modello selezionato con stepwise

Coefficienti del Modello selezionato con stepwise

Statistiche del Modello selezionato con stepwise

Correzione del Modello selezionato con stepwise

Residui del Modello selezionato con stepwise (Corretto)

Coefficienti del Modello selezionato con stepwise (Corretto)

Statistiche del Modello selezionato con stepwise (Corretto)

Residui del Modello selezionato con stepwise (Rivisto)

Coefficienti del Modello selezionato con stepwise (Rivisto)

Statistiche del Modello selezionato con stepwise (Rivisto)

Confronto tra modelli: BIC e R² Adjusted

Analisi della qualità del modello

3. Previsioni e risultati

4. Visualizzazioni

NOTE

Fumatrici ⁽¹⁾

Cranio ⁽²⁾

Tipo.parto ⁽³⁾

Ospedale ⁽⁴⁾