1. CARICAMENTO E PREPARAZIONE DEI DATI

# -----------------------------------------------------------------------------
# LIBRERIE
# -----------------------------------------------------------------------------
library(tidyverse)  
library(car)         
library(MASS)        
library(broom)       
library(lmtest)      
library(caret)       
library(GGally)      

# -----------------------------------------------------------------------------
# 1. CARICAMENTO E PREPARAZIONE DEI DATI
# -----------------------------------------------------------------------------
df <- read.csv("neonati.csv")

# Conversione variabili qualitative in factor
df$Fumatrici  <- factor(df$Fumatrici,  levels = c(0, 1),
                        labels = c("Non Fumatrice", "Fumatrice"))
df$Tipo.parto <- as.factor(df$Tipo.parto)
df$Ospedale   <- as.factor(df$Ospedale)
df$Sesso      <- as.factor(df$Sesso)

# Ispezione dataset
str(df)
## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : Factor w/ 2 levels "Non Fumatrice",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Ospedale    : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
##  $ Sesso       : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...
summary(df)
##    Anni.madre     N.gravidanze             Fumatrici      Gestazione   
##  Min.   : 0.00   Min.   : 0.0000   Non Fumatrice:2396   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   Fumatrice    : 104   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000                        Median :39.00  
##  Mean   :28.16   Mean   : 0.9812                        Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000                        3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000                        Max.   :43.00  
##       Peso        Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   : 830   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :3300   Median :500.0   Median :340              osp3:835           
##  Mean   :3284   Mean   :494.7   Mean   :340                                 
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :4930   Max.   :565.0   Max.   :390

Considerazioni:

  1. La riga 1380 ha l’osservazione corrispondente ad Anni.madre pari a 0 (chiaramente risulta essere un errore).
  2. La riga 1152 ha l’osservazione relativa agli anni pari a 1.
  3. Considerato che si tratta di due sole osservazioni, posso anche considerare l’opzione di rimuoverle.
df <- df[df$Anni.madre > 1,]

2. ANALISI ESPLORATIVA (EDA)

# --- 2.1 Distribuzione del peso e ricerca outlier ---
ggplot(df, aes(y = Peso)) +
  geom_boxplot(fill = "lightblue") +
  theme_minimal() +
  labs(title = "Distribuzione del Peso Neonatale", y = "Peso (grammi)")

Possiamo notare che ci sono più outlier nella parte sinistra della distribuzione che nella parte destra.

# --- 2.2 Matrice di correlazione tra variabili numeriche ---
# Permette di anticipare problemi di multicollinearità
df |>
  dplyr::select(Peso, Lunghezza, Cranio, Gestazione, Anni.madre, N.gravidanze) |>
  ggpairs(
    title = "Matrice di Correlazione - Variabili Numeriche",
    upper = list(continuous = wrap("cor", size = 3)),
    lower = list(continuous = wrap("points", alpha = 0.2, size = 0.5)),
    diag  = list(continuous = wrap("densityDiag", fill = "lightblue"))
  )

Nella matrice di correlazione possiamo notare che: Il peso, ovvero la nostra variabile target, ha una correlazione elevata con la lunghezza e il diametro del cranio del neonato e la gestazione. Un’altra correlazione interessante è quella tra la gestazione e la lunghezza del cranio. Da tenerne conto per un’eventuale multicollinearità.

