Progetto: Un Modello Statistico per Prevedere il Peso dei Neonati

Svolgimento

1. Analisi e Modellizzazione

1.1. Analisi Preliminare

Carico il dataset

dati <- read.csv("neonati.csv", sep = ",", stringsAsFactors = T) 
attach(dati)
head(dati,4)

Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	0	42	3380	490	325	Nat	osp3	M
21	2	0	39	3150	490	345	Nat	osp1	F
34	3	0	38	3640	500	375	Nat	osp2	M
28	1	0	41	3690	515	365	Nat	osp2	M

la classificazione delle variabili è

summary_tbl <- tibble(
  variabile = c("Anni.madre", "N.gravidanze", "Fumatrici", "Gestazione", "Peso",
                "Lunghezza", "Cranio", "Tipo.parto", "Ospedale", "Sesso"),
  tipo = c("Quantitativa continua","Quantitativa discreta",
           "Qualitativa dicotomica", "Quantitativa continua (riportata come discreta)",
           "Quantitativa continua", "Quantitativa continua",
           "Quantitativa continua", "Qualitativa nominale (dicotomica)", 
           "Qualitativa nominale (politomica)", "Qualitativa dicotomica")
)

kable(summary_tbl, col.names = c("Variabile","Tipo statistico"))

Variabile	Tipo statistico
Anni.madre	Quantitativa continua
N.gravidanze	Quantitativa discreta
Fumatrici	Qualitativa dicotomica
Gestazione	Quantitativa continua (riportata come discreta)
Peso	Quantitativa continua
Lunghezza	Quantitativa continua
Cranio	Quantitativa continua
Tipo.parto	Qualitativa nominale (dicotomica)
Ospedale	Qualitativa nominale (politomica)
Sesso	Qualitativa dicotomica

Ora, per ogni variabile, facciamo una statistica descrittiva

info_indici <- function(var){
  tibble(
    Nome_var  = as.character(ensym(var)),
    Media     = round(mean(var),2),
    Mediana   = round(median(var),2),
    Sd        = round(sd(var),2),
    Skewness  = round(skewness(var),2),
    Kurtosis  = round(kurtosis(var) - 3,2),
  )
}

freq_table <- function(dati, var, freq_cum = "no"){
  ft <- dati %>%
    group_by(!!ensym(var)) %>%
    summarise(count_class = n(), .groups = "drop") %>%
    ungroup() %>%
    mutate(count_tot = sum(count_class),
           freq_rel = round(count_class / count_tot,3))
  if(freq_cum == "yes"){
    ft <- ft %>%
      mutate(freq_cum = cumsum(freq_rel))
  }
  return(ft) 
}

info_indici(Anni.madre)

Nome_var	Media	Mediana	Sd	Skewness	Kurtosis
Anni.madre	28.16	28	5.27	0.04	0.38

distribuzione quasi perfettamente simmetrica e leggermente leptocurtica.

freq_table(dati, N.gravidanze, freq_cum = "yes") 

N.gravidanze	count_class	count_tot	freq_rel	freq_cum
0	1096	2500	0.438	0.438
1	818	2500	0.327	0.765
2	340	2500	0.136	0.901
3	150	2500	0.060	0.961
4	48	2500	0.019	0.980
5	21	2500	0.008	0.988
6	11	2500	0.004	0.992
7	1	2500	0.000	0.992
8	8	2500	0.003	0.995
9	2	2500	0.001	0.996
10	3	2500	0.001	0.997
11	1	2500	0.000	0.997
12	1	2500	0.000	0.997

più dei tre quarti delle madri sono alla prima o alla seconda gravidanza.

freq_table(dati, Fumatrici) 

Fumatrici	count_class	count_tot	freq_rel
0	2396	2500	0.958
1	104	2500	0.042

solamente il \(4.2%\) delle madri in gravidanze fumano.

info_indici(Gestazione)

Nome_var	Media	Mediana	Sd	Skewness	Kurtosis
Gestazione	38.98	39	1.87	-2.07	8.26

Media e mediana sono molto simili, quindi non c’è una grande distorsione, ma si osserva sia una forte asimmetria negativa che una forte leptocurticità per la distribuzione.

info_indici(Peso) 

Nome_var	Media	Mediana	Sd	Skewness	Kurtosis
Peso	3284.08	3300	525.04	-0.65	2.03

anche qui media e mediana molto vicini, ma si osserva una distribuzione moderatamente asimmetrica (negativa) e una leptocurticità.

info_indici(Lunghezza)

Nome_var	Media	Mediana	Sd	Skewness	Kurtosis
Lunghezza	494.69	500	26.32	-1.51	6.49

la distribuzione della lunghezza presenta una forte asimmetria negativa e un’elevata leptocurtosi, analogamente a quanto osservato per la variabile Gestazione. Questo risultato è coerente con l’ipotesi che le due variabili possano essere correlate, in quanto entrambe legate allo sviluppo prenatale.

info_indici(Cranio)

Nome_var	Media	Mediana	Sd	Skewness	Kurtosis
Cranio	340.03	340	16.43	-0.79	2.95

anche la variabile Cranio ha una asimmetria negativa (ma in questo caso è moderata) ed è leptocurtica. Media e mediana praticamente coincidenti.

freq_table(dati, Tipo.parto) 

Tipo.parto	count_class	count_tot	freq_rel
Ces	728	2500	0.291
Nat	1772	2500	0.709

poco meno del 30% dei neonati sono nati con un parto cesareo.

freq_table(dati, Ospedale)

Ospedale	count_class	count_tot	freq_rel
osp1	816	2500	0.326
osp2	849	2500	0.340
osp3	835	2500	0.334

il dataset contiene le informazioni sulle nascite di neonati in 3 ospedali. I record sembrano essere sostanzialmente equiripartiti su di essi.

