(In fondi risposta al feedback del Tutor)

Dettagli del Progetto

1) Raccolta dei Dati e Struttura del Dataset

DATASET: Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali.

OBIETTIVO: L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature

knitr::opts_chunk$set(echo = TRUE)

# Carico Dataset
dataset <- read.csv("neonati.csv")
# Anni.madre con un minimo di zero non è plausibile.
dati <- subset(dataset, Anni.madre > 12)
n <- nrow(dati)
attach(dati)

#Carico Librerie
library(lmtest)
library(TeachingDemos)
library(car)
library(knitr)
library(visreg)
library(ggplot2)
# consigliati da ChatGPT per fare tabelle con summary in HTML
library(broom)
library(knitr)

2) Analisi e Modellizzazione

2.1) Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

# Funzione analisi Variabili Continue
DescribeVar <- function(var,LabelVar="Variabile") {
  cat("\n", LabelVar, "\n", rep("-", nchar(LabelVar)), "\n", sep="")
  cat("Numero di NA:", sum(is.na(var)), "\n")
  # outlier IQR
  Q1 <- quantile(var, 0.25, na.rm=TRUE)
  Q3 <- quantile(var, 0.75, na.rm=TRUE)
  med <- median(var, na.rm=TRUE)
  IQR_val <- Q3 - Q1
  outlier_idx <- which(var < Q1 - 1.5*IQR_val | var > Q3 + 1.5*IQR_val)
  cat("Outlier (IQR):", length(outlier_idx), "su", sum(!is.na(var)), "(", round(100*length(outlier_idx)/sum(!is.na(var)), 2), "% )\n")
  # grafici
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar), add = TRUE)
  par(mfrow=c(1,2))
  hist(var, breaks=30, main=paste("Istogramma:", LabelVar), xlab=LabelVar)
  abline(v=c(Q1, med, Q3), col="red", lwd=2, lty=2)
  boxplot(var, main=paste("Boxplot:", LabelVar), ylab=LabelVar)
}

# Funzione analisi Variabili categoriali
DescribeCat <- function(var,LabelVar="Variabile") {
  print(knitr::kable(table(var), caption = LabelVar))
}



# Variabili Continue
DescribeVar(Peso,LabelVar="Peso (g)")
## 
## Peso (g)
## --------
## Numero di NA: 0 
## Outlier (IQR): 69 su 2498 ( 2.76 % )

DescribeVar(Lunghezza,LabelVar="Lunghezza (mm)")
## 
## Lunghezza (mm)
## --------------
## Numero di NA: 0 
## Outlier (IQR): 59 su 2498 ( 2.36 % )

DescribeVar(Cranio,LabelVar="Diametro Cranio (mm)")
## 
## Diametro Cranio (mm)
## --------------------
## Numero di NA: 0 
## Outlier (IQR): 48 su 2498 ( 1.92 % )

DescribeVar(Anni.madre,LabelVar="Anni madre (anni)")
## 
## Anni madre (anni)
## -----------------
## Numero di NA: 0 
## Outlier (IQR): 11 su 2498 ( 0.44 % )

# Variabili Categoriali
Fumatrici_2 <- factor(Fumatrici,levels = c(0, 1),labels = c("No", "Sì"))
Tipo.parto_2 <- factor(Tipo.parto,levels = c("Nat", "Ces"),labels = c("Naturale", "Cesareo"))
DescribeCat(Fumatrici_2,LabelVar="Fumatrici")
## 
## 
## Table: Fumatrici
## 
## |var | Freq|
## |:---|----:|
## |No  | 2394|
## |Sì  |  104|
DescribeCat(Tipo.parto_2,LabelVar="Tipo di parto")
## 
## 
## Table: Tipo di parto
## 
## |var      | Freq|
## |:--------|----:|
## |Naturale | 1770|
## |Cesareo  |  728|
DescribeCat(Ospedale,LabelVar="Ospedale")
## 
## 
## Table: Ospedale
## 
## |var  | Freq|
## |:----|----:|
## |osp1 |  816|
## |osp2 |  848|
## |osp3 |  834|
DescribeCat(Sesso,LabelVar="Sesso")
## 
## 
## Table: Sesso
## 
## |var | Freq|
## |:---|----:|
## |F   | 1255|
## |M   | 1243|
DescribeCat(N.gravidanze,LabelVar="Numero gravidanze")
## 
## 
## Table: Numero gravidanze
## 
## |var | Freq|
## |:---|----:|
## |0   | 1095|
## |1   |  817|
## |2   |  340|
## |3   |  150|
## |4   |   48|
## |5   |   21|
## |6   |   11|
## |7   |    1|
## |8   |    8|
## |9   |    2|
## |10  |    3|
## |11  |    1|
## |12  |    1|
DescribeCat(Gestazione,LabelVar="Gestazione")
## 
## 
## Table: Gestazione
## 
## |var | Freq|
## |:---|----:|
## |25  |    1|
## |26  |    1|
## |27  |    2|
## |28  |    4|
## |29  |    3|
## |30  |    5|
## |31  |    8|
## |32  |    9|
## |33  |   18|
## |34  |   16|
## |35  |   33|
## |36  |   62|
## |37  |  192|
## |38  |  437|
## |39  |  580|
## |40  |  741|
## |41  |  328|
## |42  |   56|
## |43  |    2|