# --- 2.3 Distribuzione del peso per variabili categoriali ---
ggplot(df, aes(x = Sesso, y = Peso, fill = Sesso)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Peso per Sesso")

ggplot(df, aes(x = Fumatrici, y = Peso, fill = Fumatrici)) +
  geom_boxplot() +
  scale_fill_manual(values = c("Non Fumatrice" = "#2c3e50", "Fumatrice" = "#e74c3c")) +
  theme_minimal() +
  labs(title = "Peso per Abitudine al Fumo")

ggplot(df, aes(x = Ospedale, y = Peso, fill = Ospedale)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Peso per Ospedale")

3. TEST DI IPOTESI

# --- 3.1 Chi-quadro: parti cesarei diversi tra ospedali? ---
tabella_parti <- table(df$Ospedale, df$Tipo.parto)
print(tabella_parti)
##       
##        Ces Nat
##   osp1 242 574
##   osp2 254 594
##   osp3 232 602
# verifica delle frequenze attese affinché siano >= 5
# (assunzione del test chi-quadro)
cat("\n--- Frequenze attese (devono essere tutte >= 5) ---\n")
## 
## --- Frequenze attese (devono essere tutte >= 5) ---
print(chisq.test(tabella_parti)$expected)
##       
##             Ces      Nat
##   osp1 237.8094 578.1906
##   osp2 247.1353 600.8647
##   osp3 243.0552 590.9448
# Se tutte >= 5, possiamo usare il chi-quadro; altrimenti usiamo il test di Fisher
risultato_chi <- chisq.test(tabella_parti)
print(risultato_chi)
## 
##  Pearson's Chi-squared test
## 
## data:  tabella_parti
## X-squared = 1.083, df = 2, p-value = 0.5819

COMMENTO: Abbiamo p-value > 0.05 che ci porta a non rifiutare l’ipotesi nulla H0. I tre ospedali hanno proporzioni di cesarei statisticamente equivalenti.

# Visualizzazione proporzioni per ospedale
df |>
  count(Ospedale, Tipo.parto) |>
  group_by(Ospedale) |>
  mutate(prop = n / sum(n)) |>
  ggplot(aes(x = Ospedale, y = prop, fill = Tipo.parto)) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal() +
  labs(title = "Proporzione Tipi di Parto per Ospedale",
       y = "Proporzione", fill = "Tipo di Parto")

# --- 3.2 T-test a un campione: confronto con valori di riferimento ---
# Valori di riferimento clinici letteratura scientifica: Peso = 3300g, Lunghezza = 500mm (50cm)
# H0: la media campionaria è uguale al valore di riferimento

cat("\n--- T-test: Peso (mu = 3300g) ---\n")
## 
## --- T-test: Peso (mu = 3300g) ---
test_peso <- t.test(df$Peso, mu = 3300)
print(test_peso)
## 
##  One Sample t-test
## 
## data:  df$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

COMMENTO: Non significativo. La media campionaria di 3284g è compatibile con il riferimento clinico di 3300g. Il campione è rappresentativo della popolazione di riferimento per questa misura.

cat("\n--- T-test: Lunghezza (mu = 500mm) ---\n")
## 
## --- T-test: Lunghezza (mu = 500mm) ---
test_lunghezza <- t.test(df$Lunghezza, mu = 500)
print(test_lunghezza)
## 
##  One Sample t-test
## 
## data:  df$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6628 495.7287
## sample estimates:
## mean of x 
##  494.6958

COMMENTO: Il test in questo caso è significativo, in quanto non abbiamo elementi sufficienti per accettare l’ipotesi H0. La media campionaria è 494.7mm contro il riferimento di 500mm.

# --- 3.3 T-test a due campioni: differenze antropometriche per sesso ---
# H0: le medie sono uguali tra maschi e femmine

cat("\n--- T-test: Peso per Sesso ---\n")
## 
## --- T-test: Peso per Sesso ---
print(t.test(Peso ~ Sesso, data = df))
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496
cat("\n--- T-test: Lunghezza per Sesso ---\n")
## 
## --- T-test: Lunghezza per Sesso ---
print(t.test(Lunghezza ~ Sesso, data = df))
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750
cat("\n--- T-test: Cranio per Sesso ---\n")
## 
## --- T-test: Cranio per Sesso ---
print(t.test(Cranio ~ Sesso, data = df))
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

COMMENTO: I test sono tutti significativi. I maschi sono sistematicamente più grandi in tutte le misure. La differenza di peso di 247g è clinicamente rilevante oltre che statisticamente significativa.

4. MODELLO DI REGRESSIONE

# --- 4.1 Modello completo con interazione teoricamente motivata ---
# L'interazione Fumatrici*Gestazione testa se l'effetto delle settimane
# di gestazione sul peso cambia in funzione dello status di fumo
modello_full <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici * Gestazione +
                     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso,
                   data = df)
summary(modello_full)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici * Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1122.52  -181.96   -14.95   162.19  2611.87 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -6753.4971   142.3219 -47.452  < 2e-16 ***
## Anni.madre                        0.7781     1.1468   0.678   0.4976    
## N.gravidanze                     11.4666     4.6689   2.456   0.0141 *  
## FumatriciFumatrice              830.2103   756.1258   1.098   0.2723    
## Gestazione                       33.1481     3.8533   8.603  < 2e-16 ***
## Lunghezza                        10.2864     0.3009  34.185  < 2e-16 ***
## Cranio                           10.4670     0.4263  24.555  < 2e-16 ***
## Tipo.partoNat                    29.8321    12.0910   2.467   0.0137 *  
## Ospedaleosp2                    -10.7584    13.4494  -0.800   0.4238    
## Ospedaleosp3                     28.5942    13.5080   2.117   0.0344 *  
## SessoM                           78.1833    11.1987   6.981 3.73e-12 ***
## FumatriciFumatrice:Gestazione   -21.9184    19.2473  -1.139   0.2549    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2486 degrees of freedom
## Multiple R-squared:  0.729,  Adjusted R-squared:  0.7278 
## F-statistic: 608.1 on 11 and 2486 DF,  p-value: < 2.2e-16
# --- 4.2 Test esplicito dell'interazione con F-test ---
# Confronto diretto tra modello con e senza interazione
modello_no_interazione <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
                               Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso,
                             data = df)