freq_table(dati, Sesso)

Sesso	count_class	count_tot	freq_rel
F	1256	2500	0.502
M	1244	2500	0.498

il numero di neonate e maggiore dei neonati, seppur di poco.

Facciamo ora delle osservazioni grafiche delle distribuzioni:

p1 <- ggplot(dati, aes(x = "", y = Anni.madre)) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Età madre", y = "Anni", x = "") +
  theme_minimal()

p2 <- ggplot(dati, aes(x = N.gravidanze)) +
  geom_histogram(fill = "lightblue", color = "black", bins = 13) +
  labs(title = "Numero gravidanze", x = "N.gravidanze") +
  theme_minimal()

p3 <- ggplot(dati, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
  geom_boxplot() +
  labs(x = "Fumo materno",
       y = "Peso alla nascita (g)",
       fill = "Fumo materno",
       title = "Fumatrici e non") +
  scale_fill_manual(values = c("0" = "steelblue", "1" = "tomato"),
                    labels = c("Non fumatrice", "Fumatrice")) +
  theme_minimal()

p4 <- ggplot(dati, aes(x = Gestazione)) +
  geom_histogram(fill = "lightblue", color = "black", bins = 19) +
  labs(title = "Durata gestazione", x = "Settimane") +
  theme_minimal()

p5 <- ggplot(dati, aes(x = "", y = Peso)) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Peso alla nascita", y = "Grammi", x = "") +
  theme_minimal()

p6 <- ggplot(dati, aes(x = "", y = Lunghezza)) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Lunghezza", y = "mm", x = "") +
  theme_minimal()

p7 <- ggplot(dati, aes(x = "", y = Cranio)) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Diametro cranio", y = "mm", x = "") +
  theme_minimal()

(p1 | p2 | p4) /
(p5 | p6 | p7) /
(p3)

per la variabilità Anni.madre sembrano esserci dei record anomali in quanto riportano un età della madre di pochi anni, identifichiamoli

dati[dati$Anni.madre < 10,]

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1152	1	1	0	41	3250	490	350	Nat	osp2	F
1380	0	0	0	39	3060	490	330	Nat	osp3	M

sono 2 e probabilmente sono frutto di una errata trascrizione. Decido di eliminarli dal campione

dati <- dati %>%
  filter(dati$Anni.madre >= 10)
suppressMessages(attach(dati))

Saggiamo ora alcune ipotesi che ci vengono richieste:

In alcuni ospedali si fanno più parti cesarei?

ballon_osp_parto <- dati %>%
  count(Ospedale, Tipo.parto)   

ggplot(ballon_osp_parto, aes(x = Tipo.parto, y = Ospedale)) +
  geom_point(aes(size = n), shape = 21, color="black", fill="skyblue") +
  geom_text(aes(label = n), color="black", size=3) +  
  scale_size_area(max_size = 15) +
  theme_minimal() +
  labs(title="Balloon Plot Ospedale vs Tipo di parto", size="Count")

tab_osp_parto <- table(Ospedale, Tipo.parto)
chisq.test(tab_osp_parto)

## 
##  Pearson's Chi-squared test
## 
## data:  tab_osp_parto
## X-squared = 1.083, df = 2, p-value = 0.5819

Poiché il p-value è \(0.58\) (maggiore del solito \(5\%\) della soglia di significatività), non rifiuto l’ipotesi nulla concludendo che non ci sono evidenze statistiche per affermare che in alcuni ospedali avvengono più parti cesarei rispetto ad altri.

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

Da una ricerca è stato trovato che il peso medio alla nascita dei neonati è \(\mu_{Peso} = 3300 \; g\), mentre la loro lunghezza al momento del parto è \(\mu_{Lung} = 500 \; mm\), per cui

t.test_Peso <- t.test(Peso, mu = 3300)
cat("p-value Peso:", round(t.test_Peso$p.value,3), "\n")

## p-value Peso: 0.132

t.test_Lunghezza <-t.test(Lunghezza, mu = 500)
cat("p-value Lunghezza:", round(t.test_Lunghezza$p.value,3), "\n")

## p-value Lunghezza: 0

concludiamo che nel caso del Peso tramite un t-test si ottiene un p-value di circa \(0.13\) che è maggiore del valore di significatività del 5%, si accetta allora l’ipotesi nulla “il peso medio dei neonati del campione non è diverso dal peso medio dei neonati della popolazione”, mentre per la Lunghezza si ottiene un p-value nullo, quindi si scarta l’ipotesi nulla e si può affermare che “la lunghezza dei neonati del campione è significativamente diverso dalla lunghezza media della popolazione”.

Le misure antropometriche sono significativamente diverse tra i due sessi?

Peso_Sesso_ttest <- t.test(Peso~Sesso, data = dati)
Lunghezza_Sesso_ttest <- t.test(Lunghezza~Sesso, data = dati)
Cranio_Sesso_ttest <- t.test(Cranio~Sesso, data = dati)

cat(" p-value Peso/Sesso:", Peso_Sesso_ttest$p.value, "\n",
    "p-value Lunghezza/Sesso:", Lunghezza_Sesso_ttest$p.value, "\n",
    "p-value Cranio/Sesso:", Cranio_Sesso_ttest$p.value, "\n")

##  p-value Peso/Sesso: 7.275684e-33 
##  p-value Lunghezza/Sesso: 2.232328e-21 
##  p-value Cranio/Sesso: 1.414019e-13

Per tutti e tre i t-test si osserva un p-value nullo, per cui si rifiuta l’ipotesi nulla e si puo affermare che tutte e tre le misure antropometriche differiscono significativamente tra neonati maschi e neonate femmine.

1.2. Creazione del Modello di Regressione

Prima di tutto osserviamo la matrice di correlazione

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)
}

