neonati

knitr::opts_chunk$set(echo = TRUE)
# sfondo bianco per tutte le figure prodotte da knitr
knitr::opts_chunk$set(dev.args = list(bg = "white"))

# colori di testo/assi/titoli ben visibili
par(col = "black", col.axis = "black", col.lab = "black", col.main = "black")

# Librerie 
library(janitor)    # clean_names

## Warning: package 'janitor' was built under R version 4.3.3

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(dplyr)      # tabelle sintetiche

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(knitr)      # kable

## Warning: package 'knitr' was built under R version 4.3.3

library(moments)    # skewness/kurtosis
library(lmtest)     # bptest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)        # outlierTest, vif

## Warning: package 'car' was built under R version 4.3.3

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

df <- read.csv("/Users/chiaratombolini/Desktop/RSTUDIO/neonati.csv", stringsAsFactors = FALSE)
df <- janitor::clean_names(df)

# Variabili categoriche come fattori
df$tipo_parto <- as.factor(df$tipo_parto)
df$ospedale   <- as.factor(df$ospedale)
df$sesso      <- as.factor(df$sesso)
df$fumatrici  <- as.factor(df$fumatrici)

# Struttura + statistiche sintetiche 
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 "0","1": 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         peso     
##  Min.   : 0.00   Min.   : 0.0000   0:2396    Min.   :25.00   Min.   : 830  
##  1st Qu.:25.00   1st Qu.: 0.0000   1: 104    1st Qu.:38.00   1st Qu.:2990  
##  Median :28.00   Median : 1.0000             Median :39.00   Median :3300  
##  Mean   :28.16   Mean   : 0.9812             Mean   :38.98   Mean   :3284  
##  3rd Qu.:32.00   3rd Qu.: 1.0000             3rd Qu.:40.00   3rd Qu.:3620  
##  Max.   :46.00   Max.   :12.0000             Max.   :43.00   Max.   :4930  
##    lunghezza         cranio    tipo_parto ospedale   sesso   
##  Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :500.0   Median :340              osp3:835           
##  Mean   :494.7   Mean   :340                                 
##  3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :565.0   Max.   :390

# Distribuzioni (unità coerenti: peso=g, lunghezza=mm, età=anni, gestazione=settimane)
par(mfrow=c(2,2))
hist(df$peso,       main="Distribuzione Peso",       xlab="Peso (g)")
hist(df$lunghezza,  main="Distribuzione Lunghezza",  xlab="Lunghezza (mm)")
hist(df$anni_madre, main="Distribuzione Età madre",  xlab="Età (anni)")
hist(df$gestazione, main="Distribuzione Settimane",  xlab="Settimane di gestazione")

par(mfrow=c(1,1))
# Outlier per peso e lunghezza 
outlier_iqr <- function(x) {
  x <- x[!is.na(x)]
  q1 <- quantile(x, 0.25); q3 <- quantile(x, 0.75)
  iqr <- q3 - q1
  low <- q1 - 1.5 * iqr; high <- q3 + 1.5 * iqr
  list(bounds=c(low=low, high=high))
}
oi_peso      <- outlier_iqr(df$peso)
oi_lunghezza <- outlier_iqr(df$lunghezza)

# Boxplot
boxplot(df$peso,      main="Boxplot Peso",      ylab="Peso (g)")

boxplot(df$lunghezza, main="Boxplot Lunghezza", ylab="Lunghezza (mm)")

# Soglie IQR stampate in tabella
kable(data.frame(
  variabile=c("peso","lunghezza"),
  low=c(oi_peso$bounds["low"], oi_lunghezza$bounds["low"]),
  high=c(oi_peso$bounds["high"], oi_lunghezza$bounds["high"])
), digits=1, caption="Soglie IQR (flag outlier potenziali)")

Soglie IQR (flag outlier potenziali)

variabile	low	high
peso	NA	NA
lunghezza	NA	NA

#(a) Tipo di parto ~ Ospedale (test di indipendenza)
tab_parti <- table(df$ospedale, df$tipo_parto)
tab_parti

