library(ggplot2)
library(dplyr)
library(corrplot)
library(kableExtra)
library(car)
library(caret)
data <- read.csv("neonati.csv", stringsAsFactors = FALSE)

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

print(na_count)
##   Anni.madre N.gravidanze    Fumatrici   Gestazione         Peso    Lunghezza 
##            0            0            0            0            0            0 
##       Cranio   Tipo.parto     Ospedale        Sesso 
##            0            0            0            0
data <- data %>% distinct()

data <- data %>% filter(Anni.madre > 10)

data$Fumatrici <- as.factor(data$Fumatrici)
data$Tipo.parto <- as.factor(data$Tipo.parto)
data$Ospedale <- as.factor(data$Ospedale)
data$Sesso <- as.factor(data$Sesso)

kable(summary(data))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
Min. :13.00 Min. : 0.0000 0:2394 Min. :25.00 Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1255
1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1770 osp2:848 M:1243
Median :28.00 Median : 1.0000 NA Median :39.00 Median :3300 Median :500.0 Median :340 NA osp3:834 NA
Mean :28.19 Mean : 0.9816 NA Mean :38.98 Mean :3284 Mean :494.7 Mean :340 NA NA NA
3rd Qu.:32.00 3rd Qu.: 1.0000 NA 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350 NA NA NA
Max. :46.00 Max. :12.0000 NA Max. :43.00 Max. :4930 Max. :565.0 Max. :390 NA NA NA

Nella prima fase del progetto ho importato il dataset e verificato la sua struttura. In seguito ho cercato eventuali dati mancanti, eliminato eventuali dupplicati e i dati per le madri con età inferiore a 10 anni, assumendo che sia un errore di registrazione dei dati. Infine ho convertito in fattori le variabili categoriche. Il dataset è composto da 2498 registrazioni.

par(mfrow = c(2, 3))
hist(data$Anni.madre, main = "Età Madre", xlab = "Età (anni)", col = "lightblue", border = "black")
hist(data$N.gravidanze, main = "Numero Gravidanze", xlab = "Numero", col = "lightgreen", border = "black")
hist(data$Gestazione, main = "Settimane di Gestazione", xlab = "Settimane", col = "lightpink", border = "black")
hist(data$Peso, main = "Peso Neonatale", xlab = "Peso (g)", col = "lightcoral", border = "black")
hist(data$Lunghezza, main = "Lunghezza", xlab = "Lunghezza (mm)", col = "lightyellow", border = "black")
hist(data$Cranio, main = "Diametro Cranio", xlab = "Diametro (mm)", col = "lightgray", border = "black")

I grafici a istogramma mostrano la distribuzione delle variabili continue. La maggior parte delle variabili, come l’età materna e la durata della gestazione, presenta una distribuzione concentrata intorno alla media, con qualche variabilità. Tuttavia, il peso neonatale evidenzia una distribuzione leggermente asimmetrica con possibili outlier visibili ai valori estremamente bassi e alti

par(mfrow = c(1, 3))
boxplot(data$Peso, main = "Peso", col = "lightblue")
boxplot(data$Lunghezza, main = "Lunghezza", col = "lightgreen")
boxplot(data$Cranio, main = "Diametro Cranio", col = "lightpink")

### Il boxplot del peso neonatale evidenzia la presenza di alcuni outlier sia nella parte inferiore che superiore. Questi valori estremi potrebbero essere casi di nascite premature o neonati macrosomici. Anche le variabili lunghezza e cranio mostrano qualche outlier, ma in misura minore rispetto al peso.

par(mfrow = c(1, 2))
barplot(table(data$Tipo.parto), main = "Tipo di Parto", col = c("lightblue", "lightcoral"))
barplot(table(data$Ospedale), main = "Distribuzione per Ospedale", col = c("lightblue", "lightgreen", "lightpink"))

I grafici a barre ci mostrano una marcata disparità tra il parto nat e il parto ces in favore del primo. La distribuzione degli ospedali invece è equilibrata.