pairs(dati, upper.panel = panel.smooth, lower.panel = panel.cor) 

Come ci si aspettava, le variabili Gestazione, Lunghezza e Cranio sono fortemente correlate positivamente con la variabile peso (soprattutto la variabile Lunghezza). Mentre le variabili Anni.madre e N.gravidanze sono debolmente correlate la peso del neonato.

Facciamo ora un focus per le variabili qualitative. Partiamo dalla variabile Fumatrici e testiamo l’ipotesi nulla secondo cui il peso medio dei neonati alla nascita non differisce tra madri fumatrici e non fumatrici. Per farlo utilizziamo un t-test per campioni indipendenti:

ttest_Peso.Fum <- t.test(Peso~Fumatrici)
ttest_Peso.Fum$p.value

## [1] 0.3022785

non si rifiuta l’ipotesi nulla visto che il pvalue è maggiore della soglia di significatività del \(5 \%\). Si conclude che non ci sono evidenze statisticamente significative per affermare che il peso dei neonati differisca tra madri fumatrici e non fumatrici.

Passiamo alla variabile Tipo.parto e facciamo un t-test

ttest_Peso.TipParto <- t.test(Peso~Tipo.parto)
ttest_Peso.TipParto$p.value

## [1] 0.8916349

anche qui non si rifiuta l’ipotesi nulla (pvalue \(> 5 \%\)), non ci sono evidenze statisticamente significative per affermare che il peso dei neonati differisca tra parti cesarei e naturali.

Per la variabile Ospedale conduciamo test ANOVA a una via, si vuole confrontare il peso medio dei neonati nei tre ospedali.

anova_osp <- aov(Peso~Ospedale)
summary(anova_osp)

##               Df    Sum Sq Mean Sq F value Pr(>F)
## Ospedale       2    948538  474269    1.72  0.179
## Residuals   2495 687888603  275707

Il test fornisce un p-value di \(0.18\), maggiore della soglia di significatività fissato al \(5\%\). Pertanto non si rifiuta l’ipotesi nulla e si conclude che non vi sono evidenze statistiche per affermare che il peso medio differisca tra i tre ospedali.

Infine, passiamo alla variabile Sesso e con un t-test verifichiamo se il peso medio dei neonati di sesso maschile è statisticamente differente con il peso dei neonati di sesso femminile

ttest_Peso.Sesso <- t.test(Peso~Sesso)
ttest_Peso.Sesso$p.value

## [1] 7.275684e-33

allora si rifiuta l’ipotesi nulla ed è possibile affermare che statisticamente il peso dei neonati maschi e femmine differiscono significativamente alla nascita.

Splittiamo il dataset in \(train\_data\), campione su cui verrano construiti modelli, e \(test\_data\), campione che permetterà di misurare la bonta del best_model

set.seed(123)
train_index <- createDataPartition(dati$Peso, p = 0.9, list = FALSE)
train_data <- dati[train_index, ]
test_data  <- dati[-train_index, ]

Creiamo ora un modello di regrassione lineare full, quindi che tenga conto di tutte le variabili del campione

mod_full <- lm(Peso~., data = train_data)
summary(mod_full)

## 
## Call:
## lm(formula = Peso ~ ., data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1119.18  -184.60   -15.72   163.59  2603.16 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6695.2255   151.1718 -44.289  < 2e-16 ***
## Anni.madre        0.9697     1.2252   0.792   0.4287    
## N.gravidanze      9.7108     4.9632   1.957   0.0505 .  
## Fumatrici       -33.4629    28.6038  -1.170   0.2422    
## Gestazione       32.9770     4.0643   8.114 7.99e-16 ***
## Lunghezza        10.2558     0.3214  31.914  < 2e-16 ***
## Cranio           10.3769     0.4522  22.947  < 2e-16 ***
## Tipo.partoNat    27.0790    12.9186   2.096   0.0362 *  
## Ospedaleosp2    -14.5641    14.4189  -1.010   0.3126    
## Ospedaleosp3     26.6961    14.5162   1.839   0.0660 .  
## SessoM           69.8884    11.9911   5.828 6.41e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278.5 on 2240 degrees of freedom
## Multiple R-squared:  0.7199, Adjusted R-squared:  0.7187 
## F-statistic: 575.8 on 10 and 2240 DF,  p-value: < 2.2e-16

anova(mod_full)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Anni.madre	1	222127.46	222127.46	2.8635499	0.0907464
N.gravidanze	1	21274.89	21274.89	0.2742646	0.6005371
Fumatrici	1	239522.71	239522.71	3.0878002	0.0790180
Gestazione	1	218426416.26	218426416.26	2815.8379464	0.0000000
Lunghezza	1	181850044.01	181850044.01	2344.3146816	0.0000000
Cranio	1	42170949.16	42170949.16	543.6455944	0.0000000
Tipo.parto	1	359509.01	359509.01	4.6346002	0.0314399
Ospedale	2	692854.57	346427.29	4.4659575	0.0115962
Sesso	1	2635039.85	2635039.85	33.9695415	0.0000000
Residuals	2240	173758285.00	77570.66	NA	NA

vif(mod_full)

##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.198384  1        1.094707
## N.gravidanze 1.198202  1        1.094624
## Fumatrici    1.007891  1        1.003938
## Gestazione   1.686855  1        1.298790
## Lunghezza    2.069725  1        1.438654
## Cranio       1.614286  1        1.270546
## Tipo.parto   1.003632  1        1.001815
## Ospedale     1.004007  2        1.001000
## Sesso        1.043084  1        1.021315

Valutiamo le stime sui coefficienti associati ai regressori:

• Anni.madre: Per ogni anno in più della madre il bambino alla nascita ha un peso di \(1 g\) in più, praticamente la variabile non è influente (come anche il p-value ci mostra).