##       
##        Ces Nat
##   osp1 242 574
##   osp2 254 595
##   osp3 232 603

chi_res <- chisq.test(tab_parti)
chi_res

## 
##  Pearson's Chi-squared test
## 
## data:  tab_parti
## X-squared = 1.0972, df = 2, p-value = 0.5778

#commento 
if (chi_res$p.value > 0.05) {
  cat("Esito: p =", signif(chi_res$p.value,3), "→ non si rifiuta H0: nel campione non c'è evidenza di **associazione** tra tipo di parto e ospedale.\n")
} else {
  cat("Esito: p =", signif(chi_res$p.value,3), "→ si rifiuta H0: nel campione c'è evidenza di **associazione** tra tipo di parto e ospedale.\n")
}

## Esito: p = 0.578 → non si rifiuta H0: nel campione non c'è evidenza di **associazione** tra tipo di parto e ospedale.

#(b) Media del peso e della lunghezza vs popolazione
mu_peso       <- 3300   # g
mu_lunghezza  <- 500    # mm

t_peso_mu <- t.test(df$peso,      mu = mu_peso)
t_lung_mu <- t.test(df$lunghezza, mu = mu_lunghezza)

t_peso_mu

## 
##  One Sample t-test
## 
## data:  df$peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

t_lung_mu

## 
##  One Sample t-test
## 
## data:  df$lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

Commento alle distribuzioni e ai boxplot: Le distribuzioni delle variabili quantitative mostrano un comportamento complessivamente coerente con quanto atteso in un campione neonatale. Peso: la distribuzione risulta leggermente asimmetrica verso sinistra, con la maggior parte dei valori concentrati tra 2800 e 3600 grammi. Il boxplot evidenzia pochi valori inferiori considerabili outlier, riconducibili a neonati pretermine o di basso peso alla nascita. Lunghezza: presenta un andamento simile, con una forma quasi normale ma una leggera coda sinistra, coerente con la presenza di alcuni soggetti più piccoli rispetto alla media. Anche qui gli outlier inferiori coincidono con i pesi più bassi, suggerendo una coerenza interna del dataset. Età della madre: l’istogramma mostra una distribuzione concentrata tra 25 e 35 anni, con alcuni casi isolati oltre i 40. Il boxplot rivela la presenza di pochi valori anomali elevati, che potrebbero essere dovuti a errori di digitazione o a età materne particolarmente alte; questi casi andrebbero eventualmente verificati o imputati. Settimane di gestazione: la distribuzione è fortemente concentrata tra 37 e 41 settimane, come atteso per nascite a termine, con una coda verso valori inferiori (prematurità). Gli outlier identificati rappresentano neonati nati molto prima del termine fisiologico.

Fonti peso e lunghezza: Peso medio neonati a termine: ≈ 3300 g (Fonte: Istituto Superiore di Sanità, Rapporto 2023). Lunghezza media neonati a termine: ≈ 500 mm (50 cm) (Fonte: Ministero della Salute, “Crescita e sviluppo del neonato”, 2023).

-Commento alla media del peso e della lunghezza vs popolazione:

cat("Peso medio campione:", round(mean(df$peso, na.rm=TRUE),1), "g; riferimento:", mu_peso, "g.\n")

## Peso medio campione: 3284.1 g; riferimento: 3300 g.

if (t_peso_mu$p.value > 0.05) {
  cat("→ Differenza non statisticamente significativa.\n")
} else {
  cat("→ Differenza statisticamente significativa (p =", signif(t_peso_mu$p.value,3), ").\n")
}

## → Differenza non statisticamente significativa.

cat("Lunghezza media campione:", round(mean(df$lunghezza, na.rm=TRUE),1), "mm; riferimento:", mu_lunghezza, "mm.\n")

## Lunghezza media campione: 494.7 mm; riferimento: 500 mm.

if (t_lung_mu$p.value > 0.05) {
  cat("→ Differenza non statisticamente significativa.\n")
} else {
  cat("→ Differenza statisticamente significativa (p =", signif(t_lung_mu$p.value,3), ").\n")
}

## → Differenza statisticamente significativa (p = 1.81e-23 ).