cor_matrix <- cor(data[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")])
corrplot(cor_matrix, method = "color")

### La matrice di correlazione evidenzia relazioni lineari positive tra il peso neonatale e variabili come la lunghezza e il diametro craniale, suggerendo che queste variabili antropometriche potrebbero essere predittori rilevanti. L’età materna e il numero di gravidanze mostrano correlazioni più deboli, suggerendo che il loro impatto diretto sul peso potrebbe essere minore rispetto ad altri fattori.

par(mfrow = c(1, 2))
boxplot(Peso ~ Fumatrici, data = data, main = "Peso e Fumo Materno", col = c("lightblue", "lightcoral"))
plot(data$Gestazione, data$Peso, main = "Peso vs Gestazione", xlab = "Settimane di Gestazione", ylab = "Peso (g)", col = "blue", pch = 16)
abline(lm(Peso ~ Gestazione, data = data), col = "red")

### Fumo materno: I boxplot indicano una differenza visibile nel peso tra neonati di madri fumatrici e non fumatrici. I bambini di madri fumatrici tendono ad avere un peso inferiore, supportando l’ipotesi di un impatto negativo del fumo sulla crescita fetale.

Settimane di gestazione: Il grafico di dispersione mostra una chiara tendenza positiva tra settimane di gestazione e peso alla nascita. La linea di regressione suggerisce che un prolungamento della gestazione è associato a un aumento del peso neonatale, confermando l’importanza di questa variabile nel modello predittivo.

cont_table_ospedale_parto <- table(data$Ospedale, data$Tipo.parto)

test_chi_ospedale_parto <- chisq.test(cont_table_ospedale_parto)

chi_square_results <- data.frame(
  Statistic = round(test_chi_ospedale_parto$statistic, 2),
  P_value = round(test_chi_ospedale_parto$p.value, 2),
  DF = test_chi_ospedale_parto$parameter
)

kable(chi_square_results, 
      col.names = c("Chi-squared Statistic", "P-value", "Degrees of Freedom"), 
      caption = "Risultati del Test Chi-quadrato")
Risultati del Test Chi-quadrato
Chi-squared Statistic P-value Degrees of Freedom
X-squared 1.08 0.58 2

Poiché il p-value è maggiore di 0.05, non rifiutiamo l’ipotesi nulla. Questo significa che non ci sono prove sufficienti per affermare che ci siano differenze significative nel numero di parti cesarei tra i diversi ospedali.

In base ai dati a disposizione, possiamo concludere che la distribuzione dei parti cesarei sembra indipendente dall’ospedale. Non ci sono evidenze che in alcuni ospedali si facciano significativamente più parti cesarei rispetto ad altri.

media_peso_campione <- mean(data$Peso)
media_lunghezza_campione <- mean(data$Lunghezza)

media_peso_popolazione <- 3300 #dato preso da un articolo dell'ospedale pediatrico bambino gesù
media_lunghezza_popolazione <- 50  #dato preso da un articolo dell'ospeale pediatrico bambino gesù

test_peso <- t.test(data$Peso, mu = media_peso_popolazione)
test_lunghezza <- t.test(data$Lunghezza, mu = media_lunghezza_popolazione)

results <- data.frame(
  Test = c("Peso", "Lunghezza"),
  Statistic_t = round(c(test_peso$statistic, test_lunghezza$statistic), 2),
  P_value = round(c(test_peso$p.value, test_lunghezza$p.value), 2),
  CI_lower = round(c(test_peso$conf.int[1], test_lunghezza$conf.int[1]), 2),
  CI_upper = round(c(test_peso$conf.int[2], test_lunghezza$conf.int[2]), 2)
)

kable(results, col.names = c("Test", "Statistic t", "P-value", "CI Lower", "CI Upper"), 
      caption = "Risultati dei Test t")
Risultati dei Test t
Test Statistic t P-value CI Lower CI Upper
Peso -1.51 0.13 3263.58 3304.79
Lunghezza 844.17 0.00 493.66 495.73

Per il peso: Non c’è una differenza significativa rispetto alla popolazione.

Per la lunghezza: La differenza rispetto alla popolazione è altamente significativa.

test_peso <- t.test(data$Peso ~ data$Sesso)
test_lunghezza <- t.test(data$Lunghezza ~ data$Sesso)
test_cranio <- t.test(data$Cranio ~ data$Sesso)

t_test_results <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  `t-Statistic` = round(c(test_peso$statistic, test_lunghezza$statistic, test_cranio$statistic), 2),
  `P-value` = round(c(test_peso$p.value, test_lunghezza$p.value, test_cranio$p.value), 2),
  `Gradi di Libertà` = c(test_peso$parameter, test_lunghezza$parameter, test_cranio$parameter)
)

kable(t_test_results, 
      col.names = c("Variabile", "t-Statistic", "P-value", "Gradi di Libertà"), 
      caption = "Risultati dei Test t per le Variabili Peso, Lunghezza e Cranio in relazione al Sesso")
Risultati dei Test t per le Variabili Peso, Lunghezza e Cranio in relazione al Sesso
Variabile t-Statistic P-value Gradi di Libertà
Peso -12.11 0 2488.674
Lunghezza -9.58 0 2457.301
Cranio -7.44 0 2489.389

In tutti e tre i test (Peso, Lunghezza, e Cranio), la differenza tra i sessi è statisticamente significativa, con le femmine che tendono ad avere valori medi inferiori per ciascuna delle misure rispetto ai maschi. Pertanto, possiamo rifiutare l’ipotesi nulla di assenza di differenze tra i gruppi per tutte e tre le variabili.

data_standardizzato <- data
data_standardizzato[c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")] <- 
  scale(data[c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")])
set.seed(123)

index <- createDataPartition(data_standardizzato$Peso, p = 0.8, list = FALSE)

trainData <- data_standardizzato[index, ]
testData <- data_standardizzato[-index, ]

kable(dim(trainData))
x
2001
10 
kable(dim(testData))
x
497
10 
modello <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio, data = trainData)

summary_tab <- summary(modello)$coefficients %>% 
  as.data.frame() %>%
  mutate(Termini = rownames(.)) %>%
  select(Termini, Estimate, `Std. Error`, `t value`, `Pr(>|t|)`)

kable(summary_tab, 
      digits = 3, 
      col.names = c("Termini", "Coefficiente", "Errore Std.", "t-Statistic", "P-value"),
      caption = "Risultati del Modello Lineare")
Risultati del Modello Lineare
Termini Coefficiente Errore Std. t-Statistic P-value
(Intercept) (Intercept) 3289.900 6.422 512.252 0.000
Anni.madre Anni.madre 7.111 6.845 1.039 0.299
N.gravidanze N.gravidanze 10.676 6.822 1.565 0.118
Fumatrici1 Fumatrici1 -33.468 30.929 -1.082 0.279
Gestazione Gestazione 65.871 8.109 8.124 0.000
Lunghezza Lunghezza 265.047 9.004 29.436 0.000
Cranio Cranio 177.659 8.000 22.206 0.000 
predizioni <- predict(modello, newdata = testData)
residui_test <- testData$Peso - predizioni
RMSE_test <- sqrt(mean(residui_test^2))
R2_test <- 1 - sum(residui_test^2) / sum((testData$Peso - mean(testData$Peso))^2)

performance_metrics <- data.frame(
  Metrica = c("R²", "RMSE"),
  Valore = c(round(R2_test, 3), round(RMSE_test, 2))
)

kable(performance_metrics, 
      digits = 3,
      col.names = c("Metrica", "Valore"),
      caption = "Metriche di Performance del Modello")
Metriche di Performance del Modello
Metrica Valore
R² 0.742
RMSE 263.070 
modello <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio, data = trainData)
modello2 <- lm(Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio, data = trainData)
modello3 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio, data = trainData)

anova_res <- anova(modello, modello2, modello3)

BIC_modello1 <- BIC(modello)
BIC_modello2 <- BIC(modello2)
BIC_modello3 <- BIC(modello3)

BIC_values <- data.frame(
  Modello = c("Modello 1", "Modello 2", "Modello 3"),
  BIC = c(round(BIC_modello1, 2), round(BIC_modello2, 2), round(BIC_modello3, 2))
)

kable(BIC_values, 
      digits = 2, 
      col.names = c("Modello", "BIC"), 
      caption = "BIC dei Modelli")
BIC dei Modelli
Modello BIC
Modello 1 28296.48
Modello 2 28289.97
Modello 3 28283.55 
predizioni <- predict(modello3, newdata = testData)

residui_test <- testData$Peso - predizioni
RMSE_test <- sqrt(mean(residui_test^2))
R2_test <- 1 - sum(residui_test^2) / sum((testData$Peso - mean(testData$Peso))^2)

performance <- data.frame(
  Metric = c("R²", "RMSE"),
  Value = c(round(R2_test, 3), round(RMSE_test, 2))
)

kable(performance, col.names = c("Metric", "Value"), caption = "Performance sul Test Set")
Performance sul Test Set
Metric Value
R² 0.743
RMSE 262.800 
modello_interazioni <- lm(Peso ~ (N.gravidanze + Gestazione + Lunghezza + Cranio)^2, data = trainData)

anova_res_interazioni <- anova(modello3, modello_interazioni)

BIC_modello3 <- BIC(modello3)
BIC_modello_interazioni <- BIC(modello_interazioni)

predizioni_interazioni <- predict(modello_interazioni, newdata = testData)
residui_test_interazioni <- testData$Peso - predizioni_interazioni
RMSE_interazioni <- sqrt(mean(residui_test_interazioni^2))
R2_interazioni <- 1 - sum(residui_test_interazioni^2) / sum((testData$Peso - mean(testData$Peso))^2)

BIC_values <- data.frame(
  Modello = c("Modello 3", "Modello con interazioni"),
  BIC = c(round(BIC_modello3, 2), round(BIC_modello_interazioni, 2))
)

performance_test <- data.frame(
  Metrica = c("R²", "RMSE"),
  Valore = c(round(R2_interazioni, 3), round(RMSE_interazioni, 2))
)

kable(BIC_values, 
      digits = 2, 
      col.names = c("Modello", "BIC"), 
      caption = "Confronto BIC tra i modelli")
Confronto BIC tra i modelli
Modello BIC
Modello 3 28283.55
Modello con interazioni 28281.46 
kable(performance_test, 
      digits = 2, 
      col.names = c("Metrica", "Valore"), 
      caption = "Performance del modello con interazioni sul test set")
Performance del modello con interazioni sul test set
Metrica Valore
R² 0.74
RMSE 262.46 
modello_interazioni_ridotto <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Gestazione:Cranio, data = trainData)