cat("\n--- F-test: il termine di interazione Fumatrici:Gestazione è significativo? ---\n")
## 
## --- F-test: il termine di interazione Fumatrici:Gestazione è significativo? ---
print(anova(modello_no_interazione, modello_full))
## Analysis of Variance Table
## 
## Model 1: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici * Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2487 186743194                           
## 2   2486 186645831  1     97363 1.2968 0.2549

COMMENTO: Effettuando un F-test considerando il modello senza interazione e quello con l’interazione, possiamo giungere alla conclusione che l’interazione risulta essere non significativa.

# --- 4.3 Controllo multicollinearità (VIF) ---
# VIF > 5: potenziale problema; VIF > 10: problema serio
cat("\n--- Variance Inflation Factor (VIF) ---\n")
## 
## --- Variance Inflation Factor (VIF) ---
print(vif(modello_no_interazione))
##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.190241  1        1.090982
## N.gravidanze 1.189278  1        1.090540
## Fumatrici    1.007426  1        1.003706
## Gestazione   1.695675  1        1.302181
## Lunghezza    2.086879  1        1.444604
## Cranio       1.631049  1        1.277125
## Tipo.parto   1.004227  1        1.002111
## Ospedale     1.004267  2        1.001065
## Sesso        1.040743  1        1.020168

COMMENTO: Non è presente un problema di multicollinearità.

# --- 4.4 Selezione del modello ottimale con stepAIC ---
modello_aic <- stepAIC(modello_full, direction = "both", trace = FALSE)
cat("\n--- Modello ottimale (AIC) ---\n")
## 
## --- Modello ottimale (AIC) ---
summary(modello_aic)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.07  -181.71   -16.66   161.08  2619.57 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.9252   136.0257 -49.314  < 2e-16 ***
## N.gravidanze     12.3360     4.3344   2.846  0.00446 ** 
## Gestazione       32.0386     3.7925   8.448  < 2e-16 ***
## Lunghezza        10.3059     0.3006  34.286  < 2e-16 ***
## Cranio           10.4920     0.4257  24.648  < 2e-16 ***
## Tipo.partoNat    29.4080    12.0875   2.433  0.01505 *  
## Ospedaleosp2    -10.8939    13.4447  -0.810  0.41786    
## Ospedaleosp3     28.7917    13.4969   2.133  0.03301 *  
## SessoM           77.4657    11.1842   6.926 5.48e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2489 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 835.7 on 8 and 2489 DF,  p-value: < 2.2e-16
# Confronto anche con selezione BIC
modello_bic <- stepAIC(modello_full, direction = "both",
                       k = log(nrow(df)), trace = FALSE)