• N.gravidanze: Per ogni gravidanze precendeti avuta dalla madre si ha alla nascita un peso medio aggiuntivo di \(9.7 g\).

• Fumatrici: Bambini nati da madri fumatrici hanno in media un peso di \(33 g\) in meno rispetto a bambini nati da madri non fumatrici.

• Gestazione: Per ogni settimana di gestazione c’è un aumento medio del peso del neonato di \(33 g\).

• Lunghezza: Per ogni millimetro di lunghezza in più del neonato c’è un aumento medio del peso di \(10.3 g\).

• Cranio: Per ogni millimetro di diametro in più del neonato c’è un aumento medio del peso di \(10.4 g\).

• Tipo.parto: Neonati nati da parti cesarei pesano in media \(29.33 g\) in meno rispetto ai neonati nati da parto naturale.

• Ospedale: Neonati nati nell’ospedale 2 pesano in media \(14.6 g\) in meno rispetto a quelli nati nell’ospedale 1, mentre quelli nati nell’ospedale 3 pesano in media \(26.7 g\) in più rispetto all’ospedale 1.

• Sesso: I neonati pesano in media \(69 g\) in più rispetto alle neonate.

Per il modello full si ha \(R^2_{adj} = 0.7187\), buon valore, ma migliorabile. Inoltre i VIF sono bel al di sotto della soglia indicativa di \(5\), quindi c’è una bassa multicollinearità tra le variabili.

Creiamo un modello scartando le variabili Anni.madre e Fumatrici

mod2 <- update(mod_full,
               ~.-Anni.madre -Fumatrici, 
               data = train_data)
summary(mod2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso, data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1106.94  -185.61   -17.06   163.51  2612.11 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6662.3723   145.6151 -45.753  < 2e-16 ***
## N.gravidanze     10.8955     4.5917   2.373   0.0177 *  
## Gestazione       32.3482     4.0348   8.017 1.72e-15 ***
## Lunghezza        10.2709     0.3211  31.988  < 2e-16 ***
## Cranio           10.4022     0.4515  23.038  < 2e-16 ***
## Tipo.partoNat    26.9261    12.9164   2.085   0.0372 *  
## Ospedaleosp2    -14.0815    14.4145  -0.977   0.3287    
## Ospedaleosp3     27.4640    14.5050   1.893   0.0584 .  
## SessoM           69.7974    11.9890   5.822 6.66e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278.5 on 2242 degrees of freedom
## Multiple R-squared:  0.7197, Adjusted R-squared:  0.7187 
## F-statistic: 719.4 on 8 and 2242 DF,  p-value: < 2.2e-16

il valore \(R^2_{adj}\) rimane uguale, c’è un miglioramente di significatività della variabile N.gravidanze e in generale ci sono variazioni sugli stimatori del modello.

Costruiamo un modello che tanga conto solo di variabili molto significative, come: N.gravidaze, Gestazione, Lunghezza, Cranio, Ospedale e Sesso

mod3 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Ospedale + Sesso,
           data = train_data)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Ospedale + Sesso, data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1124.78  -187.54   -17.32   166.07  2613.37 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6641.1593   145.3674 -45.685  < 2e-16 ***
## N.gravidanze    10.6146     4.5931   2.311   0.0209 *  
## Gestazione      32.4044     4.0377   8.025 1.62e-15 ***
## Lunghezza       10.2410     0.3210  31.903  < 2e-16 ***
## Cranio          10.4334     0.4516  23.102  < 2e-16 ***
## Ospedaleosp2   -14.2468    14.4251  -0.988   0.3234    
## Ospedaleosp3    27.8103    14.5149   1.916   0.0555 .  
## SessoM          69.8230    11.9979   5.820 6.74e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278.7 on 2243 degrees of freedom
## Multiple R-squared:  0.7191, Adjusted R-squared:  0.7182 
## F-statistic: 820.4 on 7 and 2243 DF,  p-value: < 2.2e-16

Facciamo ora delle considerazioni generali sulle variabili significative che abbiamo deciso di tenere. Dagli scatterplot si osservano andamenti non lineari e inoltre sappiamo che le settimane di gestazione influiscono positivamente sulle variabili Lunghezza e Cranio.

Costruisco un modello che tenga conto delle interazioni sopra mezionati

mod4 <- lm(Peso~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso 
           + Gestazione:Cranio, 
           data = train_data)
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Cranio, data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1130.16  -185.10   -15.34   168.34  2693.02 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        296.02782 1163.04478   0.255  0.79911    
## N.gravidanze        11.72767    4.56490   2.569  0.01026 *  
## Gestazione        -152.39216   31.04740  -4.908 9.84e-07 ***
## Lunghezza           10.45764    0.32136  32.542  < 2e-16 ***
## Cranio             -11.34019    3.65222  -3.105  0.00193 ** 
## SessoM              63.77135   11.97031   5.327 1.10e-07 ***
## Gestazione:Cranio    0.57045    0.09489   6.011 2.14e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277 on 2244 degrees of freedom
## Multiple R-squared:  0.7225, Adjusted R-squared:  0.7217 
## F-statistic: 973.6 on 6 and 2244 DF,  p-value: < 2.2e-16

\(R^2_{adj}\) migliorato rispetto ai modelli precedenti, è stato osservato che l’interazione Gestazione:Lunghezza porta ad un deterioramento della variabile Lunghezza.

Costruisco un modello con termini quadratici per le variabili Gestazione e Lunghezza