2.1) Analisi Preliminare

Inoltre si saggeranno le seguenti ipotesi con i test adatti:

  1. in alcuni ospedali si fanno più parti cesarei

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

  3. Le misure antropometriche sono significativamente diverse tra i due sessi

# 1. ----------- test -> in alcuni ospedali si fanno più parti cesarei?
# creare una tabella incrociata tipo di parto vs Ospedale e creare una matrice dei valori attesi 
tab <- table(Tipo.parto_2, Ospedale)
test.indipendenza <- chisq.test(tab)
kable(round(test.indipendenza$expected, 2), caption = "Frequenze attese (test Chi-quadro)")
Frequenze attese (test Chi-quadro)
osp1 osp2 osp3
Naturale 578.19 600.86 590.94
Cesareo 237.81 247.14 243.06 
kable(data.frame(statistic = round(unname(test.indipendenza$statistic), 2),df = unname(test.indipendenza$parameter),p_value = signif(test.indipendenza$p.value, 3)), caption = "Risultato test Chi-quadro (Pearson)")
Risultato test Chi-quadro (Pearson)
statistic df p_value
1.08 2 0.582

CONCLUSIONE 2.1.1: Il test Chi2 Pearson ci dice che non si può statisticamente rifiutare l’ipotesi di indipendenza (pvalue>0.05) suggerendo che la distribuzione di frequenze di parti naturali/cesarei non differisce in modo significativo tra i tre ospedali.

——– 2. test -> quanto media peso/lunghezza significativamente uguali a quelli della pop. mondiale? Ho trovato questo studio come riferimento https://www.scielo.br/j/spmj/a/5HBjGPShrJ4nYpPgQRDKt7m/?format=pdf&lang=en (Table 1) peso medio pop mondiale: mean= 3215.4 g sd= 488.5 g lughezza media pop mondiale: mean= 493 mm sd= 22 mm

——-Test t, H0: Peso = Peso mondiale ——–

BW_mondiale <- rnorm(10000000,mean=3215.4,sd=488.5)
plot(density(Peso),col=1,lwd=2,xlab="Peso alla nascita (g)", main="Confronto densità peso alla nascita")
lines(density(BW_mondiale),col=4,lwd=2)
legend("topright",legend=c("Pop. campione","Pop. mondiale"),col=c(1, 4),lwd=2,bty="n")

tt_peso <- t.test(Peso, mu=3215.4)
kable(data.frame(mean_sample = round(mean(Peso), 2),mu0 = 3215.4,t = round(unname(tt_peso$statistic), 2),df=round(unname(tt_peso$parameter), 0),p_value = signif(tt_peso$p.value, 3)), caption = "t-test Peso vs riferimento mondiale")
t-test Peso vs riferimento mondiale
mean_sample mu0 t df p_value
3284.18 3215.4 6.55 2497 0

——-Test t, H0: Lunghezza = Lunghezza mondiale ——–

BL_mondiale <- rnorm(10000000,mean=493,sd=22)
plot(density(Lunghezza),col=1,lwd=2, xlab = "Lunghezza (mm)",main = "Confronto densità lunghezza alla nascita")
lines(density(BL_mondiale),col =4,lwd=2)
legend("topright",legend =c("Pop. campione","Pop. mondiale"),col=c(1, 4),lwd=2,bty="n")

tt_lunghezza <- t.test(Lunghezza, mu=493)
kable(data.frame(mean_sample = round(mean(Lunghezza), 2),mu0 = 493,t = round(unname(tt_lunghezza$statistic), 2),df = round(unname(tt_lunghezza$parameter), 0),p_value = signif(tt_lunghezza$p.value, 3)), caption = "t-test Lunghezza vs riferimento mondiale")
t-test Lunghezza vs riferimento mondiale
mean_sample mu0 t df p_value
494.7 493 3.22 2497 0.0013