cat("\n--- Modello ottimale (BIC) ---\n")
## 
## --- Modello ottimale (BIC) ---
summary(modello_bic)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16
# Confronto AIC e BIC tra i due modelli selezionati
cat("\n--- Confronto AIC/BIC tra modello AIC-ottimale e BIC-ottimale ---\n")
## 
## --- Confronto AIC/BIC tra modello AIC-ottimale e BIC-ottimale ---
AIC(modello_aic, modello_bic)
##             df      AIC
## modello_aic 10 35143.29
## modello_bic  7 35152.89
BIC(modello_aic, modello_bic)
##             df      BIC
## modello_aic 10 35201.52
## modello_bic  7 35193.65

COMMENTO: I due modelli concordano sul nucleo (le 4 variabili biologiche + sesso) ma divergono sui predittori ospedalieri. La differenza di R^2 risulta essere minima. Inoltre risulta essere interessante il fatto che la variabile Fumatrici è stata eliminata, considerata l’assunzione teorica.

# Scelta del modello finale (di default usiamo quello selezionato da AIC)
modello_ottimale <- modello_aic

5. DIAGNOSTICA DEL MODELLO

# --- 5.1 RMSE in-sample 
rmse_insample <- sqrt(mean(modello_ottimale$residuals^2))
cat("\nRMSE in-sample:", round(rmse_insample, 2), "grammi\n")
## 
## RMSE in-sample: 273.51 grammi
# --- 5.2 RMSE con 10-fold cross-validation 
set.seed(123)
ctrl <- trainControl(method = "cv", number = 10)
modello_cv <- train(formula(modello_ottimale),
                    data    = df,
                    method  = "lm",
                    trControl = ctrl)
cat("\nRMSE 10-fold CV:", round(modello_cv$results$RMSE, 2), "grammi\n")
## 
## RMSE 10-fold CV: 274.24 grammi
cat("R² 10-fold CV:  ", round(modello_cv$results$Rsquared, 4), "\n")
## R² 10-fold CV:   0.7266

COMMENTO: Considerati i risultati ottenuti dal modello in-sample e utilizzando la cross validation per testarlo out-of-sample, possiamo dire che la differenza minima è trascurabile e il modello non fa overfitting.

# --- 5.3 Test formali sulle assunzioni della regressione ---