mod5 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso 
           + I(Gestazione^2) + I(Lunghezza^2), 
           data = train_data)
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2) + I(Lunghezza^2), data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1189.11  -184.02   -12.92   169.23  1334.89 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.111e+03  9.508e+02  -2.220   0.0265 *  
## N.gravidanze     1.283e+01  4.488e+00   2.858   0.0043 ** 
## Gestazione       3.545e+02  6.574e+01   5.393 7.67e-08 ***
## Lunghezza       -3.444e+01  4.236e+00  -8.130 7.00e-16 ***
## Cranio           1.034e+01  4.434e-01  23.323  < 2e-16 ***
## SessoM           6.570e+01  1.176e+01   5.587 2.59e-08 ***
## I(Gestazione^2) -4.106e+00  8.655e-01  -4.744 2.22e-06 ***
## I(Lunghezza^2)   4.607e-02  4.346e-03  10.598  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.2 on 2243 degrees of freedom
## Multiple R-squared:  0.732,  Adjusted R-squared:  0.7312 
## F-statistic: 875.4 on 7 and 2243 DF,  p-value: < 2.2e-16

\(R^2_{adj}\) migliorato sensisibilmente rispetto al modello mod4. È stato notato che un termine quadratico sulla variabile Cranio porta ad un suo deterioramento.

Costruiamo un modello che tenga conto di termini quadratici e di interazione

mod6 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso 
           + Gestazione:Lunghezza + I(Lunghezza^2), 
           data = train_data)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Lunghezza + I(Lunghezza^2), data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1211.34  -185.76   -11.57   170.08  1321.56 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -2.298e+03  9.527e+02  -2.412  0.01594 *  
## N.gravidanze          1.267e+01  4.485e+00   2.825  0.00478 ** 
## Gestazione            2.764e+02  4.647e+01   5.948 3.14e-09 ***
## Lunghezza            -2.756e+01  3.492e+00  -7.892 4.60e-15 ***
## Cranio                1.024e+01  4.453e-01  22.987  < 2e-16 ***
## SessoM                6.672e+01  1.176e+01   5.673 1.59e-08 ***
## I(Lunghezza^2)        5.838e-02  6.034e-03   9.675  < 2e-16 ***
## Gestazione:Lunghezza -4.873e-01  9.674e-02  -5.037 5.10e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.1 on 2243 degrees of freedom
## Multiple R-squared:  0.7324, Adjusted R-squared:  0.7315 
## F-statistic: 876.9 on 7 and 2243 DF,  p-value: < 2.2e-16

si osserva un \(R^2_{adj}\) praticamente invariato rispetto al modello mod5.

Costruisco un modello con termine cubico per la variabile Lunghezza

mod7 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso 
           + Gestazione:Lunghezza + I(Lunghezza^2) + I(Lunghezza^3), 
           data = train_data)
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Lunghezza + I(Lunghezza^2) + I(Lunghezza^3), 
##     data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1213.55  -185.78   -13.56   167.04  1316.06 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.796e+04  4.542e+03   3.954 7.92e-05 ***
## N.gravidanze          1.290e+01  4.465e+00   2.889  0.00391 ** 
## Gestazione            2.964e+02  4.647e+01   6.378 2.18e-10 ***
## Lunghezza            -1.671e+02  3.080e+01  -5.426 6.38e-08 ***
## Cranio                1.016e+01  4.437e-01  22.908  < 2e-16 ***
## SessoM                6.619e+01  1.171e+01   5.652 1.79e-08 ***
## I(Lunghezza^2)        3.717e-01  6.895e-02   5.390 7.79e-08 ***
## I(Lunghezza^3)       -2.289e-04  5.018e-05  -4.561 5.38e-06 ***
## Gestazione:Lunghezza -5.341e-01  9.686e-02  -5.514 3.90e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.9 on 2242 degrees of freedom
## Multiple R-squared:  0.7348, Adjusted R-squared:  0.7339 
## F-statistic: 776.7 on 8 and 2242 DF,  p-value: < 2.2e-16

ancora una volta si osserva leggero miglioramente rispetto al modello mod6.

Infine, proviamo a costruire un modello lineare generalizzato che trasforma la variabile risposta applicando una trasformazione Box-Cox al modello mod7. Prima di tutto cerco il valore ottimale lambda per la varaibile Peso

bc <- boxcox(mod7, lambda = seq(-3, 3, 0.001), plotit = FALSE)
max_index <- which.max(bc$y)
lambda_opt <- bc$x[max_index]

cat("lambda ottimale:", lambda_opt)

## lambda ottimale: 0.581

if(lambda_opt == 0){
  Peso_bc <- log(train_data$Peso)
} else {
  Peso_bc <- (train_data$Peso^lambda_opt - 1)/lambda_opt
}
  
mod_lambda <- lm(Peso_bc ~ N.gravidanze + Gestazione + Lunghezza + Cranio 
                 + Sesso + Gestazione:Lunghezza + I(Lunghezza^2) + I(Lunghezza^3), 
                 data = train_data)
summary(mod_lambda)

## 
## Call:
## lm(formula = Peso_bc ~ N.gravidanze + Gestazione + Lunghezza + 
##     Cranio + Sesso + Gestazione:Lunghezza + I(Lunghezza^2) + 
##     I(Lunghezza^3), data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.750  -6.103  -0.292   5.754  47.185 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           6.305e+02  1.528e+02   4.126 3.82e-05 ***
## N.gravidanze          4.109e-01  1.502e-01   2.735  0.00629 ** 
## Gestazione            1.350e+01  1.564e+00   8.636  < 2e-16 ***
## Lunghezza            -6.106e+00  1.036e+00  -5.892 4.39e-09 ***
## Cranio                3.449e-01  1.493e-02  23.103  < 2e-16 ***
## SessoM                2.212e+00  3.940e-01   5.613 2.24e-08 ***
## I(Lunghezza^2)        1.450e-02  2.320e-03   6.250 4.91e-10 ***
## I(Lunghezza^3)       -9.348e-06  1.688e-06  -5.536 3.45e-08 ***
## Gestazione:Lunghezza -2.508e-02  3.259e-03  -7.696 2.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.114 on 2242 degrees of freedom
## Multiple R-squared:  0.7556, Adjusted R-squared:  0.7547 
## F-statistic: 866.3 on 8 and 2242 DF,  p-value: < 2.2e-16