anova_res_interazioni_ridotto <- anova(modello_interazioni, modello_interazioni_ridotto)

BIC_modello_interazioni_ridotto <- BIC(modello_interazioni_ridotto)

predizioni_interazioni_ridotto <- predict(modello_interazioni_ridotto, newdata = testData)
residui_test_interazioni_ridotto <- testData$Peso - predizioni_interazioni_ridotto
RMSE_interazioni_ridotto <- sqrt(mean(residui_test_interazioni_ridotto^2))
R2_interazioni_ridotto <- 1 - sum(residui_test_interazioni_ridotto^2) / sum((testData$Peso - mean(testData$Peso))^2)

BIC_values_ridotto <- data.frame(
  Modello = c("Modello con interazioni", "Modello con interazioni ridotte"),
  BIC = c(round(BIC_modello_interazioni, 2), round(BIC_modello_interazioni_ridotto, 2))
)

performance_test_ridotto <- data.frame(
  Metrica = c("R²", "RMSE"),
  Valore = c(round(R2_interazioni_ridotto, 3), round(RMSE_interazioni_ridotto, 2))
)

kable(BIC_values_ridotto, 
      digits = 2, 
      col.names = c("Modello", "BIC"), 
      caption = "Confronto BIC tra il modello completo e quello ridotto")