-Misure antropometriche diverse tra i sessi:

t_peso_sex <- t.test(peso      ~ sesso, data = df)   # peso M vs F
t_lung_sex <- t.test(lunghezza ~ sesso, data = df)   # lunghezza M vs F
t_cran_sex <- t.test(cranio    ~ sesso, data = df)   # cranio M vs F

t_peso_sex; t_lung_sex; t_cran_sex

## 
##  Welch Two Sample t-test
## 
## data:  peso by sesso
## t = -12.106, df = 2490.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.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

## 
##  Welch Two Sample t-test
## 
## data:  lunghezza by sesso
## t = -9.582, df = 2459.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.929470  -7.876273
## sample estimates:
## mean in group F mean in group M 
##        489.7643        499.6672

## 
##  Welch Two Sample t-test
## 
## data:  cranio by sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M 
##        337.6330        342.4486

Commento:

msg <- function(tt, nome){
  if (tt$p.value < 0.05) paste0(nome, ": differenza **significativa** tra i sessi (p = ", signif(tt$p.value,3), ").\n")
  else paste0(nome, ": differenza **non significativa** (p = ", signif(tt$p.value,3), ").\n")
}
cat(msg(t_peso_sex,"Peso"))

## Peso: differenza **significativa** tra i sessi (p = 8.03e-33).

cat(msg(t_lung_sex,"Lunghezza"))

## Lunghezza: differenza **significativa** tra i sessi (p = 2.24e-21).

cat(msg(t_cran_sex,"Cranio"))

## Cranio: differenza **significativa** tra i sessi (p = 1.72e-13).

Indici di forma

sk <- moments::skewness(df$peso, na.rm = TRUE)
ku <- moments::kurtosis(df$peso, na.rm = TRUE) - 3  # curtosi in eccesso
kable(data.frame(
  skewness = round(sk,2),
  excess_kurtosis = round(ku,2)
), caption="Asimmetria e curtosi (in eccesso) del peso")

Asimmetria e curtosi (in eccesso) del peso

skewness	excess_kurtosis
-0.65	2.03

Commento:

if (sk < 0) cat("Asimmetria negativa: coda a sinistra (pesi più bassi relativamente comuni).\n") else if (sk > 0) cat("Asimmetria positiva: coda a destra.\n") else cat("Distribuzione simmetrica.\n")

## Asimmetria negativa: coda a sinistra (pesi più bassi relativamente comuni).

if (ku > 0) cat("Curtosi in eccesso > 0: code più pesanti della normale (possibili valori estremi).\n") else if (ku < 0) cat("Curtosi in eccesso < 0: code più leggere della normale.\n")

## Curtosi in eccesso > 0: code più pesanti della normale (possibili valori estremi).

Relazione bivariata: peso ~ gestazione e Correlazioni tra numeriche

plot(df$gestazione, df$peso,
     xlab = "Settimane di gestazione",
     ylab = "Peso (g)",
     main = "Relazione Peso ~ Settimane Gestazione",
     pch = 19, col = "darkblue")
abline(lm(peso ~ gestazione, data = df), col = "red")

# Correlazioni tra variabili numeriche
library(dplyr)
library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

num_df <- select(df, where(is.numeric))
cor_mat <- round(cor(num_df, use = "pairwise.complete.obs"), 2)

kable(cor_mat, caption = "Matrice di correlazione tra variabili numeriche") %>%
  kable_styling(full_width = F, position = "center", bootstrap_options = c("striped", "hover"))

Matrice di correlazione tra variabili numeriche

	anni_madre	n_gravidanze	gestazione	peso	lunghezza	cranio
anni_madre	1.00	0.38	-0.14	-0.02	-0.06	0.02
n_gravidanze	0.38	1.00	-0.10	0.00	-0.06	0.04
gestazione	-0.14	-0.10	1.00	0.59	0.62	0.46
peso	-0.02	0.00	0.59	1.00	0.80	0.70
lunghezza	-0.06	-0.06	0.62	0.80	1.00	0.60
cranio	0.02	0.04	0.46	0.70	0.60	1.00