# Normalità dei residui (Shapiro-Wilk; affidabile per n <= 5000)
cat("\n--- Test di Shapiro-Wilk (normalità dei residui) ---\n")
## 
## --- Test di Shapiro-Wilk (normalità dei residui) ---
print(shapiro.test(residuals(modello_ottimale)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modello_ottimale)
## W = 0.97414, p-value < 2.2e-16
# Omoschedasticità (Breusch-Pagan)
# H0: varianza dei residui costante (omoschedasticità)
cat("\n--- Test di Breusch-Pagan (omoschedasticità) ---\n")
## 
## --- Test di Breusch-Pagan (omoschedasticità) ---
print(bptest(modello_ottimale))
## 
##  studentized Breusch-Pagan test
## 
## data:  modello_ottimale
## BP = 91.829, df = 8, p-value < 2.2e-16
# --- 5.4 Grafici diagnostici standard ---
par(mfrow = c(2, 2))
plot(modello_ottimale)

par(mfrow = c(1, 1))

COMMENTO NORMALITÀ: Considerando il p-value del test BP andiamo a “rifiutare” l’ipotesi H0 di normalità. Dal grafico Q-Q dei residui possiamo notare una leggera coda destra. Possiamo concludere che si tratta di una deviazione tollerabile.

COMMENTO OMOSCHEDASTICITÀ: Considerando il p-value del test BP andiamo a “rifiutare” l’ipotesi H0 di omoschedasticità. Dal grafico della Scale-Location emerge un aumento della variabilità all’aumentare del peso dei neonati.

# --- 5.5 Valori influenti: distanza di Cook ---
soglia_cook <- 4 / nrow(df)
plot(cooks.distance(modello_ottimale), type = "h",
     main = "Distanza di Cook (Valori Influenti)",
     ylab = "Distanza di Cook", xlab = "Osservazione")
abline(h = soglia_cook, col = "red", lty = 2)
legend("topright", legend = paste("Soglia =", round(soglia_cook, 5)),
       col = "red", lty = 2)

# Ispezione dei valori influenti
influenti_idx <- which(cooks.distance(modello_ottimale) > soglia_cook)
cat("\nNumero di osservazioni influenti (Cook > soglia):", length(influenti_idx), "\n")
## 
## Numero di osservazioni influenti (Cook > soglia): 116
cat("--- Dettaglio osservazioni influenti ---\n")
## --- Dettaglio osservazioni influenti ---
print(df[influenti_idx, ])
##      Anni.madre N.gravidanze     Fumatrici Gestazione Peso Lunghezza Cranio
## 13           36            5 Non Fumatrice         38 3060       455    325
## 34           27            0 Non Fumatrice         39 3150       480    382
## 119          31            0 Non Fumatrice         40 3410       550    372
## 130          30            2 Non Fumatrice         39 4240       485    352
## 134          38            6 Non Fumatrice         37 3950       500    350
## 140          29            1     Fumatrice         41 4420       530    362
## 146          24            1 Non Fumatrice         40 3820       500    320
## 155          30            0 Non Fumatrice         36 3610       410    330
## 161          35            9 Non Fumatrice         42 3760       540    348
## 220          23            1 Non Fumatrice         40 3520       445    363
## 274          28            0 Non Fumatrice         39 4320       530    350
## 295          18            0 Non Fumatrice         40 1850       460    305
## 310          40            3 Non Fumatrice         28 1560       420    379
## 329          25            1 Non Fumatrice         40 4560       540    340
## 375          33            4 Non Fumatrice         38 4270       510    353
## 377          31            1 Non Fumatrice         39 3400       460    325
## 390          38            0 Non Fumatrice         40 3700       470    320
## 455          24            0 Non Fumatrice         38 2920       505    351
## 472          25            2 Non Fumatrice         41 3990       495    335
## 478          32            3 Non Fumatrice         42 4310       505    365
## 516          40            8 Non Fumatrice         38 3520       470    341
## 582          30            7 Non Fumatrice         35 2220       470    316
## 615          27            1 Non Fumatrice         35 2100       440    345
## 616          29            0 Non Fumatrice         42 3540       540    368
## 632          21            0 Non Fumatrice         37 2750       510    333
## 633          30            0 Non Fumatrice         40 3860       480    338
## 648          35            1 Non Fumatrice         38 2850       500    360
## 656          38            3 Non Fumatrice         41 2320       450    330
## 657          32            3 Non Fumatrice         41 2970       520    318
## 684          30            1 Non Fumatrice         39 3000       475    390
## 791          30            1 Non Fumatrice         41 4440       510    335
## 805          30            2 Non Fumatrice         29 1190       360    272
## 828          32            6 Non Fumatrice         40 4200       510    350
## 890          33            2 Non Fumatrice         39 3660       490    318
## 908          35            2 Non Fumatrice         39 3300       460    320
## 928          25            0 Non Fumatrice         28  830       310    254
## 950          25            0 Non Fumatrice         41 2980       500    356
## 1008         32            5 Non Fumatrice         36 2400       470    333
## 1014         17            0 Non Fumatrice         37 2050       390    295
## 1036         26            0 Non Fumatrice         40 4330       500    355
## 1130         33           11 Non Fumatrice         43 3400       475    360
## 1132         21            0 Non Fumatrice         39 3950       500    330
## 1137         26            0 Non Fumatrice         40 3000       520    340
## 1181         30            1 Non Fumatrice         36 4070       500    373
## 1188         21            0 Non Fumatrice         40 4140       550    320
## 1192         24            2 Non Fumatrice         40 2560       490    323
## 1194         41            1 Non Fumatrice         39 2160       450    325
## 1212         26            0 Non Fumatrice         34 2050       460    325
## 1230         29            1 Non Fumatrice         41 4010       470    353
## 1253         35            4 Non Fumatrice         39 3450       523    365
## 1267         24            0 Non Fumatrice         39 4260       545    343
## 1268         28            1 Non Fumatrice         40 3790       460    332
## 1273         32            1 Non Fumatrice         33 2040       480    307
## 1287         35            3 Non Fumatrice         40 2660       485    335
## 1293         30            3 Non Fumatrice         38 4600       485    380
## 1306         23            0 Non Fumatrice         41 4900       510    352
## 1317         28            2 Non Fumatrice         39 3650       500    315
## 1341         32            3 Non Fumatrice         36 3780       480    349
## 1395         30            2 Non Fumatrice         39 3790       505    304
## 1399         42            2 Non Fumatrice         38 2560       525    349
## 1402         35            1 Non Fumatrice         39 2660       460    364
## 1426         32            4     Fumatrice         39 2250       460    319
## 1429         24            4 Non Fumatrice         29 1280       390    355
## 1433         31            4 Non Fumatrice         41 4810       530    364
## 1450         36            8 Non Fumatrice         41 3730       480    335
## 1472         37            4 Non Fumatrice         37 3720       480    355
## 1499         30            1 Non Fumatrice         39 2380       495    325
## 1505         30            8 Non Fumatrice         39 2860       490    337
## 1541         30            0 Non Fumatrice         38 4540       530    343
## 1551         35            1 Non Fumatrice         38 4370       315    374
## 1553         30            4 Non Fumatrice         35 4520       520    360
## 1556         37            0 Non Fumatrice         41 2420       490    300
## 1588         22            0 Non Fumatrice         42 3100       510    361
## 1593         41            3 Non Fumatrice         35 1500       420    304
## 1619         31            0 Non Fumatrice         31  990       340    278
## 1635         32            2 Non Fumatrice         39 3430       445    322
## 1639         39            3 Non Fumatrice         40 4760       550    365
## 1694         23            1 Non Fumatrice         36 3850       460    334
## 1702         21            0 Non Fumatrice         41 3000       495    357
## 1712         28            0 Non Fumatrice         39 3800       520    300
## 1718         34            4 Non Fumatrice         42 2660       500    320
## 1780         25            2 Non Fumatrice         25  900       325    253
## 1838         34            3 Non Fumatrice         40 4580       515    360
## 1856         33            0 Non Fumatrice         36 3110       465    300
## 1868         29            0 Non Fumatrice         40 3470       525    390
## 1893         22            1 Non Fumatrice         34 3030       470    312
## 1915         27            1 Non Fumatrice         38 2480       500    315
## 1920         26            0     Fumatrice         39 4930       550    350
## 1926         29            0 Non Fumatrice         41 3840       485    358
## 1937         36            1 Non Fumatrice         41 2940       500    358
## 1962         38            4 Non Fumatrice         38 4370       530    340
## 1963         27            0 Non Fumatrice         42 4700       540    362
## 2023         27            1 Non Fumatrice         39 4650       510    354
## 2040         27            1 Non Fumatrice         38 3240       410    359
## 2076         23            0 Non Fumatrice         40 4720       540    360
## 2086         26            8 Non Fumatrice         40 3250       500    355
## 2114         36            0 Non Fumatrice         31 1180       355    270
## 2115         35            1 Non Fumatrice         32 1890       500    309
## 2120         32            0 Non Fumatrice         27 1140       370    267
## 2123         32            0 Non Fumatrice         39 3950       520    330
## 2135         33            3 Non Fumatrice         39 3950       480    345
## 2175         37            8 Non Fumatrice         28  930       355    235
## 2195         38            1 Non Fumatrice         40 3980       480    335
## 2204         33            0 Non Fumatrice         39 4000       490    360
## 2219         37            1 Non Fumatrice         39 2500       490    352
## 2221         35           10 Non Fumatrice         39 2950       495    335
## 2225         27            0 Non Fumatrice         35 3140       465    290
## 2287         33            3 Non Fumatrice         40 3640       470    337
## 2315         24            0 Non Fumatrice         42 2800       520    340
## 2317         25            6 Non Fumatrice         38 2530       460    340
## 2392         28            1 Non Fumatrice         40 4720       540    355
## 2421         20            0 Non Fumatrice         41 4160       500    358
## 2422         33           10 Non Fumatrice         40 3090       485    353
## 2437         28            1 Non Fumatrice         27  980       320    265
## 2452         28            0 Non Fumatrice         26  930       345    245
## 2471         34           10 Non Fumatrice         38 2880       470    345
##      Tipo.parto Ospedale Sesso
## 13          Ces     osp1     F
## 34          Nat     osp1     F
## 119         Nat     osp2     M
## 130         Nat     osp2     M
## 134         Nat     osp3     M
## 140         Ces     osp2     F
## 146         Ces     osp2     F
## 155         Nat     osp1     M
## 161         Nat     osp2     F
## 220         Nat     osp1     F
## 274         Ces     osp3     M
## 295         Nat     osp3     F
## 310         Nat     osp3     F
## 329         Nat     osp1     M
## 375         Nat     osp2     M
## 377         Nat     osp1     M
## 390         Nat     osp1     M
## 455         Ces     osp3     M
## 472         Nat     osp3     M
## 478         Nat     osp1     M
## 516         Nat     osp3     M
## 582         Nat     osp3     M
## 615         Ces     osp2     F
## 616         Nat     osp1     M
## 632         Nat     osp1     M
## 633         Nat     osp1     M
## 648         Nat     osp1     F
## 656         Nat     osp2     F
## 657         Ces     osp2     M
## 684         Ces     osp2     F
## 791         Nat     osp3     M
## 805         Nat     osp2     F
## 828         Nat     osp1     M
## 890         Nat     osp2     F
## 908         Ces     osp1     M
## 928         Nat     osp1     F
## 950         Nat     osp3     M
## 1008        Nat     osp1     F
## 1014        Nat     osp2     F
## 1036        Nat     osp3     F
## 1130        Nat     osp1     M
## 1132        Nat     osp1     F
## 1137        Ces     osp3     M
## 1181        Nat     osp2     M
## 1188        Ces     osp1     M
## 1192        Nat     osp1     M
## 1194        Ces     osp3     M
## 1212        Nat     osp1     F
## 1230        Nat     osp2     M
## 1253        Nat     osp1     M
## 1267        Ces     osp2     F
## 1268        Nat     osp2     F
## 1273        Ces     osp1     F
## 1287        Ces     osp3     M
## 1293        Nat     osp1     M
## 1306        Nat     osp2     F
## 1317        Nat     osp1     M
## 1341        Nat     osp3     M
## 1395        Ces     osp3     M
## 1399        Ces     osp2     M
## 1402        Ces     osp1     F
## 1426        Nat     osp2     F
## 1429        Nat     osp1     F
## 1433        Ces     osp3     M
## 1450        Nat     osp3     M
## 1472        Ces     osp1     F
## 1499        Nat     osp1     F
## 1505        Ces     osp2     F
## 1541        Ces     osp3     M
## 1551        Nat     osp3     F
## 1553        Nat     osp2     F
## 1556        Ces     osp1     M
## 1588        Nat     osp2     F
## 1593        Nat     osp1     M
## 1619        Ces     osp2     F
## 1635        Ces     osp1     F
## 1639        Nat     osp2     F
## 1694        Ces     osp3     F
## 1702        Ces     osp1     M
## 1712        Nat     osp3     F
## 1718        Nat     osp2     F
## 1780        Nat     osp3     F
## 1838        Nat     osp2     M
## 1856        Ces     osp3     M
## 1868        Nat     osp1     M
## 1893        Nat     osp2     F
## 1915        Ces     osp2     M
## 1920        Ces     osp2     F
## 1926        Ces     osp2     F
## 1937        Nat     osp3     M
## 1962        Nat     osp3     F
## 1963        Nat     osp3     F
## 2023        Nat     osp2     M
## 2040        Ces     osp1     F
## 2076        Nat     osp3     F
## 2086        Nat     osp2     M
## 2114        Nat     osp3     F
## 2115        Nat     osp2     F
## 2120        Nat     osp3     F
## 2123        Ces     osp1     F
## 2135        Nat     osp3     M
## 2175        Nat     osp1     F
## 2195        Nat     osp1     F
## 2204        Ces     osp3     F
## 2219        Ces     osp1     M
## 2221        Nat     osp1     F
## 2225        Nat     osp2     F
## 2287        Nat     osp3     F
## 2315        Ces     osp2     M
## 2317        Nat     osp1     F
## 2392        Nat     osp3     M
## 2421        Ces     osp1     M
## 2422        Nat     osp3     M
## 2437        Nat     osp1     M
## 2452        Ces     osp3     F
## 2471        Ces     osp2     M

COMMENTO: 116 osservazioni (circa il 4.6% del campione) superano la soglia di Cook. Tra le prime si trovano casi come l’osservazione 140 (madre fumatrice, neonato di 4420 g a 41 settimane) e la 130 (4240 g con lunghezza pari a 485 mm), che presentano combinazioni di valori piuttosto insolite. Si tratta di punti che il modello riesce a rappresentare con più difficoltà e che possono avere un certo peso nella stima complessiva.

6. VISUALIZZAZIONI PRINCIPALI

# --- 6.1 Impatto gestazione e fumo sul peso ---
ggplot(df, aes(x = Gestazione, y = Peso, color = Fumatrici)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() +
  scale_color_manual(values = c("Non Fumatrice" = "#2c3e50",
                                "Fumatrice"     = "#e74c3c")) +
  labs(
    title    = "Impatto della Gestazione e del Fumo sul Peso Neonatale",
    subtitle = "Pendenze diverse indicano un effetto di interazione",
    x        = "Settimane di Gestazione",
    y        = "Peso del Neonato (g)",
    color    = "Status Madre"
  )

COMMENTO: A parità di settimane di gestazione, i neonati di madri fumatrici tendono a pesare meno. In entrambi i gruppi all’aumentare della gestazione -> aumenta il peso.

# --- 6.3 Coefficienti del modello con intervalli di confidenza ---
tidy(modello_ottimale, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  ggplot(aes(x = reorder(term, estimate), y = estimate,
             ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  coord_flip() +
  theme_minimal() +
  labs(title = "Coefficienti del Modello con IC al 95%",
       x = NULL, y = "Stima (grammi)")

7. PREVISIONE

# Neonata femmina, madre alla 3a gravidanza, 39 settimane di gestazione.

nuovo_caso <- data.frame(
  Anni.madre = 30,
  N.gravidanze = 3,
  Fumatrici = factor("Non Fumatrice",
                     levels = levels(df$Fumatrici)),
  Gestazione = 39,
  Lunghezza  = round(mean(df$Lunghezza)),    # media campionaria
  Cranio     = round(mean(df$Cranio)),        # media campionaria
  Tipo.parto = factor("Nat",
                      levels = levels(df$Tipo.parto)),
  Ospedale   = factor("osp1",
                      levels = levels(df$Ospedale)),
  Sesso      = factor("F",
                      levels = levels(df$Sesso))
)

previsione <- predict(modello_ottimale,
                      newdata  = nuovo_caso,
                      interval = "prediction",
                      level    = 0.95)

cat("\n=== PREVISIONE PESO NEONATO ===\n")
## 
## === PREVISIONE PESO NEONATO ===
cat("Stima puntuale: ", round(previsione[1], 0), "grammi\n")
## Stima puntuale:  3277 grammi
cat("Intervallo di previsione al 95%: [",
    round(previsione[2], 0), ",",
    round(previsione[3], 0), "] grammi\n")
## Intervallo di previsione al 95%: [ 2739 , 3815 ] grammi
cat("(L'intervallo di PREVISIONE, più ampio dell'IC sulla media,",
    "riflette l'incertezza sulla singola nascita futura)\n")
## (L'intervallo di PREVISIONE, più ampio dell'IC sulla media, riflette l'incertezza sulla singola nascita futura)