Confronto BIC tra il modello completo e quello ridotto
Modello BIC
Modello con interazioni 28281.46
Modello con interazioni ridotte 28250.78 
kable(performance_test_ridotto, 
      digits = 2, 
      col.names = c("Metrica", "Valore"), 
      caption = "Performance del modello con interazioni ridotte sul test set")
Performance del modello con interazioni ridotte sul test set
Metrica Valore
R² 0.74
RMSE 262.63

Dopo aver standardizzato le variabili continue del dataset e diviso lo stesso in train e test set. Ho creato vari modelli per analizzarne le prestazioni e decidere col quale proseguire per l’analisi degli outliers. Le prestazioni migliori in base ai dati raccolti con il BIC e l’anova vengono fornite dal modello_interazioni_ridotto.

shapiro_test <- shapiro.test(residuals(modello_interazioni_ridotto))

shapiro_statistic <- round(shapiro_test$statistic, 3)
shapiro_p_value <- (shapiro_test$p.value)

shapiro_results <- data.frame(
  Statistica = shapiro_statistic,
  P_value = shapiro_p_value
)

kable(shapiro_results, 
      digits = 3, 
      col.names = c("Statistica", "P-value"), 
      caption = "Risultati del test di Shapiro-Wilk sui residui del modello con interazioni ridotte")