CONCLUSIONE 2.1.2: I t-test indicano che le medie di Peso e Lunghezza osservate nel campione differiscono in modo statisticamente significativo dai rispettivi valori di riferimento riportati nello studio utilizzato (p-value<0.05), almeno entro l’incertezza campionaria del nostro campione.

#—— 3. misure antropometriche diverse tra i due sessi? ———

#oldpar <- par(no.readonly = TRUE)
#on.exit(par(oldpar), add = TRUE)

# t-test Peso, Lunghezza e Cranio
tt_peso <- t.test(Peso ~ Sesso, alternative = "two.sided", conf.level = 0.95)
tt_lung <- t.test(Lunghezza ~ Sesso, alternative = "two.sided", conf.level = 0.95)
tt_cranio <- t.test(Cranio ~ Sesso, alternative = "two.sided", conf.level = 0.95)

tt_tab <- data.frame(Variabile = c("Peso (g)", "Lunghezza (mm)", "Cranio (mm)"),
  t_statistic = round(c(tt_peso$statistic,tt_lung$statistic,tt_cranio$statistic), 3),
  df = round(c(tt_peso$parameter,tt_lung$parameter,tt_cranio$parameter), 1),
  p_value = signif(c(tt_peso$p.value,tt_lung$p.value,tt_cranio$p.value), 3),
  CI_low = round(c(tt_peso$conf.int[1],tt_lung$conf.int[1],tt_cranio$conf.int[1]), 2),
  CI_high = round(c(tt_peso$conf.int[2],tt_lung$conf.int[2],tt_cranio$conf.int[2]), 2),check.names = FALSE)

kable(tt_tab,caption = "Confronto tra Maschi e Femmine (t-test a due campioni)")
Confronto tra Maschi e Femmine (t-test a due campioni)
Variabile t_statistic df p_value CI_low CI_high
Peso (g) -12.115 2488.7 0 -287.48 -207.38
Lunghezza (mm) -9.582 2457.3 0 -11.94 -7.88
Cranio (mm) -7.437 2489.4 0 -6.11 -3.56 
par(mfrow=c(1,3))
boxplot(Peso ~ Sesso)
boxplot(Lunghezza ~ Sesso)
boxplot(Cranio ~ Sesso)

CONCLUSIONE 2.1.3: In tutti i casi i t test restituiscono p-value<0.05, così da rifiutare H0 (uguaglianza delle medie tra maschi e femmine). In media, Peso, Lunghezza e Diametro del cranio risultano maggiori nei neonati di sesso maschile rispetto a quelli di sesso femminile.

2.2) Creazione del Modello di Regressione

Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

2.3) Selezione del Modello Ottimale

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

# Voglio solo valori numerici ed escludo variabili non generalizzabili per previsione  Peso, come ospedale e tipo.parto
dati_mod <- subset(dati, select = -c(Ospedale, Tipo.parto))
dati_num <- dati_mod[,sapply(dati_mod, is.numeric)]

# covarianza =sum((Peso-mean(Peso))*(Gestazione-mean(Gestazione))/(n-1)
covarianza_PesoGest <- cov(Peso,Gestazione) 
#coefficente di correlazione lineare di Pearson
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 = cex.cor * r)}

par(mfrow=c(1,1))
pairs(dati_num,upper.panel = panel.smooth,lower.panel = panel.cor)

kable(round(cor(dati_num),2))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
Anni.madre 1.00 0.38 0.01 -0.13 -0.02 -0.06 0.02
N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06 0.04
Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02 -0.01
Gestazione -0.13 -0.10 0.03 1.00 0.59 0.62 0.46
Peso -0.02 0.00 -0.02 0.59 1.00 0.80 0.70
Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00 0.60
Cranio 0.02 0.04 -0.01 0.46 0.70 0.60 1.00

Dalla matrice di correlazione noto delle correlazioni positive tra peso e gestazione, lunghezza e diametro del cranio, ed una lieve correlazione positiva tra anni della madre e numero di gravidanze, che però è abbastanza ovvia come correlazione e non inerente al nostro scopo. Per il resto non abbiamo correlazioni evidenti.