-Il grafico mostra una relazione positiva e lineare tra il peso del neonato e le settimane di gestazione: all’aumentare della durata della gravidanza, il peso cresce in modo pressoché costante. La retta di regressione evidenzia una tendenza crescente chiara, coerente con l’andamento fisiologico dello sviluppo fetale. La dispersione dei punti attorno alla linea rossa indica che, pur essendoci una variabilità individuale, il trend generale è ben definito: i neonati nati prima della 35ª settimana hanno pesi sensibilmente inferiori rispetto a quelli nati a termine (≥ 37 settimane).

-La matrice di correlazione mostra una relazione positiva marcata tra il peso e le principali misure antropometriche (lunghezza e circonferenza cranica), confermando che i neonati più lunghi e con cranio più ampio tendono anche a pesare di più. Si osserva inoltre una correlazione positiva moderata tra il peso e le settimane di gestazione, coerente con l’aumento fisiologico del peso al progredire della gravidanza. Le variabili età della madre e numero di gravidanze risultano invece debolmente correlate con il peso neonatale, suggerendo che il loro effetto diretto è limitato o mediato da altri fattori. -Modello regressione multipla

mod_full <- lm(
  peso ~ anni_madre + n_gravidanze + fumatrici +
    gestazione + lunghezza + cranio +
    tipo_parto + ospedale + sesso,
  data = df
)
summary(mod_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 
## -1124.40  -181.66   -14.42   160.91  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6738.4762   141.3087 -47.686  < 2e-16 ***
## anni_madre        0.8921     1.1323   0.788   0.4308    
## n_gravidanze     11.2665     4.6608   2.417   0.0157 *  
## fumatrici1      -30.1631    27.5386  -1.095   0.2735    
## gestazione       32.5696     3.8187   8.529  < 2e-16 ***
## lunghezza        10.2945     0.3007  34.236  < 2e-16 ***
## cranio           10.4707     0.4260  24.578  < 2e-16 ***
## tipo_partoNat    29.5254    12.0844   2.443   0.0146 *  
## ospedaleosp2    -11.2095    13.4379  -0.834   0.4043    
## ospedaleosp3     28.0958    13.4957   2.082   0.0375 *  
## sessoM           77.5409    11.1776   6.937 5.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 669.2 on 10 and 2489 DF,  p-value: < 2.2e-16

# Stepwise 
mod_step <- step(mod_full, direction = "both", trace = FALSE)
summary(mod_step)

## 
## Call:
## lm(formula = peso ~ n_gravidanze + gestazione + lunghezza + cranio + 
##     tipo_parto + ospedale + sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## n_gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## cranio           10.4922     0.4254  24.661  < 2e-16 ***
## tipo_partoNat    29.2803    12.0817   2.424  0.01544 *  
## ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## sessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16

# Confronto informativo: AIC e BIC 
AIC(mod_full, mod_step)

##          df      AIC
## mod_full 12 35171.95
## mod_step 10 35169.79

BIC(mod_full, mod_step)

##          df      BIC
## mod_full 12 35241.84
## mod_step 10 35228.03

v_vals <- car::vif(mod_step)
# ---- Analisi di multicollinearità ----
library(car)
library(knitr)
library(kableExtra)

v <- car::vif(mod_step)

# Costruzione tabella  con GVIF e VIF corretti
as_vif_table <- function(vif_obj) {
  if (is.matrix(vif_obj)) {
    out <- data.frame(
      Predittore = rownames(vif_obj),
      Df         = vif_obj[, "Df"],
      GVIF       = round(vif_obj[, "GVIF"], 3),
      GVIF_adj   = round(vif_obj[, "GVIF^(1/(2*Df))"], 3)
    )
  } else {
    out <- data.frame(
      Predittore = names(vif_obj),
      VIF        = round(as.numeric(vif_obj), 3)
    )
  }
  rownames(out) <- NULL
  out
}

vif_table <- as_vif_table(v)

kable(vif_table, caption = "Analisi di multicollinearità: VIF e GVIF") %>%
  kable_styling(full_width = F, position = "center",
                bootstrap_options = c("striped", "hover"))

Analisi di multicollinearità: VIF e GVIF

Predittore	Df	GVIF	GVIF_adj
n_gravidanze	1	1.025	1.013
gestazione	1	1.670	1.292
lunghezza	1	2.082	1.443
cranio	1	1.627	1.275
tipo_parto	1	1.004	1.002
ospedale	2	1.003	1.001
sesso	1	1.040	1.020

# Metriche principali
R2   <- summary(mod_step)$r.squared
res  <- residuals(mod_step)
RMSE <- sqrt(mean(res^2, na.rm = TRUE))

R2; RMSE  # ~0.729; ~273 g attesi in linea con il tuo testo

## [1] 0.7286934

## [1] 273.4227

# Diagnostica standard
par(mfrow = c(2, 2))
plot(mod_step)

par(mfrow = c(1, 1))

# Test diagnostici
lmtest::bptest(mod_step)                  # eteroschedasticità

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_step
## BP = 91.768, df = 8, p-value < 2.2e-16

shapiro.test(residuals(mod_step))         # normalità residui

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_step)
## W = 0.97408, p-value < 2.2e-16