Risultati del test di Shapiro-Wilk sui residui del modello con interazioni ridotte
Statistica P-value
W 0.972 0

Dal p value ottenuto col shapiro test possiamo affermare che i residui non seguano una distribuzione normale.

plot(fitted(modello_interazioni_ridotto), residuals(modello_interazioni_ridotto))
abline(h = 0, col = "red")

boxplot(residuals(modello_interazioni_ridotto), main = "Boxplot dei residui")

cooks.distance <- cooks.distance(modello_interazioni_ridotto)
plot(cooks.distance, type = "h", main = "Distanze di Cook")
abline(h = 4 / length(cooks.distance), col = "red")

outliers_influenti <- which(cooks.distance(modello_interazioni_ridotto) > (4 / length(cooks.distance(modello_interazioni_ridotto))))
data_senza_outliers <- trainData[-outliers_influenti, ]

modello_senza_outliers <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Gestazione:Cranio, data = data_senza_outliers)

summary_modello_senza_outliers <- summary(modello_senza_outliers)
summary_table <- data.frame(
  Coefficienti = rownames(summary_modello_senza_outliers$coefficients),
  Estimati = round(summary_modello_senza_outliers$coefficients[, 1], 3),
  Errori_Standard = round(summary_modello_senza_outliers$coefficients[, 2], 3),
  Statistiche_t = round(summary_modello_senza_outliers$coefficients[, 3], 3),
  P_value = round(summary_modello_senza_outliers$coefficients[, 4], 3)
)

kable(summary_table, caption = "Sommario del modello senza outliers", digits = 3)
Sommario del modello senza outliers
Coefficienti Estimati Errori_Standard Statistiche_t P_value
(Intercept) (Intercept) 3270.944 5.735 570.365 0.000
N.gravidanze N.gravidanze 18.642 6.221 2.996 0.003
Gestazione Gestazione 76.562 7.983 9.590 0.000
Lunghezza Lunghezza 301.516 8.661 34.812 0.000
Cranio Cranio 173.870 7.577 22.947 0.000
Gestazione:Cranio Gestazione:Cranio 19.751 3.467 5.697 0.000 
shapiro_test_senza_outliers <- shapiro.test(residuals(modello_senza_outliers))

shapiro_statistic_senza_outliers <- round(shapiro_test_senza_outliers$statistic, 3)
shapiro_p_value_senza_outliers <- round(shapiro_test_senza_outliers$p.value, 3)