sensibile miglioramento a seguito della trasformazione.

1.3. Selezione del Modello Ottimale

Tramite le metrice AIC e BIC vogliamo selezionare il modello migliore

AIC(mod_full, mod2, mod3, mod4, mod5, mod6, mod7, mod_lambda)

	df	AIC
mod_full	12	31744.92
mod2	10	31742.94
mod3	9	31745.30
mod4	8	31716.31
mod5	9	31639.27
mod6	9	31636.43
mod7	10	31617.64
mod_lambda	10	16347.43

BIC(mod_full, mod2, mod3, mod4, mod5, mod6, mod7, mod_lambda)

	df	BIC
mod_full	12	31813.55
mod2	10	31800.14
mod3	9	31796.78
mod4	8	31762.06
mod5	9	31690.75
mod6	9	31687.90
mod7	10	31674.83
mod_lambda	10	16404.62

entrambe le metriche suggeriscono che il modello migliore sia mod_lambda. Sarà il nostro modello di riferimento.

1.4. Analisi della Qualità del Modello

Valutiamo la capacità predittiva del modello sul test_data

pred_test <- predict(mod_lambda, newdata = test_data)
pred_test <- (lambda_opt * pred_test + 1)^(1/lambda_opt)
rmse_test <- rmse(test_data$Peso, pred_test)

calcoliamo la RSS (Residual Sum of Squares) e TSS (Total Sum of Squares)

rss <- sum((test_data$Peso - pred_test)^2)
tss <- sum((test_data$Peso - mean(test_data$Peso))^2)

da cui calcoliamo \(R^2\) (coefficiente di determinazione) e \(R^2_{adj}\) sul test_data

R2_test <- 1 - rss/tss

n <- nrow(test_data)
p <- length(coef(mod_lambda)) - 1 
R2_adj_test <- 1 - (1 - R2_test) * (n - 1) / (n - p - 1)

cat(" R^2 adj sul test set:", round(R2_adj_test, 3),"\n",
    "RMSE sul test set:", round(rmse_test, 0))

##  R^2 adj sul test set: 0.801 
##  RMSE sul test set: 231

un valore \(R^2_{adj} \simeq 80 \%\) significa che il modello riesca a spiegare circa il \(80 \%\) della variabilità del peso neonatale su test set, per cui è un buon modello predittivo soprattutto se si pensa che riguarda dati bioloci. Un valore di \(RMSE \simeq 231 g\) sta ad indicare che in media l’errore sulla previzione del peso è di circa \(231 g\), con un errore medio rispetto al peso medio \(< 10 \%\).

Facciamo ora un’analisi sui residui. Poiché la variabile target Peso risultava moderatamente asimmetrica, questa asimmetria potrebbe riflettersi sui residui e compromettere alcune delle condizioni necessarie per il modello lineare. Facciamo osservazioni grafiche:

plot(mod_lambda, cex = 0.25)

• Il primo grafico mostra effettivamente una distribuzione senza particolari pattern attorno alla media nulla.

• Nel secondo grafico vengono messi in relazione i residui con i quantili di una distribuzione normale. I punti si dispongono sulla bisettrice tranne che per le code, quindi i residui non seguono perfettamente una distribuzione normale (nelle code).

• Nel terzo grafico, utilie per testare l’eteroschedasticità, non si osservano pattern evidenti. I punti si distribuiscono grossomodo orizontalmente attorno ad un valore di y che indica una varianza abbastanza costante.

• Nel quarto graifico, si visualizzano i potenziali valori influenti, ovvero quei valori molto leverage o molto outliers o una combinazioned dei due. In generale quando alcuni residui mostrano una distanza di Cook che supera le linee trattegiate vengono considerate “residui di osservazioni potenzialmente influenti sulle stime di regressioni”, qui se ne osserva uno problematico, osserviamolo

cooksd <- cooks.distance(mod_lambda)
train_data[which(cooksd > 1),]

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1549	35	1	0	38	4370	315	374	Nat	osp3	F

Peso ben oltre la media, lunghezza molto piccola, non sembra errore di trascrizione ma un record di un neonato con possibili patologie, ma non ne possiamo essere sicuri.

Si vuole osservare come la presenza di questo record influenzi il modello

mod_lambda_no_outlier <- update(mod_lambda, subset = -which(cooksd > 1))
summary(mod_lambda_no_outlier)

## 
## Call:
## lm(formula = Peso_bc ~ N.gravidanze + Gestazione + Lunghezza + 
##     Cranio + Sesso + Gestazione:Lunghezza + I(Lunghezza^2) + 
##     I(Lunghezza^3), data = train_data, subset = -which(cooksd > 
##     1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.977  -6.069  -0.313   5.754  49.193 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.485e+02  1.659e+02   2.101 0.035767 *  
## N.gravidanze          4.095e-01  1.497e-01   2.736 0.006264 ** 
## Gestazione            8.227e+00  1.987e+00   4.141 3.59e-05 ***
## Lunghezza            -3.741e+00  1.171e+00  -3.194 0.001422 ** 
## Cranio                3.380e-01  1.496e-02  22.602  < 2e-16 ***
## SessoM                2.194e+00  3.925e-01   5.590 2.54e-08 ***
## I(Lunghezza^2)        9.242e-03  2.617e-03   3.531 0.000423 ***
## I(Lunghezza^3)       -6.061e-06  1.849e-06  -3.278 0.001060 ** 
## Gestazione:Lunghezza -1.427e-02  4.114e-03  -3.467 0.000536 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.079 on 2241 degrees of freedom
## Multiple R-squared:  0.7572, Adjusted R-squared:  0.7563 
## F-statistic: 873.4 on 8 and 2241 DF,  p-value: < 2.2e-16