################### modello lineare completo ######################
mod_all <- lm(Peso ~., data=dati_mod)
kable(tidy(mod_all, conf.int = TRUE),digits = 3,caption = "Modello Lineare con tutti i parametri")
Modello Lineare con tutti i parametri
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6712.240 141.334 -47.492 0.000 -6989.385 -6435.096
Anni.madre 0.880 1.149 0.766 0.444 -1.373 3.134
N.gravidanze 11.379 4.677 2.433 0.015 2.208 20.549
Fumatrici -30.396 27.608 -1.101 0.271 -84.533 23.741
Gestazione 32.947 3.829 8.605 0.000 25.439 40.455
Lunghezza 10.232 0.301 33.979 0.000 9.641 10.822
Cranio 10.520 0.427 24.633 0.000 9.682 11.357
SessoM 78.079 11.213 6.963 0.000 56.091 100.067 
kable(glance(mod_all),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.727 0.726 274.707 948.29 0 7 -17568.53 35155.07 35207.48 187905214 2490 2498 
################## (A) rimuovo solo Anni.madre (Pr(>|t|): 0.444) #################
mod1 <- update(mod_all,~.-Anni.madre)
kable(tidy(mod1, conf.int = TRUE),digits = 3,caption = "Modello iniziale (escludendo Anni.madre)")
Modello iniziale (escludendo Anni.madre)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6682.264 135.798 -49.207 0.000 -6948.553 -6415.975
N.gravidanze 12.700 4.347 2.921 0.004 4.176 21.224
Fumatrici -30.573 27.605 -1.108 0.268 -84.703 23.558
Gestazione 32.644 3.808 8.573 0.000 25.177 40.111
Lunghezza 10.231 0.301 33.979 0.000 9.640 10.821
Cranio 10.537 0.426 24.707 0.000 9.700 11.373
SessoM 78.160 11.212 6.971 0.000 56.174 100.145 
kable(glance(mod1),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.727 0.726 274.684 1106.424 0 6 -17568.83 35153.66 35200.24 187949505 2491 2498 
################## (B) rimuovo solo Fumatrici  (Pr(>|t|): 0.26818) ##################
mod2 <- update(mod_all,~.-Fumatrici)
kable(tidy(mod2, conf.int = TRUE),digits = 3,caption = "Modello iniziale (escludendo Fumatrici)")
Modello iniziale (escludendo Fumatrici)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6712.065 141.340 -47.489 0.000 -6989.221 -6434.910
Anni.madre 0.891 1.149 0.775 0.438 -1.362 3.144
N.gravidanze 11.120 4.671 2.381 0.017 1.961 20.280
Gestazione 32.691 3.822 8.554 0.000 25.197 40.186
Lunghezza 10.246 0.301 34.058 0.000 9.656 10.836
Cranio 10.524 0.427 24.642 0.000 9.687 11.361
SessoM 77.900 11.212 6.948 0.000 55.913 99.887 
kable(glance(mod2),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.727 0.726 274.719 1106.042 0 6 -17569.14 35154.28 35200.87 187996688 2491 2498 
################## (C) rimuovo sia Fumatrici che  Anni.madre ##################
mod3 <- update(mod1,~.-Fumatrici)
kable(tidy(mod3, conf.int = TRUE),digits = 3,caption = "Modello iniziale (escludendo sia Anni.madre che Fumatrici)")
Modello iniziale (escludendo sia Anni.madre che Fumatrici)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6681.725 135.804 -49.201 0.000 -6948.025 -6415.426
N.gravidanze 12.455 4.342 2.869 0.004 3.942 20.969
Gestazione 32.383 3.801 8.520 0.000 24.930 39.836
Lunghezza 10.245 0.301 34.059 0.000 9.656 10.835
Cranio 10.541 0.426 24.717 0.000 9.705 11.377
SessoM 77.981 11.211 6.956 0.000 55.997 99.965 
kable(glance(mod3),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.727 0.726 274.697 1327.343 0 5 -17569.44 35152.89 35193.65 188042054 2492 2498 
### (D) modello ulteriormente ridotto escludendo N.gravidanze per valutare un modello più semplice ###
mod4 <- update(mod3,~.-N.gravidanze)
kable(tidy(mod4, conf.int = TRUE),digits = 3,caption = "Modello iniziale (escludendo sia Anni.madre che Fumatrici che N.gravidanze)")
Modello iniziale (escludendo sia Anni.madre che Fumatrici che N.gravidanze)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6651.673 135.595 -49.055 0 -6917.564 -6385.782
Gestazione 31.326 3.788 8.269 0 23.897 38.755
Lunghezza 10.202 0.301 33.909 0 9.612 10.792
Cranio 10.671 0.425 25.126 0 9.838 11.503
SessoM 79.103 11.220 7.050 0 57.100 101.105 
kable(glance(mod4),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.726 0.726 275.095 1652.329 0 4 -17573.56 35159.12 35194.06 188663107 2493 2498 
# In tutti i casi Adjusted R-squared non cambia significativamente (~0.7265) 

################### 2.3) Selezione del Modello Ottimale ##################
kable(BIC(mod_all,mod1,mod2,mod3,mod4))
df BIC
mod_all 9 35207.48
mod1 8 35200.24
mod2 8 35200.87
mod3 7 35193.65
mod4 6 35194.06 
kable(AIC(mod_all,mod1,mod2,mod3,mod4))
df AIC
mod_all 9 35155.07
mod1 8 35153.66
mod2 8 35154.28
mod3 7 35152.89
mod4 6 35159.12 
kable(anova(mod_all,mod1,mod2,mod3,mod4))
Res.Df RSS Df Sum of Sq F Pr(>F)
2490 187905214 NA NA NA NA
2491 187949505 -1 -44291.73 0.5869257 0.4436830
2491 187996688 0 -47182.11 NA NA
2492 188042054 -1 -45366.21 0.6011640 0.4382079
2493 188663107 -1 -621053.31 8.2298022 0.0041555

Il confronto tra modelli tramite il BIC (in accordo con ANOVA e AIC) mostra una preferenza per il modello più semplice, che esclude i parametri Fumatrici e Anni.madre, con una differenza di oltre 6 rispetto ai modelli alternativi, indicando che è il modello migliore da considerare secondo il principio del rasoio di Occam.

# Selezione automatica delle covarianze tramite procedura stepwise con criterio BIC
stepwise.mod_best<-MASS::stepAIC(mod_all,direction="both",k=log(n))

Start: AIC=28110.64 Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Sesso

           Df Sum of Sq       RSS   AIC

Step: AIC=28103.4 Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Sesso

           Df Sum of Sq       RSS   AIC

Step: AIC=28096.81 Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso

           Df Sum of Sq       RSS   AIC

188042054 28097 - N.gravidanze 1 621053 188663107 28097 + Fumatrici 1 92548 187949505 28103 + Anni.madre 1 45366 187996688 28104 - Sesso 1 3650790 191692844 28137 - Gestazione 1 5477493 193519547 28161 - Cranio 1 46098547 234140601 28637 - Lunghezza 1 87532691 275574744 29044

# Output del modello selezionato (formula finale e stime)
kable(tidy(stepwise.mod_best, conf.int = TRUE),digits = 3,caption = "BEST MODEL")
BEST MODEL
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -6681.725 135.804 -49.201 0.000 -6948.025 -6415.426
N.gravidanze 12.455 4.342 2.869 0.004 3.942 20.969
Gestazione 32.383 3.801 8.520 0.000 24.930 39.836
Lunghezza 10.245 0.301 34.059 0.000 9.656 10.835
Cranio 10.541 0.426 24.717 0.000 9.705 11.377
SessoM 77.981 11.211 6.956 0.000 55.997 99.965 
kable(glance(stepwise.mod_best),digits = 3,caption = "Statistiche globali del BEST MODEL")
Statistiche globali del BEST MODEL
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.727 0.726 274.697 1327.343 0 5 -17569.44 35152.89 35193.65 188042054 2492 2498

COMMENTI ALLE STIME DEI COEFF beta PER IL BEST MODEL

  • A parità delle altre variabili, il N.gravidanze è associato a un aumento medio del peso alla nascita di circa 12.5g per ogni gravidanza aggiuntiva.
  • La durata della Gestazione mostra ogni settimana in più è associata a un incremento medio del peso di ~32.4g.
  • La Lunghezza alla nascita è fortemente associata al Peso, con un aumento medio di circa 10.2g per ogni mm aggiuntivo.
  • Un incremento di 1 mm del Cranio corrisponde a un aumento medio del peso di ~10.5g.
  • Infine, a parità delle altre condizioni, i neonati di sesso maschile presentano un peso medio alla nascita superiore di ~78g rispetto alle femmine.

Test per considerare regressioni non-lineari (excursus formale)

Per completezza, verifico se semplici estensioni del modello lineare (usando termini quadratici o interazioni) producano un miglioramento sostanziale (in accordo con il principio di marginalità). Nella scelta finale, seppur mi baso sul confronto tramite BIC, cerco di dare priorità alla semplicità del modello (rasoio di Occam) e all’interpretabilità inferenziale.

# Modello con termini quadratici 
mod_quad <- update(stepwise.mod_best, ~ . + I(Gestazione^2) + I(Lunghezza^2))
kable(tidy(mod_quad, conf.int = TRUE),digits = 3,caption = "Modello con termini quadratici per Gestazione e Lunghezza")
Modello con termini quadratici per Gestazione e Lunghezza
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -2359.379 905.751 -2.605 0.009 -4135.481 -583.278
N.gravidanze 14.466 4.249 3.405 0.001 6.134 22.798
Gestazione 336.175 62.739 5.358 0.000 213.149 459.201
Lunghezza -32.123 4.039 -7.953 0.000 -40.043 -24.203
Cranio 10.447 0.419 24.909 0.000 9.625 11.270
SessoM 72.599 11.007 6.596 0.000 51.015 94.183
I(Gestazione^2) -3.868 0.825 -4.689 0.000 -5.486 -2.251
I(Lunghezza^2) 0.044 0.004 10.545 0.000 0.036 0.052 
kable(glance(mod_quad),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.739 0.739 268.587 1008.397 0 7 -17512.25 35042.5 35094.91 179625530 2490 2498 
# Modello con una sola interazione 
mod_int <- update(stepwise.mod_best, ~ . + Cranio * Lunghezza)
kable(tidy(mod_int, conf.int = TRUE),digits = 3,caption = "Modello con termini interazione Cranio * Lunghezza")
Modello con termini interazione Cranio * Lunghezza
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -1803.011 1018.128 -1.771 0.077 -3799.476 193.454
N.gravidanze 12.933 4.323 2.991 0.003 4.455 21.410
Gestazione 38.149 3.967 9.616 0.000 30.370 45.929
Lunghezza -0.306 2.203 -0.139 0.890 -4.626 4.014
Cranio -4.755 3.192 -1.490 0.136 -11.014 1.505
SessoM 73.240 11.204 6.537 0.000 51.270 95.210
Lunghezza:Cranio 0.032 0.007 4.835 0.000 0.019 0.044 
kable(glance(mod_int),digits = 3,caption = "Statistiche globali del modello")
Statistiche globali del modello
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.73 0.729 273.472 1119.946 0 6 -17557.78 35131.56 35178.14 186293990 2491 2498 
# Confronto BIC: modello lineare vs estensioni
kable(BIC(stepwise.mod_best, mod_quad, mod_int))
df BIC
stepwise.mod_best 7 35193.65
mod_quad 9 35094.91
mod_int 8 35178.14

COMMENTI ALL’ACCETTAZIONE O MENO DELL’AUMENTO DELLA COMPLESSITÀ

  1. confronto con BIC: il confronto tramite BIC indica una preferenza per il modello con termini quadratici “mod_quad.” Quindi, dal punto di vista statistico, una lieve curvatura può descrivere meglio alcune regioni dei dati.

  2. miglioramento del fit: rispetto a stepwise.mod_best, c’è un incremento di Adjusted R-squared di 0.012 ed una riduzione del Residual standard error è di circa 6 g. Questo è un miglioramento contenuto rispetto alla scala “Peso”.

  3. complessità dell’interpretazione: un altro punto importante riguarda la perdita di immediatezza nell’interpretazione dei coefficienti principali. Infatti, in presenza dei termini quadratici, i coefficienti lineari di Gestazione e Lunghezza cambiano segno: questo non significa che “sparisce” l’effetto, ma che l’effetto non è più costante e va interpretato come variabile al variare di Gestazione e Lunghezza. In pratica il modello diventa meno leggibile e più complicato da spiegare.

  4. confronto con i grafici: visivamente le relazioni appaiono monotone e senza punti di flesso evidenti nei grafici in “pairs”. Questo non nega che una lieve curvatura possa aiutare il fit (come infatti suggerisce il BIC), ma indica che la non-linearità non è così marcata da rendere indispensabile un modello complesso per raccontare i risultati.

Conclusione e scelta finale: dato che l’obiettivo principale del progetto è interpretare in modo semplice e comunicabile i fattori associati alla variabile risposta “Peso”, mantengo come best model quello lineare “stepwise.mod_best”, privilegiando semplicità (rasoio di Occam) e interpretabilità inferenziale rispetto a un miglioramento di fit contenuto rispetto alla scala del “Peso” e ottenuto al prezzo di una maggiore complessità interpretativa.

2.4) Analisi della Qualità del Modello

Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R2 e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.

par(mfrow=c(2,2))
plot(stepwise.mod_best)

- Residuals vs Fitted: la nube dei punti è complessivamente centrata intorno a 0, ma si nota una leggera struttura per valori bassi dei fitted (possibile non-linearità locale o sottogruppo di osservazioni).
- Q–Q Residuals: i punti sono abbastanza allineati nella parte centrale, mentre nelle code (soprattutto superiore) si osservano deviazioni evidenti dalla normalità, coerenti con la presenza di residui estremi (e.g., 1551, 155, 1306).
- Scale-Location: la dispersione non è perfettamente costante lungo i fitted (andamento non del tutto piatto), suggerendo eteroschedasticità, coerentemente con il test di Breusch–Pagan che rifiuta H0.
- Residuals vs Leverage: è presente almeno un’osservazione con residuo standardizzato molto alto e leverage relativamente elevato (1551), associata a Cook’s distance molto grande (vicina alla curva di riferimento 1). È quindi un candidato “influente” e va verificato.

#Usiamo un metodo numerico per indagare valori di Leverage e Outliers

#Leverage (hat-values)
par(mfrow=c(1,1))
lev <- hatvalues(stepwise.mod_best)
plot(lev, ylab="Leverage (hat-values)",main="hat values")
p <- length(coef(stepwise.mod_best))
soglia <- 2*p/n
abline(h=soglia, col=2)

lev_oltresoglia <- lev[lev>soglia]
# Vediamo che percentuali di valori sono sopra la soglia
infl <- which(lev > soglia)
nlev <- length(infl)
perclev <- round(100*nlev/n, 1)
cat("% Leverage sopra la soglia =", perclev, "%\n")
## % Leverage sopra la soglia = 6.1 %
cat("Num Leverage sopra la soglia =", nlev, "\n")
## Num Leverage sopra la soglia = 152
# Nota: una quota ~6% sopra soglia indica presenza di punti con leverage relativamente alto, che non implica influenza sulle stime.

# Outlier nei residui 
plot(rstudent(stepwise.mod_best),main="Outlier nei residui ")
abline(h=c(-2,2), col=2)

outls <- outlierTest(stepwise.mod_best)
# Estrazione della matrice 
outls_mat <- outls$rstudent
outls_tab <- as.data.frame(round(outls_mat, 3))
kable(outls_tab, caption = "Outlier nei residui")
Outlier nei residui
round(outls_mat, 3)
1551 10.046
155 5.025
1306 4.825 
# 3 outliers: 1306, 155, 1551 (estremamente anomalo con rstudent~10)


#distanza di cook
cookD <- cooks.distance(stepwise.mod_best)
plot(cookD,main = paste("Max CookDistance =", round(max(cookD), 3)),ylab="Cook Distance")

# Un valore massimo ~0.83 è elevato (vicino alla curva di riferimento 1 nel grafico diagnostico), quindi l'osservazione corrispondente merita verifica. 


# Test diagnostici sui residui
bptest(stepwise.mod_best)
## 
##  studentized Breusch-Pagan test
## 
## data:  stepwise.mod_best
## BP = 90.297, df = 5, p-value < 2.2e-16
# bptest: rifiuto H0 -> evidenza di eteroschedasticità (varianza non costante)

dwtest(stepwise.mod_best)
## 
##  Durbin-Watson test
## 
## data:  stepwise.mod_best
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0
# dwtest: non rifiuto H0 -> nessuna evidenza di autocorrelazione positiva dei residui

shapiro.test(residuals(stepwise.mod_best))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(stepwise.mod_best)
## W = 0.97414, p-value < 2.2e-16
# shapiro.test: rifiuto H0 -> residui formalmente non normali con test molto sensibile a n grande 

# Check visivo (supporto, non decisionale)
plot(density(residuals(stepwise.mod_best)), main="Densità dei residui")

# Deviazioni soprattutto nelle code, coerentemente con il Q-Q plot.


# CASO ANOMALO VISUALIZZATO in Peso~4370 g Lunghezza~315 mm 
plot(Peso, Lunghezza, pch=20, xlab="Peso (g)", ylab="Lunghezza (mm)", main="Peso vs Lunghezza")
# in basso a destra
x_out <- 4370
y_out <- 315
symbols(x_out, y_out,circles = 100,inches = FALSE,add = TRUE,fg = "red",lwd = 2)
text(x_out, y_out, labels = "Leverage", pos = 3, col = "red")

Commento finale: il modello appare complessivamente adeguato. Dalla diagnostica emerge evidenza di eteroschedasticità e deviazioni dalla normalità soprattutto nelle code, coerenti con pochi residui estremi. Le stime beta dei coefficienti restano una sintesi utile delle relazioni principali, ma i risultati inferenziali (p-value e intervalli di confidenza) vanno interpretati con cautela. Inoltre, alcune osservazioni (in particolare 1551) mostrano residui estremi e potenziale influenza (Cook’s distance elevata), per cui è opportuno verificarne la plausibilità e valutare l’impatto sulle stime (per esempio si potrebbe confrontare il modello con e senza questi punti).

3) Previsioni e Risultati

Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

# Dall'esempio ora capisco che probabilmente il numero di gravidanze è un parametro importante per il fit.

L39 <- median(Lunghezza[Gestazione == 39])
C39 <- median(Cranio[Gestazione == 39])

#Predizione: N.gravidanze=3,Sesso=F,Gestazione=39,Lunghezza=500 mm,Cranio=340 mm
caso_test <- data.frame(N.gravidanze=3,Sesso="F",Gestazione=39,Lunghezza=L39,Cranio=C39)
#PESO (g):
predict(stepwise.mod_best,newdata=caso_test,interval="prediction", level=0.95)
##        fit      lwr      upr
## 1 3325.219 2786.036 3864.402

4) Visualizzazioni

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the visreg package.
##   Please report the issue at <https://github.com/pbreheny/visreg/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
##   Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Conclusioni

Nel complesso, il modello fornisce una descrizione interpretabile e statisticamente solida dei principali fattori associati al Peso dei neonati alla nascita, per il campione (statisticamente scarso rispetto quello mondiale) analizzato.

risposta al feedback del Tutor

Ciao Andrea! Ok usare la programmazione, ma nell’analisi descrittiva non serve fare alcune cose un tanto al chilo. Ad esempio, saggiare la normalità per ogni variabile.

Ho ridotto l’analisi descrittiva finalizzata all’identificazione di outlier o anomalie

potevi usare solo il t-test (lo z test è didattico)

Fatto

tipo di parto e ospedale potevano essere escluse a priori “logicamente” dal modello in quanto non hanno molto senso in ottica previsionale

Ok grazie. Sinceramente non ho escluso a priori queste variabili in quanto l’ospedale potrebbe essere correlato ad una diversa locazione delle nascite, che presumo influirebbe sulle caratteristiche fisiche del bambino, mentre il tipo di parto potrebbe ci sta, in quanto, anche nel caso mostrasse una qualche correlazione potrebbe essere indiretta, e direttamente correlata alla gestazione.

considerare l’effetto quadratico (o d’interazione) senza considerare uno o più effetti principali è una questione molto spinosa su cui la gente ci si scanna nei forum. In ottica di statistica pura e non ML/DL, direi che ci si rifà a questo: https://en.wikipedia.org/wiki/Principle_of_marginality

Chiarissimo grazie per l’intervento. Ho preso in considerazione tale effetto di marginalizzazione.

pertanto il modello scelto non andrebbe bene.

Chiaro.

Il processo di selezione del modello è comunque molto confusinario e può essere snellito di molto. Se provi così tanti modelli e così tanti regressori ce ne saranno alcuni significativi solo per effetto del caso ma senza effetto reale per poter fare inferenza. RASOIO DI OCCAM!

Ho scelto il modello lineare semplice. La capacità di individuare modelli più complessi richiede probabilmente un campione statisticamente molto più ampio.

Mancano i commenti alle stime (coefficienti beta) in termini di “incremento medio di Y al variare di X”

Commenti inseriti solo per il best fit

si chiamano “code” della distribuzione, non “ali”

Ok. Mi sono fatto confondere dal termine in inglesi wings :)

Qualità dell’analisi L’analisi non è male ma è presentata male. Migliora la leggibilità, commenta meglio i passaggi e sopratuttto smeplifica il processo di modellazione. Meno modelli e spiegati più chiaramente.

Hai ragione, sono stato molto superficiale considerando che faccio questo master nel tempo libero in cui non lavoro. Mi sono concentrato sul comprendere la teoria piuttosto che l’esposizione. Questa volta ho riguardato il testo, spero sia chiaro.

Qualità del documento Ci sono dei print di console grezzi che non stanno benissimo… sempre meglio inserire i dati almeno in una tabella semplice per sistemare l’output, tipo kable()

Ho riguardato l’output HTML. Ho visto che summary non andava bene, e mi sono fatto consigliare da chatgpt come poterlo printare meglio in HTML. Per il resto, quando possibile, ho usato kable() come consigliato.

i messaggi e i warnings della console, come per l’importazione di librerie e altro, si possono anche omettere per fare tutto più pulito

Fatto. Ho preferito mettere {r setup, include=TRUE, message=FALSE, warning=FALSE}

Arrotonda a 2 cifre decimali, o più solo quando serve

Fatto. I p-values sono stati approssimati a 3 cifre decimali.