shapiro_results_senza_outliers <- data.frame(
  Statistica = shapiro_statistic_senza_outliers,
  P_value = shapiro_p_value_senza_outliers
)

kable(shapiro_results_senza_outliers, 
      digits = 3, 
      col.names = c("Statistica", "P-value"), 
      caption = "Risultati del test di Shapiro-Wilk sui residui del modello senza outliers")
Risultati del test di Shapiro-Wilk sui residui del modello senza outliers
Statistica P-value
W 0.999 0.259

Grazie al calcolo degli outliers influenti ho creato un nuovo modello dal quale ho eliminato quest’ultimi. In questo modo ho ottenuto un risultato soddisfacente per quanto riguarda il p-value del shapiro test, i residui del nuovo modello seguono una distribuzione normale.

modelloFinale <- step(modello_senza_outliers, direction = "backward")
## Start:  AIC=21078.89
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Gestazione:Cranio
## 
##                     Df Sum of Sq       RSS   AIC
## <none>                           112404776 21079
## - N.gravidanze       1    527554 112932330 21086
## - Gestazione:Cranio  1   1906769 114311545 21109
## - Lunghezza          1  71206856 183611631 22019
summary_modello_finale <- summary(modelloFinale)
summary_table_finale <- data.frame(
  Coefficienti = rownames(summary_modello_finale$coefficients),
  Estimati = round(summary_modello_finale$coefficients[, 1], 3),
  Errori_Standard = round(summary_modello_finale$coefficients[, 2], 3),
  Statistiche_t = round(summary_modello_finale$coefficients[, 3], 3),
  P_value = round(summary_modello_finale$coefficients[, 4], 3)
)

kable(summary_table_finale, caption = "Sommario del modello finale", digits = 3)
Sommario del modello finale
Coefficienti Estimati Errori_Standard Statistiche_t P_value
(Intercept) (Intercept) 3270.944 5.735 570.365 0.000
N.gravidanze N.gravidanze 18.642 6.221 2.996 0.003
Gestazione Gestazione 76.562 7.983 9.590 0.000
Lunghezza Lunghezza 301.516 8.661 34.812 0.000
Cranio Cranio 173.870 7.577 22.947 0.000
Gestazione:Cranio Gestazione:Cranio 19.751 3.467 5.697 0.000 
shapiro_test_finale <- shapiro.test(residuals(modelloFinale))

shapiro_statistic_finale <- round(shapiro_test_finale$statistic, 3)
shapiro_p_value_finale <- round(shapiro_test_finale$p.value, 3)

shapiro_results_finale <- data.frame(
  Statistica = shapiro_statistic_finale,
  P_value = shapiro_p_value_finale
)

kable(shapiro_results_finale, 
      digits = 3, 
      col.names = c("Statistica", "P-value"), 
      caption = "Risultati del test di Shapiro-Wilk sui residui del modello finale")
Risultati del test di Shapiro-Wilk sui residui del modello finale
Statistica P-value
W 0.999 0.259 
predizioni <- predict(modelloFinale, newdata = testData)

residui_test <- testData$Peso - predizioni
RMSE_test <- sqrt(mean(residui_test^2))
R2_test <- 1 - sum(residui_test^2) / sum((testData$Peso - mean(testData$Peso))^2)

performance_results <- data.frame(
  R_squared = round(R2_test, 3),
  RMSE = round(RMSE_test, 2)
)

kable(performance_results, 
      digits = 3, 
      col.names = c("R²", "RMSE"), 
      caption = "Performance del modello finale")
Performance del modello finale
R² RMSE
0.747 260.27

Dopo aver creato il modello_senza_outliers ho applicato la tecnica del backward direction, eliminando progressivamente le variabili che non sono statisticamente significative. Il modelloFinale presenta solo le variabili n.gravidanze, gestazione, lunghezza e cranio più l’interazione gestazione:lunghezza. Proseguiremo con questo modello per le previsioni su nuovi dati.