I coefficienti sono cambiati, in alcuni casi anche di parecchio. Questo è normale perché l’outlier influenzava fortemente i coefficienti, quindi rimuovendolo la stima è diventata più “robusta”. In generale si osserva una diminuzione della significativita delle variabili.

pred_test_no_outliers <- predict(mod_lambda_no_outlier, newdata = test_data)
pred_test_no_outliers <- (lambda_opt * pred_test_no_outliers + 1)^(1/lambda_opt)
rmse_test_no <- rmse(test_data$Peso, pred_test_no_outliers)

rss_no <- sum((test_data$Peso - pred_test_no_outliers)^2)
tss_no <- sum((test_data$Peso - mean(test_data$Peso))^2)
R2_test_no <- 1 - rss_no/tss_no

n <- nrow(test_data)
p <- length(coef(mod_lambda_no_outlier)) - 1 
R2_adj_test_no <- 1 - (1 - R2_test_no) * (n - 1) / (n - p - 1)

cat(" R^2 adj sul test set:", round(R2_adj_test_no, 3),"\n",
    "RMSE sul test set:", round(rmse_test_no, 0))

##  R^2 adj sul test set: 0.801 
##  RMSE sul test set: 231

Invece \(R^2_{adj}\) e \(RMSE\) sul test_set rimangono uguali al modello precedente, quindi le capacità predittive del modello non sono cambiate.

Osserviamo i punti di leverage

lev <- hatvalues(mod_lambda_no_outlier)
plot(lev)
p = sum(lev)
soglia = 2*p/n
abline(h=soglia, col=2)

lev[lev>soglia]

##        928       1617       1778       2435       2450 
## 0.32212433 0.11164435 0.17532778 0.19147740 0.09979501

cinque non superano la soglia. Analizziamo i punti di outlier

plot(rstudent(mod_lambda_no_outlier))
abline(h = c(-2,2), col = 2)

outlierTest(mod_lambda_no_outlier)

##       rstudent unadjusted p-value Bonferroni p
## 155   5.507841         4.0487e-08   9.1095e-05
## 1397 -4.543006         5.8403e-06   1.3141e-02
## 1692  4.454743         8.8126e-06   1.9828e-02
## 1305  4.395716         1.1556e-05   2.6002e-02

Vengono osserti solo quattro outliers significativi. Verifichiamo la distanza di Cook

cook <- cooks.distance(mod_lambda_no_outlier)
plot(cook)

Si conclude che sono prenseti 5 leverage significativi, 4 outliers significativi ma nessun record ha una distanza di Cook problematico. Quindi il modello mod_lambda_no_outlier è robusto.

2. Previsioni e Risultati

Facciamo delle previsioni con il modello mod_lambda, visto che il record problematico non intacca le predizioni

nuovo_neonato1 <- data.frame(N.gravidanze = 3, Fumatrici = "0", Gestazione = 39,
                             Lunghezza = 500, Cranio = 340, Tipo.parto = "Nat",
                             Ospedale = "osp2", Sesso = "F" ) 
nuovo_neonato2 <- data.frame(N.gravidanze = 0, Fumatrici = "0", Gestazione = 39,
                             Lunghezza = 500, Cranio = 340, Tipo.parto = "Nat",
                             Ospedale = "osp2", Sesso = "F" ) 
nuovo_neonato3 <- data.frame(N.gravidanze = 0, Fumatrici = "0",Gestazione = 36, 
                             Lunghezza = 490, Cranio = 350, Tipo.parto = "Ces",
                             Ospedale = "osp3", Sesso = "M" ) 

pred_1 <- predict(mod_lambda, newdata = nuovo_neonato1) 
pred_2 <- predict(mod_lambda, newdata = nuovo_neonato2) 
pred_3 <- predict(mod_lambda, newdata = nuovo_neonato3) 

pred1 <- (lambda_opt*pred_1+1)^(1/lambda_opt) 
pred2 <- (lambda_opt*pred_2+1)^(1/lambda_opt) 
pred3 <- (lambda_opt*pred_3+1)^(1/lambda_opt) 

cat(" Predizione peso neonato 1:", round(pred1,0), "g \n",
    "Predizione peso neonato 2:", round(pred2,0), "g \n",
    "Predizione peso neonato 3:", round(pred3,0), "g")

##  Predizione peso neonato 1: 3324 g 
##  Predizione peso neonato 2: 3287 g 
##  Predizione peso neonato 3: 3229 g

Vediamo esplicitamente le previsioni su alcuni record del test_set (per cui conosciamo già il valore del peso)

pred_test_1 <- predict(mod_lambda, newdata = test_data[1,]) 
pred_test_2 <- predict(mod_lambda, newdata = test_data[2,]) 
pred_test_3 <- predict(mod_lambda, newdata = test_data[3,]) 

pred_test1 <- (lambda_opt*pred_test_1+1)^(1/lambda_opt) 
pred_test2 <- (lambda_opt*pred_test_2+1)^(1/lambda_opt) 
pred_test3 <- (lambda_opt*pred_test_3+1)^(1/lambda_opt) 

cat(" Neonato 1 test: Peso reale", round(test_data[1,]$Peso,0), "g, Peso stimato",
    round(pred_test1,0), "g, Differenza:",
    round(abs(test_data[1,]$Peso - pred_test1),0), "g \n",
    "Neonato 2 test: Peso reale", round(test_data[2,]$Peso,0), "g, Peso stimato",
    round(pred_test2,0), "g, Differenza:",
    round(abs(test_data[2,]$Peso - pred_test2),0), "g \n",
    "Neonato 3 test: Peso reale", round(test_data[3,]$Peso,0), "g, Peso stimato",
    round(pred_test3,0), "g, Differenza:",
    round(abs(test_data[3,]$Peso - pred_test3),0), "g")