car::outlierTest(mod_step)                # outlier influenti

##       rstudent unadjusted p-value Bonferroni p
## 1551 10.004278         3.9642e-23   9.9104e-20
## 155   5.046640         4.8210e-07   1.2053e-03
## 1306  4.872419         1.1716e-06   2.9291e-03

# Distanza di Cook
par(mfrow = c(1, 1), mar=c(5,4,2,1))
plot(cooks.distance(mod_step), main = "Distanza di Cook", ylab="Cook's Distance", xlab="Osservazioni")

# Leverage
lev <- hatvalues(mod_step)
par(mar=c(5,4,2,1))
plot(lev, main = "Leverage", ylab="Hat values", xlab="Osservazioni")
abline(h = 2 * mean(lev), col = "red")

Commento modello: I quattro pannelli mostrano la diagnostica del modello lineare stimato sul peso neonatale. (1) Residuals vs Fitted: i residui si distribuiscono in modo abbastanza simmetrico attorno allo zero, ma si nota una lieve eteroschedasticità (varianza dei residui che aumenta con i valori previsti). Ciò suggerisce che l’ipotesi di omoschedasticità non è pienamente soddisfatta, sebbene la violazione sembri moderata. (2) Q-Q plot dei residui: la maggior parte dei punti segue la linea di riferimento, ma nelle code si osservano deviazioni, indicando residui non perfettamente normali. Considerando però la dimensione del campione, tale scostamento è accettabile e non compromette la bontà complessiva del modello. (3) Scale-Location plot: conferma la presenza di eteroschedasticità lieve, con un aumento della dispersione per i valori più elevati di peso previsto. In tali casi, l’uso di errori standard robusti (HC) può migliorare l’affidabilità delle inferenze. (4) Residuals vs Leverage: la maggior parte dei punti mostra valori di leverage contenuti; alcuni casi (es. osservazioni 1551, 1306) presentano valori più elevati ma non mostrano influenze eccessive sulla stima complessiva, come indicato anche dai valori di distanza di Cook inferiori alla soglia critica. In sintesi, la diagnostica evidenzia un modello globalmente adeguato, con leggere violazioni delle assunzioni classiche (normalità e omoschedasticità) ma senza problemi di influenza marcata o specificazioni errate.

-nel confronto AIC/BIC il modello mod_step presenta valori più bassi rispetto a mod_full, quindi è preferibile a parità di contesto interpretativo (miglior trade-off bontà/complessità).

-nessun VIF oltre 5: la collinearità è accettabile; le misure antropometriche mostrano relazione ma non tale da compromettere la stabilità delle stime. -R² ≈ 0,729 e RMSE ≈ 273 g: il modello spiega ~73% della variabilità con un errore medio di previsione nell’ordine di ±270–280 g.