L’R-squared è uguale a 0,7597, che indica che circa il 75,97% della variabilità nel log del peso è spiegata dal modello

Le prestazioni sul testSet danno risultati soddisfacenti con R² pari a 0,747 ed RMSE pari a 260,27.

Inoltre dal modello possiamo dedurre che: Un aumento di una unità in N.gravidanze aumenta il peso di 18.64g. Un aumento di una unità in Gestazione aumenta il peso di 76.56g. Un aumento di una unità in Lunghezza aumenta il peso di 301.51g. Un aumento di una unità in Cranio aumenta il peso di 173.87g.

new_data <- data.frame(
  Anni.madre = c(30, 28, 35, 29, 31),
  N.gravidanze = c(4, 2, 1, 0, 3),
  Gestazione = c(41, 41, 38, 37, 39), 
  Lunghezza = c(480, 493, 505, 485, 473), 
  Cranio = c(322, 341, 321, 353, 363)
)

new_data_standardizzato <- new_data
new_data_standardizzato[c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")] <- 
  scale(data[c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")])

new_predizioni <- predict(modelloFinale, newdata = new_data_standardizzato)

new_data$Peso_predetto <- new_predizioni

print(new_data)
##   Anni.madre N.gravidanze Gestazione Lunghezza Cranio Peso_predetto
## 1         30            4         41       480    322      3010.562
## 2         28            2         41       493    341      3660.250
## 3         35            1         38       505    321      3470.390
## 4         29            0         37       485    353      2708.672
## 5         31            3         39       473    363      3216.423

Ho utilizzato il modelloFinale per prevedere il peso di 5 nuvoi neonati.

predizioni <- data.frame(
  Gestazione = seq(min(data$Gestazione), max(data$Gestazione), length.out = 100),
  Lunghezza = mean(data$Lunghezza),  
  Cranio = mean(data$Cranio),         
  N.gravidanze = mean(data$N.gravidanze)  
)

predizioni$Peso_Previsto <- predict(modelloFinale, newdata = predizioni)

ggplot(predizioni, aes(x = Gestazione, y = Peso_Previsto)) +
  geom_line(color = "blue") +
  labs(title = "Impatto della Gestazione sul Peso", x = "Settimane di Gestazione") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

predizioni$Lunghezza <- seq(min(data$Lunghezza), max(data$Lunghezza), length.out = 100)

predizioni$Peso_Previsto_Lunghezza <- predict(modelloFinale, newdata = predizioni)

ggplot(predizioni, aes(x = Lunghezza, y = Peso_Previsto_Lunghezza)) +
  geom_line(color = "green") +
  labs(title = "Impatto della Lunghezza sul Peso", x = "Lunghezza del Feto") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

predizioni$N.gravidanze <- seq(min(data$N.gravidanze), max(data$N.gravidanze), length.out = 100)

predizioni$Peso_Previsto_Gravidanze <- predict(modelloFinale, newdata = predizioni)

ggplot(predizioni, aes(x = N.gravidanze, y = Peso_Previsto_Gravidanze)) +
  geom_line(color = "purple") +
  labs(title = "Impatto di N.gravidanze sul Peso", x = "Numero di Gravidanze") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

residui <- residuals(modelloFinale)

ggplot(data.frame(residui), aes(x = residui)) +
  geom_histogram(bins = 30, fill = "red", alpha = 0.7) +
  labs(title = "Distribuzione dei Residui", x = "Residui", y = "Frequenza")

Dai grafici possiamo notare l’impatto delle variabili del modello sulla predizione del peso: Le variabili settimana di gestazione, lunghezza e n.gravidanze hanno tutte un andamento simile, all’aumentare della variabile aumenta anche il peso previsto del neonato.

Per qunato riguarda la distribuzione dei residui possiamo notare una distribuzione intorno allo zero con una tendenza ad avere più residui positivi. Questo potrebbe indicare che il modello tenda leggermenta a sovrastimare i valori rispetto alla previsione perfetta.