##  Neonato 1 test: Peso reale 3380 g, Peso stimato 3190 g, Differenza: 190 g 
##  Neonato 2 test: Peso reale 2950 g, Peso stimato 2675 g, Differenza: 275 g 
##  Neonato 3 test: Peso reale 3050 g, Peso stimato 2801 g, Differenza: 249 g

l’errore per questi record è confrontabile con RMSE calcolato.

Valutiamo il record del test_set che presenta RMSE più grande e osserviamolo per trarre qualche conclusione

errors <- abs(test_data$Peso - pred_test)

worst_idx <- which.max(errors)
worst_idx

## 2193 
##  221

test_data[worst_idx, ]

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
2193	38	1	0	40	3980	480	335	Nat	osp1	F

cat(" Peso predetto: ", round(pred_test[worst_idx],0), "g \n",
    "Differenza di peso nella predizione: ", round(errors[worst_idx],0), "g")

##  Peso predetto:  3063 g 
##  Differenza di peso nella predizione:  917 g

3. Visualizzazioni Grafiche

Mostriamo graficamente come varia l’andamento del Peso dei neonati con la Gestazione tra madri fumatrice e non fumatrici

ggplot(dati, aes(x = Gestazione, y = Peso, color = factor(Fumatrici))) +
  geom_jitter(width = 0.6, height = 0, alpha = 0.5, size = 1) +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  labs(x = "Settimane di gestazione",
       y = "Peso alla nascita (g)",
       color = "Fumo materno",
       title = "Peso alla nascita vs Settimane di gestazione",
       subtitle = "Confronto tra figli di madri fumatrici e non") +
  scale_color_manual(values = c("0" = "steelblue", "1" = "tomato"),
                     labels = c("Non fumatrice", "Fumatrice")) +
  theme_minimal()

come visto la classe Fumatrice è sotto rappresentata rispetto alla classe Non Fummatrice, dalle linee di regressione si osserva che i neonati nati da madri fumatrici verranno predetti con un peso minore rispetto al caso di madri non fumatrici.

Si vuole mettere a confronto la distribuzione del peso dei neonati nati tra parto cesare e naturale nei tre ospedali

ggplot(dati, aes(x = Tipo.parto, y = Peso, fill = Tipo.parto)) +
  geom_boxplot() +
  facet_wrap(~ Ospedale) +
  labs(x = "Tipo di parto",
       y = "Peso alla nascita (g)",
       fill = "Tipo di parto",
       title = "Peso alla nascita per tipo di parto e per ospedale") +
  theme_minimal()

si puo concludere che il tipo di parto non influenza significativamente il peso alla nascita, anche confrontando tra ospedali, le distribuzioni rimangono molto simili, ciò suggerisce che eventuali differenze tra ospedali sono minime o trascurabili rispetto alla variabilità individuale.

Visualizziamo la mappa delle correlazioni per le variabili continue

numeriche <- dati[, c("Gestazione", "Peso", "Lunghezza", "Cranio", "Anni.madre")]

cor_mat <- cor(numeriche, use = "complete.obs")

cor_df <- melt(cor_mat)

ggplot(cor_df, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(value, 2)), size = 3) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
                       limits = c(-1, 1), name = "Correlazione") +
  theme_minimal() +
  labs(title = "Heatmap delle correlazioni tra variabili continue")

già abbiamo fatto considerazioni sulle correlazioni.

Facciamo uno scatterplot che metta in relazione il peso osservato e il peso predetto dal modello mod_lambda

dati_pred <- data.frame(
  PesoOsservato = test_data$Peso,
  PesoPredetto = (lambda_opt*predict(mod_lambda, newdata = test_data)+1)^(1/lambda_opt))

dati_pred <- dati_pred %>%
  mutate(scarto = abs(PesoPredetto - PesoOsservato))

max_scarto <- dati_pred %>% filter(scarto == max(scarto))

ggplot(dati_pred, aes(x = PesoOsservato, y = PesoPredetto)) +
  geom_point(alpha = 0.6, size = 1.5, color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "brown", linetype = "dashed") +
  geom_segment(aes(xend = PesoOsservato, yend = PesoOsservato),
               color = "gray70", alpha = 0.5) +

  geom_point(data = max_scarto, aes(x = PesoOsservato, y = PesoPredetto),
             color = "red", size = 3) +
  
  labs(x = "Peso osservato (in g)",
       y = "Peso predetto (in g)",
       title = "Peso predetto vs Peso osservato") +
  theme_minimal()

il punto di massimo scarto sembra essere quello trovato in precedenza.

Infine, Facciamo uno scatterplot dove si mette in relazione Lunghezza e Peso facendo distinzione tra neoanti nati da madri fumatrici e non

ggplot(dati, aes(x = Lunghezza, y = Cranio)) +
  geom_point(aes(fill = factor(Fumatrici), size = Peso),
             shape = 21, color = "black", alpha = 0.2) +
  geom_smooth(aes(color = factor(Fumatrici)), method = "lm", formula = y ~ x, se = FALSE) +
  scale_fill_manual(values = c("0" = "steelblue", "1" = "tomato")) +
  scale_color_manual(values = c("0" = "steelblue", "1" = "tomato")) +
  labs(x = "Lunghezza", y = "Diametro craniale",
       fill = "Fumo", color = "Fumo", size = "Peso") +
  theme_minimal()

il rapporto tra Lunghezza e Cranio non cambia significativamente a seconda del fumo materno, il fumo materno potrebbe non alterare direttamente il rapporto lunghezza-diametro craniale.