Il grafico della Distanza di Cook mostra che la quasi totalità delle osservazioni presenta valori molto bassi, ben al di sotto della soglia comunemente usata (≈ 0.5 o 1). Ciò indica che nessuna singola osservazione esercita un’influenza eccessiva sui coefficienti stimati del modello. L’unico punto con un valore leggermente più alto non supera comunque i limiti di attenzione e non altera significativamente le stime. Nel grafico del Leverage, la linea rossa rappresenta il valore soglia di riferimento. La maggior parte dei punti si colloca ampiamente al di sotto di questa soglia, confermando che i dati non contengono casi con leva anomala (ovvero, osservazioni con combinazioni di predittori troppo distanti dalla media). Nel complesso, l’analisi diagnostica sui punti influenti conferma che il modello stimato è stabile e non dominato da osservazioni estreme o anomale, garantendo l’affidabilità delle stime ottenute.

Interpretazione dei coefficienti (a parità delle altre variabili): Gestazione: ~+32 g per ogni settimana aggiuntiva. Sesso (M vs F): ~+77 g per i maschi rispetto alle femmine (livello di riferimento). Lunghezza: ~+10,3 g per ogni mm in più. Cranio: ~+10,5 g per ogni mm in più. n_gravidanze: effetto positivo e modesto (~+12 g per gravidanza). Sintesi: il peso è spiegato soprattutto da settimane di gestazione, misure antropometriche (lunghezza, cranio) e sesso; gli altri effetti sono minori o poco stabili.

# Funzioni di supporto
moda <- function(x) { ux <- na.omit(x); if (!length(ux)) return(NA); names(sort(table(ux), decreasing = TRUE))[1] }
typ_num <- function(x) median(x, na.rm = TRUE)


newx <- data.frame(
  anni_madre   = if ("anni_madre"   %in% names(df)) typ_num(df$anni_madre)   else NA_real_,
  n_gravidanze = if ("n_gravidanze" %in% names(df)) 3                         else NA_real_,
  fumatrici    = if ("fumatrici"    %in% names(df)) factor(moda(df$fumatrici), levels = levels(df$fumatrici)) else NA,
  gestazione   = if ("gestazione"   %in% names(df)) 39                        else NA_real_,
  lunghezza    = if ("lunghezza"    %in% names(df)) typ_num(df$lunghezza)     else NA_real_,
  cranio       = if ("cranio"       %in% names(df)) typ_num(df$cranio)        else NA_real_,
  tipo_parto   = if ("tipo_parto"   %in% names(df)) factor(moda(df$tipo_parto), levels = levels(df$tipo_parto)) else NA,
  ospedale     = if ("ospedale"     %in% names(df)) factor(moda(df$ospedale), levels = levels(df$ospedale))     else NA,
  sesso        = if ("sesso"        %in% names(df)) factor("F", levels = levels(df$sesso))                      else NA
)

# Allineo alle variabili effettive del modello stepwise
nm_needed <- all.vars(formula(mod_step))[-1]
missing_cols <- setdiff(nm_needed, colnames(newx))
for (mc in missing_cols) {
  if (mc %in% names(df)) {
    newx[[mc]] <- if (is.numeric(df[[mc]])) typ_num(df[[mc]]) else factor(moda(df[[mc]]), levels = levels(df[[mc]]))
  } else {
    newx[[mc]] <- NA
  }
}
newx <- newx[, nm_needed, drop = FALSE]

# Previsioni: punto, IP (singola osservazione), IC (media attesa)
pred_point <- predict(mod_step, newdata = newx)
pred_pred  <- predict(mod_step, newdata = newx, interval = "prediction")  # singola osservazione
pred_conf  <- predict(mod_step, newdata = newx, interval = "confidence")  # media attesa

pred_point; pred_pred; pred_conf

##        1 
## 3317.222

##        fit    lwr      upr
## 1 3317.222 2779.3 3855.145

##        fit      lwr      upr
## 1 3317.222 3287.977 3346.467

Il modello stima per il profilo “tipico” un peso medio previsto di circa 3317 g. L’intervallo di confidenza (IC), compreso tra 3288 g e 3346 g, rappresenta l’incertezza sulla media attesa del peso per neonati con caratteristiche simili, mentre l’intervallo di predizione (IP), compreso tra 2779 g e 3855 g, indica il range entro cui è plausibile che cada il peso effettivo di un singolo neonato di questo tipo. La differenza tra IC e IP riflette la normale variabilità individuale: l’intervallo più ampio del IP mostra che, pur avendo una stima media piuttosto precisa, i pesi dei singoli neonati possono variare sensibilmente. Nel complesso, la previsione è coerente con il comportamento del modello (RMSE ≈ 270–280 g) e con valori realistici per un neonato a termine, confermando una buona capacità predittiva e attendibilità statistica del modello stimato.

# Modello 
m_vis <- lm(peso ~ gestazione * fumatrici, data = df)
summary(m_vis)

## 
## Call:
## lm(formula = peso ~ gestazione * fumatrici, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1609.03  -289.03   -11.54   280.97  1898.44 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -3240.771    178.861 -18.119   <2e-16 ***
## gestazione              167.495      4.585  36.534   <2e-16 ***
## fumatrici1             1343.426   1164.746   1.153    0.249    
## gestazione:fumatrici1   -36.764     29.646  -1.240    0.215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 422.9 on 2496 degrees of freedom
## Multiple R-squared:  0.352,  Adjusted R-squared:  0.3513 
## F-statistic:   452 on 3 and 2496 DF,  p-value: < 2.2e-16

# Predizioni su griglia gestazione × fumatrici
gest_vals <- seq(min(df$gestazione, na.rm = TRUE),
                 max(df$gestazione, na.rm = TRUE), by = 1)
fum_lvls  <- levels(df$fumatrici)
grid <- expand.grid(gestazione = gest_vals, fumatrici = fum_lvls)
grid$yhat <- predict(m_vis, newdata = grid)

# Grafico con due linee (una per livello di fumo)
par(mfrow = c(1,1), mar = c(4,4,2,1))
plot(NA,
     xlim = range(gest_vals),
     ylim = range(grid$yhat, na.rm = TRUE),
     xlab = "Settimane di gestazione",
     ylab = "Peso previsto (g)",
     main = "Peso previsto vs Gestazione per livello di fumo (modello semplice)")
for (fv in fum_lvls) {
  sel <- grid$fumatrici == fv
  lines(grid$gestazione[sel], grid$yhat[sel], lwd = 2)
}
legend("topleft", legend = fum_lvls, lwd = 2, bty = "n", title = "Fumatrici")

Il grafico rappresenta l’andamento del peso previsto dei neonati in funzione delle settimane di gestazione, con distinzione tra madri non fumatrici (0) e fumatrici (1). In entrambi i gruppi si osserva una relazione positiva: al crescere delle settimane di gestazione il peso aumenta, come fisiologicamente atteso. Tuttavia, la pendenza delle due linee non è identica e si incrociano attorno alla 35ª settimana: ciò suggerisce che l’effetto del fumo non è costante lungo tutto il periodo gestazionale. Nelle prime settimane di gestazione, i neonati di madri fumatrici mostrano un peso previsto inferiore rispetto a quelli di madri non fumatrici, indicando un potenziale ritardo di crescita intrauterina. Nelle ultime settimane, la differenza tende a ridursi fino quasi a invertirsi, ma il peso previsto rimane comunque leggermente inferiore per le fumatrici. Questo comportamento conferma l’interazione significativa tra fumo e gestazione: l’effetto del fumo dipende dalla durata della gravidanza e risulta più marcato nelle gestazioni più brevi.

CONCLUSIONI: Tipo di parto ~ Ospedale: nessuna evidenza di associazione (test χ² → non si rifiuta H₀). Confronto con popolazione: peso campione allineato a ~3300 g; lunghezza leggermente inferiore, differenza significativa ma piccola. Differenze tra sessi: significative per peso, lunghezza e cranio (M > F). Modello (stepwise): R² ~ 0,73; RMSE ~ 273 g → buona capacità esplicativa/predittiva. Effetti principali (β): gestazione (~+32 g/settimana), sesso (M vs F ~+77 g), lunghezza (~+10,3 g/mm), cranio (~+10,5 g/mm), n_gravidanze (piccolo effetto positivo).