1. Import del dataset

#import dei pacchetti
library(knitr)
library(moments)
library(kableExtra)
library(TeachingDemos)
library(dplyr)
library(Metrics)
library(car)
library(MASS)
library(ggplot2)
library(patchwork)
library(rgl)

df <- read.csv("neonati.csv")
df$Fumatrici <- factor(df$Fumatrici)
df$Tipo.parto <- factor(df$Tipo.parto)
df$Ospedale <- factor(df$Ospedale)
df$Sesso <- factor(df$Sesso)
attach(df)
n <- nrow(df)
print(paste("Il dataset ha",n, "righe"))

[1] "Il dataset ha 2500 righe"

2. Analisi e Modellizzazione

Analisi preliminare

Analisi descrittiva delle variabili

Sono presenti 10 variabili.

Ci sono 4 variabili qualitative:

Fumatrici: qualitativa nominale (dummy)
Tipo.parto: qualitativa nominale
Ospedale: qualitativa nominale
Sesso: qualitativa nominale

6 variabili quantitative:

Peso: la variabile risposta, quantitativa continua
Anni.madre: quantitativa discreta
N. gravidanze: quantitativa discreta
Gestazione: quantitativa continua
Lunghezza: quantitativa continua
Cranio: quantitativa continua

Statistiche descrittive delle variabili quantitative

selected_cols <- df[, c("Anni.madre","N.gravidanze","Gestazione", "Peso", "Lunghezza", "Cranio")]

stats <- data.frame(
  Mean           = apply(selected_cols, 2, mean),
  Median         = apply(selected_cols, 2, median),
  Min          = apply(selected_cols, 2, min),
  Max         = apply(selected_cols, 2, max),
  IQR             = apply(selected_cols, 2, IQR),
  Variance        = apply(selected_cols, 2, var),
  Std_dev   = apply(selected_cols, 2, sd),
  Skewness  = apply(selected_cols, 2, skewness),
  Kurtosis  = apply(selected_cols, 2, kurtosis)
  
)
stats$Kurtosis = stats$Kurtosis - 3

stats <- t(stats)
stats <- round(stats, 2)
colnames(stats) <- colnames(selected_cols)
kable(stats) %>%
  kable_styling(full_width = FALSE, position = 'left')

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Mean	28.16	0.98	38.98	3284.08	494.69	340.03
Median	28.00	1.00	39.00	3300.00	500.00	340.00
Min	0.00	0.00	25.00	830.00	310.00	235.00
Max	46.00	12.00	43.00	4930.00	565.00	390.00
IQR	7.00	1.00	2.00	630.00	30.00	20.00
Variance	27.81	1.64	3.49	275665.68	692.67	269.79
Std_dev	5.27	1.28	1.87	525.04	26.32	16.43
Skewness	0.04	2.51	-2.07	-0.65	-1.51	-0.79
Kurtosis	0.38	10.99	8.26	2.03	6.49	2.95

Dai parametri descrittivi delle variaibili quantitative si può osservare che per tutte le variabili i valori di media e mediana sono abbastanza vicini tra loro, indicando un’assenza di outlier significativi. Tuttavia, si osserva che la variabile Anni.madre ha un valore minimo di 0, il che porta a pensare che ci sia qualche errore di imputazione su alcuni dati. Un altro valore estremo sembra quello del massimo di Numero di gravidanze che risulta pari a 12. Per quanto riguarda la skewness, si osserva che la variabile degli anni è praticamente simmetrica, così come il peso; il numero di gravidanze ha un’asimmetria positiva mentre le restanti variabili hanno un’asimmetria negativa. Tutte le variabili risultano invece leptocurtiche ad eccezione della variabile anni.madre, che è praticamente mesocurtica.

Curve di densità di alcune variabili quantitative

par(mfrow = c(2, 2))

col_names = c("Anni.madre", "Peso", "Lunghezza", "Cranio")

for (col in col_names) {
  plot(density(df[[col]]), main = paste("Density function of", col))
  
  test_result <- shapiro.test(df[[col]])
  
  p_value <- format(test_result$p.value, scientific = TRUE, digits = 3)
  
  text(x = min(df[[col]], na.rm = TRUE),
       y = max(density(df[[col]], na.rm = TRUE)$y)*0.9, 
       labels = paste("Shapiro test\np-value:", p_value), 
       pos = 4, col = "blue", cex = 0.8)
}

Osservando le distribuzioni delle densità di probabilità di queste variabili si conferma quanto detto con la skewness e con la kurtosi, infatti si osservano le code e la leptocurticità. Inoltre lo Shapiro-test porta a rifiutare l’ipotesi che le distribuzioni siano normali.

Barplot di alcune variabili quantitative

par(mfrow = c(1, 3))

discr_col = c("N.gravidanze", "Gestazione", "Anni.madre")

for (col in discr_col) {
  
  x <- df[[col]]
  
  h <- barplot(table(x), main = paste("Barplot of", col), col = "lightblue")
  
  test_result <- shapiro.test(df[[col]])
  
  p_value <- format(test_result$p.value, scientific = TRUE, digits = 3)
  
  usr <- par("usr")
  text(x = usr[2], y = usr[4] * 0.9,
        labels = paste("Shapiro test\np-value:", p_value),
        pos = 2, col = "blue", cex = 1)
}

Si può osservare la simmetria della variabile anni, e il fatto che sono presenti dei dati corrispondenti ad anni inferiori a 13, ma sono in numero esiguo e non sembrano influenzare la forma della distribuzione.

Si osserva inoltre come il numero di gravidanze abbia un conteggio decrescente, ovvero la maggior parte delle donne ereno alla prima o seconda gravidanza alla nascita del bambino.

Per quanto riguarda la variabile gestazione, si osserva che la distribuzione è centrata sulle 40 settimane e ha una coda negativa ad indicare che come outlier le nascite premature avvengono maggiormente rispetto a nascite tardive.

Boxplot delle variabili quantitative

Con i seguenti grafici si mostrano i boxplot delle variabili quantitative:

Si conferma la presenza dei valori anomali della variabile anni
Si conferma la forma delle distribuzioni di numero di gravidanze e gestazione, rispettivamente con coda positiva e negativa
Si osserva che le variabili peso, lunghezza e cranio sono abbastanza simmetriche ma hanno una coda negativa in quanto i valori dell’estremità inferiore della distribuzione principale sono maggiormente dispersi rispetto a quelli superiori.

par(mfrow = c(2, 3))

num_cols <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio" )

for (col in num_cols) {
  
  x <- na.omit(df[[col]]) 
  
  boxplot(x, main = paste("Boxplot of", col), col = "lightblue")

}

Si analizzano gli outlier della variabile anni.madre

kable(df[df$Anni.madre <= "14", ])%>%
  kable_styling(full_width = FALSE, position = 'left')

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
138	13	0	38	2760	470	325	Nat	osp2	F
1075	14	1	39	3510	490	365	Nat	osp2	M
1152	1	1	41	3250	490	350	Nat	osp2	F
1380	0	0	39	3060	490	330	Nat	osp3	M
1532	14	0	39	3550	500	355	Ces	osp1	M

Per quanto poco probabili gli anni 13 e 14 potrebbero essere corretti, mentre gli anni 1 e 0 sono sicuramente errori di imputazione. Si decide di sostituire questi valori con la media degli anni delle madri presenti nel dataset.

excluded_val <- c(1152,1380) 
mean_age <- mean(df$Anni.madre[-excluded_val])

Anni.madre[excluded_val] <- mean_age

Si riosserva il boxplot Anni.madre:

boxplot(Anni.madre, main = paste("Boxplot of Anni.madre"), col="lightblue")

Distribuzioni di frequenza delle variabili qualitative

#function for tables creation
table_freq <- function(variable, n){
  
  frq_ass <- table(variable)
  frq_rel <- frq_ass/n
  
  distr_freq <- cbind(frq_ass, frq_rel)
  
  return(distr_freq)
  
  
}

qual_cols <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

tables <- list()

for (col in qual_cols) {
  
  tables[[col]] = table_freq(df[[col]], n)
  
}


kable(tables$Fumatrici, caption = "Fumatrici")%>%
  kable_styling(full_width = FALSE, position = 'left')

Fumatrici
	frq_ass	frq_rel
0	2396	0.9584
1	104	0.0416

kable(tables$Tipo.parto, caption = "Tipo parto")%>%
  kable_styling(full_width = FALSE, position = 'left')

Tipo parto
	frq_ass	frq_rel
Ces	728	0.2912
Nat	1772	0.7088

kable(tables$Ospedale, caption = "Ospedale")%>%
  kable_styling(full_width = FALSE, position = 'left')

Ospedale
	frq_ass	frq_rel
osp1	816	0.3264
osp2	849	0.3396
osp3	835	0.3340

kable(tables$Sesso, caption="Sesso")%>%
  kable_styling(full_width = FALSE, position = 'left')

Sesso
	frq_ass	frq_rel
F	1256	0.5024
M	1244	0.4976

Dalle tabelle di frequenza si può osservare che i dati sono ben distribuiti nei tre ospedali e anche tra i due sessi. Si osserva invece una preponderanza di parti naturali rispetto ai cesarei (circa il 70% contro il 30%) e una forte maggioranza (96%) di donne non fumatrici rispetto alle fumatrici.

Test di ipotesi: in alcuni ospedali si fanno più parti cesarei

H0: Non c’è differenza nel numero di cesarei effettuati nei tre ospedali

H1: C’è differenza significativa nel numero di cesarei effettuati nei tre ospedali

# Tabella di contingenza
freq_parti <- table(Ospedale, Tipo.parto)
kable(freq_parti) %>%
  kable_styling(full_width = FALSE, position = 'left')

	Ces	Nat
osp1	242	574
osp2	254	595
osp3	232	603

# Test del Chi-quadro
chi <- chisq.test(freq_parti)
print(chi)


    Pearson's Chi-squared test

data:  freq_parti
X-squared = 1.0972, df = 2, p-value = 0.5778

Il test del chi quadro mostra se la proporzione di cesarei rispetto ai naturali varia da un ospedale all’altro. Il p-value superiore a 0.05 indica l’accettazione dell’ipotesi nulla, per cui non c’è differenza significativa tra le proporzioni dei cesarei nei vari ospedali.

Test d’ipotesi: La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

H0: Non c’è differenza tra la media della variabile considerata del campione in analisi e quella della popolazione

H1: C’è una differenza significativa tra la media della variabile considerata del campione in analisi e quella della popolazione

#population values
pop_weight = 3300
pop_lungh = 500

#t-tests
t_peso <- t.test(Peso, mu = pop_weight, conf.level = 0.95)
t_lungh <- t.test(Lunghezza, mu = pop_lungh, conf.level = 0.95)

t_test_results <-
  
data.frame(
  Variabile = c("Peso", "Lunghezza"),
  media_pop = round(c(pop_weight, pop_lungh), 2),
  media_obs = round(c(t_peso$estimate,t_lungh$estimate), 2),
  t_obs = round(c(t_peso$statistic, t_lungh$statistic),2),
  p_value = c(round(t_peso$p.value, 3), format(t_lungh$p.value, scientific = T, digits = 3))
)

kable(t_test_results)%>%
  kable_styling(full_width = FALSE, position = 'left')

Variabile	media_pop	media_obs	t_obs	p_value
Peso	3300	3284.08	-1.52	0.13
Lunghezza	500	494.69	-10.08	1.81e-23

Il t-test sul peso neonatale mostra un p-value maggiore di 0.05, inoltre l’intervallo di confidenza al 95% comprende la media della popolazione. Pertanto, si ritiene accettata l’ipotesi nulla e quindi non significativa la differenza della media del campione con quella della popolazione.

Nel caso della lunghezza, invece, il t-test ha dato un valore molto prossimo a 0 e comunque <<0.05, inoltre l’intervallo di confidenza non contiene la media della popolazione. Pertanto, si rifiuta l’ipotesi nulla e risulta esserci una differenza statisticamente significativa tra la media della lunghezza neonatale del campione e quella della popolazione. In particolare, in media i valori di lunghezza neonatale sono inferiori alla media della popolazione. Infatti, osservando sia il boxplot sia la distribuzione della variabile, si può osservare che presenta molti outlier nei valori inferiori alla media ed ha un’asimmetria negativa.

Nei grafici sottostanti si mostrano le distribuzioni t di student con i limiti al 95% di confidenza e il punto in cui ricade la statistica di test.

par(mfrow = c(1, 2))

dof <- 2499
t_obs <- -1.516
alpha <- 0.05

t_val <- seq(-4, +4, length = 1000)

t_critico_dx <- qt(1 - alpha/2, dof)
t_critico_sx <- qt(alpha/2, dof)

plot(t_val, dt(t_val, dof), type = "l", lwd = 2, 
     ylab = "Densità", xlab = "t", 
     main = "Distribuzione t-student peso neonatale")

abline(v = t_critico_dx, col = "red", lwd = 2, lty = 2)
abline(v = t_critico_sx, col = "red", lwd = 2, lty = 2)

points(t_obs, dt(t_obs, dof), col = "blue", pch = 16, cex = 2)

legend("topright", legend = c("alpha = 0.05", "t osservato"), 
       col = c("red", "blue"), lty = c(2,NA),lwd = c(2, NA), pch = c(NA, 16), pt.cex = c(NA,2))


dof <- 2499
t_obs <- -10.84
alpha <- 0.05

t_val <- seq(-15, +15, length = 1000)

t_critico_dx <- qt(1 - alpha/2, dof)
t_critico_sx <- qt(alpha/2, dof)

plot(t_val, dt(t_val, dof), type = "l", lwd = 2, 
     ylab = "Densità", xlab = "t", 
     main = "Distribuzione t-student lunghezza neonatale")

abline(v = t_critico_dx, col = "red", lwd = 2, lty = 2)
abline(v = t_critico_sx, col = "red", lwd = 2, lty = 2)

points(t_obs, dt(t_obs, dof), col = "blue", pch = 16, cex = 2)

legend("topright", legend = c("alpha = 0.05", "t osservato"), 
       col = c("red", "blue"), lty = c(2,NA),lwd = c(2, NA), pch = c(NA, 16), pt.cex = c(NA,2))

Test d’ipotesi: Le misure antropometriche sono significativamente diverse tra i due sessi

H0: Non c’è differenza tra le misure antropometriche delle femmine rispetto ai maschi

H1: c’è una differenza significativa tra le misure antropometriche delle femmine rispetto ai maschi

Si visualizzano i boxplot delle variabili

par(mfrow = c(1, 3))

boxplot(Peso~Sesso,main = "Peso vs Sesso", col=c("lightblue", "lightgreen"))
boxplot(Lunghezza~Sesso, main = "Lunghezza vs Sesso", col=c("lightblue", "lightgreen"))
boxplot(Cranio~Sesso, main = "Cranio vs Sesso", col=c("lightblue", "lightgreen"))

Si effettuano i test e si riportano in tabella

t_peso <- t.test(Peso ~ Sesso, data = df)
t_lungh <- t.test(Lunghezza ~ Sesso, data = df)
t_cranio <- t.test(Cranio ~ Sesso, data = df)


t_test_results <-
  
data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  media_M = round(c(mean(Peso[Sesso == "M"]), mean(Lunghezza[Sesso == "M"]), mean(Cranio[Sesso == "M"])), 2),
  media_F = round(c(mean(Peso[Sesso == "F"]), mean(Lunghezza[Sesso == "F"]), mean(Cranio[Sesso == "F"])), 2),
  t_obs = round(c(t_peso$statistic, t_lungh$statistic, t_cranio$statistic),2),
  p_value = c(format(t_peso$p.value, scientific = T, digits = 3), format(t_lungh$p.value, scientific = T, digits = 3),format(t_cranio$p.value, scientific = T, digits = 3))
)


kable(t_test_results)%>%
  kable_styling(full_width = FALSE, position = 'left')

Variabile	media_M	media_F	t_obs	p_value
Peso	3408.22	3161.13	-12.11	8.03e-33
Lunghezza	499.67	489.76	-9.58	2.24e-21
Cranio	342.45	337.63	-7.41	1.72e-13

Per tutte le caratteristiche antropometriche il t-test ha mostrato un p-value inferiore a 0.05, tale per cui si rifiutano le ipotesi nulle ed emerge una differenza statisticamente significativa tra la categoria dei neonati maschi e dei neonati femmine. I neonati maschi risultano avere in media un peso, una lunghezza e le dimensioni del cranio maggiori rispetto a quelle delle neonate femmine.

Creazione del Modello di Regressione

Correlazioni tra le variabili

Si osservano le correlazioni tra le variabili del dataset, escludendo quelle qualitative

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr=usr))
  par(usr = c(0, 1, 0, 1))
  r <- (cor(x, y))
  txt <- format(c(r, 1), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = 1.5)
}

df_num <- dplyr::select(df, where(is.numeric))

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

Dalla matrice di correlazione si può osservare una forte correlazione positiva tra peso, lunghezza, cranio e periodo di gestazione. Tra gestazione e peso si osserva una leggera tendenza ad una relazione quadratica, così come tra peso e lunghezza.

Per le variabili qualitative si osservano i boxplot.

par(mfrow=c(2,2))

boxplot(Peso~Sesso, col=c("lightblue", "lightgreen"))
boxplot(Peso~Fumatrici, col=c("lightblue", "lightgreen"))
boxplot(Peso~Tipo.parto, col=c("lightblue", "lightgreen"))
boxplot(Peso~Ospedale, col=c("lightblue", "lightgreen", "orange"))

Per il sesso è stato già osservato precedentemente che la differenza tra le medie è signficiativa. Si effettua un t-test allo stesso modo relativamente alle altre variabili.

Variabili fumatrici e tipo.parto

H0: Non c’è differenza nel peso medio del neonato in base alla categoria della variabile (fumatrice o non, oppure parto naturale o cesareo)

H1: C’è una differenza significativa tra il peso di un neonato relativamente alla categoria della variabile

t_fum <- t.test(Peso ~ Fumatrici, data = df)
t_parto <- t.test(Peso ~ Tipo.parto, data = df)

fum_parto_res <-
  data.frame(
  Variabile = c("Fumatrici", "Tipo.parto"),
  t_obs = round(c(t_fum$statistic, t_parto$statistic),2),
  p_value = round(c(t_fum$p.value, t_parto$p.value), 3)
)


kable(fum_parto_res)%>%
  kable_styling(full_width = FALSE, position = 'left')

Variabile	t_obs	p_value
Fumatrici	1.03	0.303
Tipo.parto	-0.13	0.897

Il test mostra un p-value superiore a 0.05 per entrambe le variabili, pertanto si accettano le ipotesi nulle e non c’è differenza significativa tra il peso di un neonato figlio di una fumatrice e quello di una non fumatrice, così come tra quello di un neonato nato da un cesareo o da un parto naturale. Pertanto queste due variabili potrebbero non essere rilevanti per il modello di regressione.

Variabile Ospedale

pairwise.t.test(Peso, Ospedale, paired = F, pool.sd = T, p.adjust.method = "bonferroni")


    Pairwise comparisons using t tests with pooled SD 

data:  Peso and Ospedale 

     osp1 osp2
osp2 1.00 -   
osp3 0.33 0.33

P value adjustment method: bonferroni

I p-value sono tutti superiori a 0.05 quindi non c’è differenza significativa tra il peso dei neonati nei tre ospedali. Pertanto anche in questo caso non ci si aspetta che la variabile abbia un peso rilevante nel modello di regressione.

Modello di Regressione

Modello di regressione con tutte le variabili:

mod1<-lm(df$Peso~ . , data = df)

summary(mod1)


Call:
lm(formula = df$Peso ~ ., 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

Descrizione del peso delle variabili in questo modello:

Anni.madre: variabile non significativa, p-value > 0.05
N.gravidanze: variabile significativa ad un livello di confidenza del 95%, ogni gravidanza in più aumenta il peso di 11 grammi circa con un p-value < 0.05.
Gestazione: ha una elevata significatività. Ogni settimana di gestazione in più in media aumenta il peso di 32 g
Lunghezza e Cranio: elevata significatività. Ogni unità in più aumenta il peso di circa 10 g in media.
Fumatrici: rispetto alle non fumatrici hanno una correlazione negativa, ma la significatività non risulta molto elevata, è superiore allo 0.05.
Il peso dei maschi, come visto, è significativamente maggiore, in media di 77g, di quello delle femmine
Tipo.parto: Significativa la differenza tra naturale e cesareo (p-value <0.05): un bambino nato da un parto naturale pesa in media 30g in più di un bambino nato con cesareo.
Differenza tra ospedale2 e 1 non significativa
Differenza tra ospedale3 e 1 leggermente significativa con un peso maggiore di circa 28g in media.

Le variabili più significative sono dunque: Gestazione, Lunghezza e Cranio. Inoltre è significativa la differenza tra maschi e femmine. L’R2 aggiustato di questo modello è circa 0.73, ad indicare che è in grado di spiegare il 73% della variabilità del peso neonatale.

Selezione del modello ottimale

Si provano ad aggiungere gli effetti quadratici di gestazione, lunghezza e cranio

mod2 <- update(mod1,~.+I(Gestazione^2)+I(Lunghezza^2)+I(Cranio^2))
summary(mod2)


Call:
lm(formula = df$Peso ~ Anni.madre + N.gravidanze + Fumatrici + 
    Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + 
    Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2), data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1157.69  -180.60    -8.68   156.80  1441.87 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.360e+03  1.159e+03  -1.173  0.24085    
Anni.madre       4.493e-01  1.108e+00   0.405  0.68529    
N.gravidanze     1.377e+01  4.562e+00   3.017  0.00257 ** 
Fumatrici1      -2.377e+01  2.693e+01  -0.883  0.37752    
Gestazione       3.848e+02  6.737e+01   5.712 1.25e-08 ***
Lunghezza       -2.919e+01  4.428e+00  -6.593 5.24e-11 ***
Cranio          -5.465e+00  9.787e+00  -0.558  0.57667    
Tipo.partoNat    2.712e+01  1.182e+01   2.295  0.02180 *  
Ospedaleosp2    -1.114e+01  1.314e+01  -0.848  0.39670    
Ospedaleosp3     3.085e+01  1.320e+01   2.338  0.01949 *  
SessoM           7.198e+01  1.097e+01   6.561 6.49e-11 ***
I(Gestazione^2) -4.498e+00  8.836e-01  -5.091 3.82e-07 ***
I(Lunghezza^2)   4.074e-02  4.536e-03   8.982  < 2e-16 ***
I(Cranio^2)      2.338e-02  1.443e-02   1.620  0.10527    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 267.7 on 2486 degrees of freedom
Multiple R-squared:  0.7413,    Adjusted R-squared:   0.74 
F-statistic: 548.1 on 13 and 2486 DF,  p-value: < 2.2e-16

Lunghezza e gestazione al quadrato risultano molto significative, mentre cranio no. Inoltre perde significatività anche la forma lineare della variabile cranio. R2 aumenta da 0.72 a 0.74

Si decide di rimuovere la variabile cranio al quadrato

mod3 <- update(mod2,~.-I(Cranio^2))
summary(mod3)


Call:
lm(formula = df$Peso ~ Anni.madre + N.gravidanze + Fumatrici + 
    Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + 
    Sesso + I(Gestazione^2) + I(Lunghezza^2), data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1162.18  -178.48    -8.91   156.48  1374.43 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -2.535e+03  9.043e+02  -2.803  0.00510 ** 
Anni.madre       4.853e-01  1.109e+00   0.438  0.66160    
N.gravidanze     1.379e+01  4.563e+00   3.023  0.00253 ** 
Fumatrici1      -2.420e+01  2.693e+01  -0.899  0.36892    
Gestazione       3.447e+02  6.267e+01   5.500 4.18e-08 ***
Lunghezza       -3.216e+01  4.032e+00  -7.976 2.29e-15 ***
Cranio           1.038e+01  4.189e-01  24.782  < 2e-16 ***
Tipo.partoNat    2.740e+01  1.182e+01   2.319  0.02048 *  
Ospedaleosp2    -1.129e+01  1.314e+01  -0.859  0.39047    
Ospedaleosp3     3.020e+01  1.320e+01   2.288  0.02220 *  
SessoM           7.228e+01  1.097e+01   6.587 5.45e-11 ***
I(Gestazione^2) -3.982e+00  8.243e-01  -4.831 1.44e-06 ***
I(Lunghezza^2)   4.376e-02  4.135e-03  10.582  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 267.8 on 2487 degrees of freedom
Multiple R-squared:  0.7411,    Adjusted R-squared:  0.7398 
F-statistic: 593.1 on 12 and 2487 DF,  p-value: < 2.2e-16

In questo modo si perde leggermente in R2, che però rimane circa 0.74, ma torna la significatività della variabile cranio al primo grado.

Si rimuovono le variabili non significative

mod4 <- update(mod3,~.-Anni.madre-Ospedale-Fumatrici)
summary(mod4)


Call:
lm(formula = df$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2), 
    data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1172.27  -178.53   -11.31   163.11  1405.90 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -2.426e+03  9.049e+02  -2.681 0.007384 ** 
N.gravidanze     1.474e+01  4.245e+00   3.472 0.000525 ***
Gestazione       3.373e+02  6.264e+01   5.384 7.98e-08 ***
Lunghezza       -3.199e+01  4.034e+00  -7.931 3.25e-15 ***
Cranio           1.041e+01  4.190e-01  24.851  < 2e-16 ***
Tipo.partoNat    2.799e+01  1.183e+01   2.365 0.018088 *  
SessoM           7.260e+01  1.099e+01   6.607 4.77e-11 ***
I(Gestazione^2) -3.885e+00  8.237e-01  -4.716 2.54e-06 ***
I(Lunghezza^2)   4.358e-02  4.137e-03  10.535  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 268.2 on 2491 degrees of freedom
Multiple R-squared:  0.7398,    Adjusted R-squared:  0.739 
F-statistic: 885.4 on 8 and 2491 DF,  p-value: < 2.2e-16

Ora nel modello vi sono tutte variabili significative ed è in grado di spiegare il 74% della variabilità del peso neonatale (R2 aggiustato).

Si analizzano le interazioni tra diverse variabili nel modello:

mod5 <- update(mod4,~.+Lunghezza*Cranio+Anni.madre*Gestazione)
summary(mod5)


Call:
lm(formula = df$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2) + 
    Anni.madre + Lunghezza:Cranio + Gestazione:Anni.madre, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1165.36  -181.67   -12.13   162.44  1333.87 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -1.777e+03  1.239e+03  -1.434 0.151819    
N.gravidanze           1.390e+01  4.556e+00   3.051 0.002304 ** 
Gestazione             2.271e+02  7.115e+01   3.193 0.001428 ** 
Lunghezza             -3.011e+01  4.091e+00  -7.360 2.48e-13 ***
Cranio                 2.085e+01  5.081e+00   4.102 4.22e-05 ***
Tipo.partoNat          2.724e+01  1.182e+01   2.304 0.021309 *  
SessoM                 7.325e+01  1.098e+01   6.674 3.06e-11 ***
I(Gestazione^2)       -2.965e+00  8.913e-01  -3.327 0.000891 ***
I(Lunghezza^2)         4.907e-02  5.063e-03   9.692  < 2e-16 ***
Anni.madre            -5.240e+01  2.046e+01  -2.561 0.010507 *  
Lunghezza:Cranio      -2.139e-02  1.039e-02  -2.059 0.039597 *  
Gestazione:Anni.madre  1.365e+00  5.264e-01   2.593 0.009578 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 267.8 on 2488 degrees of freedom
Multiple R-squared:  0.741, Adjusted R-squared:  0.7398 
F-statistic:   647 on 11 and 2488 DF,  p-value: < 2.2e-16

L’interazione tra cranio e lunghezza e quella tra gestazione e anni.madre risultano significative, tuttavia l’R2 rimane praticamente invariato, pertanto si sceglie di non includere queste variabili che aumenterebbero la complessità del modello e di scegliere il mod4.

Per completezza, si considera anche il modello con solo le variabili lineari e senza quelle non significative:

mod6 <- update(mod1,~.-Anni.madre-Ospedale-Fumatrici)
summary(mod6)


Call:
lm(formula = df$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1129.31  -181.70   -16.31   161.07  2638.85 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6707.2971   135.9911 -49.322  < 2e-16 ***
N.gravidanze     12.7558     4.3366   2.941   0.0033 ** 
Gestazione       32.2713     3.7941   8.506  < 2e-16 ***
Lunghezza        10.2864     0.3007  34.207  < 2e-16 ***
Cranio           10.5057     0.4260  24.659  < 2e-16 ***
Tipo.partoNat    30.0342    12.0969   2.483   0.0131 *  
SessoM           77.9285    11.1905   6.964 4.22e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.3 on 2493 degrees of freedom
Multiple R-squared:  0.7277,    Adjusted R-squared:  0.727 
F-statistic:  1110 on 6 and 2493 DF,  p-value: < 2.2e-16

Si verifica la scelta del modello con il criterio BIC:

bic <- BIC(mod1, mod2, mod3, mod4, mod5, mod6)

kable(bic)%>%
  kable_styling(full_width = FALSE, position = 'left')

	df	BIC
mod1	12	35241.84
mod2	15	35147.92
mod3	14	35142.74
mod4	10	35123.36
mod5	13	35135.79
mod6	8	35221.75

Il modello con i BIC più basso è quello selezionato precedentemente ovvero il modello 4.

Si effettua un test dell’ANOVA tra il modello completamente lineare (mod6) e quello con i termini quadratici scelto (mod4)

anova(mod4, mod6)

Analysis of Variance Table

Model 1: df$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
    Sesso + I(Gestazione^2) + I(Lunghezza^2)
Model 2: df$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
    Sesso
  Res.Df       RSS Df Sum of Sq     F    Pr(>F)    
1   2491 179236346                                 
2   2493 187601677 -2  -8365330 58.13 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il p-value risulta molto prossimo a 0, ad indicare che l’aggiunta dei termini quadratici è molto significativa per la varianza spiegata dal modello.

Si calcolano anche gli rmse dei modelli e si riportano le informazioni nella seguente tabella

y_oss <- df$Peso

rmse <- c()
r2_adj <- c()

for (el in list(mod1, mod2, mod3, mod4, mod5, mod6)){
  rmse <- c(rmse, rmse(y_oss, fitted(el)))
  r2_adj <- c(r2_adj, summary(el)$adj.r.squared)
}

tabella <-
data.frame(
  "Modello" = c("mod1", "mod2", "mod3", "mod4", "mod5", "mod6"),
  "R2_adj" = r2_adj,
  "RMSE" = rmse,
  "BIC" = bic$BIC
)

kable(tabella) %>%
  kable_styling(full_width = FALSE, position = 'left')

Modello	R2_adj	RMSE	BIC
mod1	0.7278038	273.3222	35241.84
mod2	0.7399761	266.9799	35147.92
mod3	0.7398061	267.1208	35142.74
mod4	0.7389825	267.7584	35123.36
mod5	0.7398197	267.1676	35135.79
mod6	0.7270194	273.9355	35221.75

Il modello con RMSE più basso sarebbe il 2, tuttavia esso considera anche variabili non significative che ne aumentano la complessità. Inoltre, come nel caso dell’R2, l’aumento del valore di RMSE, in particolare nel modello 4, è di piccola entità. Si sceglie dunque il modello 4 che è anche quello con BIC minore.

Si verifica la presenza di multicollinearità nel modello:

vif_4 <- vif(mod4)

vif_df <- data.frame(VIF = vif_4)


kable(vif_df) %>%
  kable_styling(full_width = FALSE, position = 'left')

	VIF
N.gravidanze	1.026094
Gestazione	475.891404
Lunghezza	391.454121
Cranio	1.645147
Tipo.parto	1.003897
Sesso	1.048541
I(Gestazione^2)	453.552580
I(Lunghezza^2)	372.926478

I VIF mostrano multicollinearità tra le variabili e i loro quadrati, come è normale che avvenga. Non c’è però multicollinearità tra le variabili di primo grado, infatti, se si calcolano i VIF per il primo modello, ovvero quello che non contiene i termini quadratici, si può osservare come tutte le variabili abbiano valori inferiori a 5.

vif_1 <- round(vif(mod1), 3)

kable(vif_1)%>%
  kable_styling(full_width = FALSE, position = 'left')

	GVIF	Df	GVIF^(1/(2*Df))
Anni.madre	1.187	1	1.090
N.gravidanze	1.186	1	1.089
Fumatrici	1.007	1	1.004
Gestazione	1.696	1	1.302
Lunghezza	2.086	1	1.444
Cranio	1.631	1	1.277
Tipo.parto	1.004	1	1.002
Ospedale	1.004	2	1.001
Sesso	1.041	1	1.020

Analisi della qualità del modello

Analisi dei residui

par(mfrow=c(2,2))
plot(mod4)

I residui sono distribuiti in modo abbastanza randomico attorno alla media.

Il Q-Q plot invece mostra una deviazione dalla normalità a valori alti e a valori bassi

La varianza sembra rimanere costante

Dal grafico dei residui il valore 1551 sembra essere particolarmente influente e superare la distanza di cook.

Si effettuano i test statistici

#shapiro-wilk per la normalità
res_sh <- shapiro.test(residuals(mod4))

#breusch-pagan per l'omoschedasticità
res_bp <- lmtest::bptest(mod4) 

#darwin-watson èer la non correlazione
res_dw <- lmtest::dwtest(mod4) 


res_test_results <-
  
data.frame(
  Test = c("Shapiro-Wilk", "Breusch-Pagan", "Durbin-Watson"),
  Caratteristica_testata = c("Normalità", "Omoschedasticità", "Autocorrelazione"),
  p_value = c(format(res_sh$p.value, scientific = T, digits = 3), format(res_bp$p.value, scientific = T, digits = 3), round(res_dw$p.value, 3))
)

kable(res_test_results)%>%
  kable_styling(full_width = FALSE, position = 'left')

Test	Caratteristica_testata	p_value
Shapiro-Wilk	Normalità	8.11e-13
Breusch-Pagan	Omoschedasticità	5.79e-18
Durbin-Watson	Autocorrelazione	0.098

Shapiro-wilk p-value molto basso: si rifiuta l’ipotesi di normalità

Breusch-Pagan p-value molto basso: si rifiuta l’ipotesi di omoschedasticità

Darwin-Watson p-value leggermente sopra a 0.05: si accetta l’ipotesi di autocorrelazione dei residui.

Tutti e tre i test dei residui non risultano passati.

Si indagano i valori di leverage

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

Si evidenziano diversi valori leverage al di sopra della soglia.

Si indagano gli outliers

plot(rstudent(mod4))
abline(h=c(-2,2))

a <- outlierTest(mod4)
kable(data.frame("Bonf p-val" = format(a$bonf.p, scientific=T, digits = 3))) %>%
  kable_styling(full_width = FALSE, position = 'left')

	Bonf.p.val
1551	3.36e-06
1306	2.13e-03
155	2.17e-02
1399	2.90e-02
1694	2.94e-02

Ci sono diversi outliers, in particolare 5 hanno un p-value significativo.

Si indaga quali valori, tra leverage e outliers, superano la distanza di cook

cook<-cooks.distance(mod4)

plot(cook,ylim = c(0,2))

text(x = 500, y = 1.8,
    labels = paste("Max cook distance:", round(max(cook), 2)),
    , col = "blue", cex = 0.8)

Un valore risulta avere una distanza di cook molto alta, superiore a 1. Si identifica nel dataset la riga corrispondente.

kable(df[which.max(cook), ]) %>%
  kable_styling(full_width = FALSE, position = 'left')

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

Modello senza valore anomalo

Si rimuove la riga identificata che ha la distanza di cook superiore a 1 e si prova a rifare il modello di regressione

df_2 <- df[-1551, ]

mod1_n<-lm(df_2$Peso~ . , data = df_2)
mod2_n <- update(mod1_n,~.+I(Gestazione^2)+I(Lunghezza^2)+I(Cranio^2))
mod3_n <- update(mod2_n,~.-I(Cranio^2))
mod4_n <- update(mod3_n,~.-Anni.madre-Ospedale-Fumatrici)
mod5_n <- update(mod4_n,~.+Lunghezza*Cranio+Anni.madre*Gestazione)
mod6_n <- update(mod1_n,~.-Anni.madre-Ospedale-Fumatrici)

y_oss <- df_2$Peso

rmse_n <- c()
r2_adj_n <- c()

for (el in list(mod1_n, mod2_n, mod3_n, mod4_n, mod5_n, mod6_n)){
  rmse_n <- c(rmse_n, rmse(y_oss, fitted(el)))
  r2_adj_n <- c(r2_adj_n, summary(el)$adj.r.squared)
}

mod_no_outlier <-
data.frame(
"Modello" = c("mod1_n", "mod2_n","mod3_n","mod4_n","mod5_n","mod6_n"),
"R2_adj" = r2_adj_n,
"RMSE" = rmse_n,
"BIC" = BIC(mod1_n, mod2_n, mod3_n, mod4_n, mod5_n, mod6_n)$BIC
)

kable(mod_no_outlier) %>%
  kable_styling(full_width = FALSE, position = 'left')

Modello	R2_adj	RMSE	BIC
mod1_n	0.7378210	268.0693	35130.78
mod2_n	0.7437090	264.8823	35094.48
mod3_n	0.7430333	265.2845	35094.24
mod4_n	0.7423658	265.8425	35073.45
mod5_n	0.7434381	265.1288	35083.48
mod6_n	0.7372214	268.5913	35109.21

In questi modelli, in cui è stato rimosso l’outlier, valgono gli stessi ragionamenti dei precedenti per cui il modello migliore è quello con i termini quadratici e l’assenza delle variabili non significative (mod4_n). Il valore dell’R2 però, rispetto al mod4 aumenta e quello di RMSE diminuisce.

r2_4 <- r2_adj[4]
r2_4n <- r2_adj_n[4]
RMSE_4 <- rmse[4]
RMSE_4n <- rmse_n[4]

kable(data.frame("modello" = c("mod4", "mod4_n"),
           "R2_adj" = c(r2_4, r2_4n),
           "RSME" = c(RMSE_4, RMSE_4n))) %>%
  kable_styling(full_width = FALSE, position = 'left')

modello	R2_adj	RSME
mod4	0.7389825	267.7584
mod4_n	0.7423658	265.8425

Si riverificano i VIF per la multicollinearità e si confermano i valori alti per i termini quadratici.

vif_4 <- vif(mod4_n)

vif_df <- data.frame(VIF = vif_4)


kable(vif_df) %>%
  kable_styling(full_width = FALSE, position = 'left')

	VIF
N.gravidanze	1.026105
Gestazione	537.633834
Lunghezza	491.627198
Cranio	1.674221
Tipo.parto	1.003751
Sesso	1.048273
I(Gestazione^2)	509.008354
I(Lunghezza^2)	464.302191

Si calcola la distanza di cook per questo modello:

cook_n<-cooks.distance(mod4_n)

plot(cook_n, ylim = c(0,1))

text(x = 500, y = 0.8,
    labels = paste("Max cook distance:", round(max(cook_n), 2)),
    , col = "blue", cex = 0.8)

Questa volta nessun valore anomalo influenza le stime del modello di regressione.

Si effettuano nuovamente i test stastici sui residui e si confrontano con i risultati precedenti

res_sh_n <- shapiro.test(residuals(mod4_n))
res_bp_n <- lmtest::bptest(mod4_n) 
res_dw_n <- lmtest::dwtest(mod4_n)


res_test_results_n <-
  
data.frame(
  Test = c("Shapiro-Wilk", "Breusch-Pagan", "Durbin-Watson"),
  Caratteristica_testata = c("Normalità", "Omoschedasticità", "Autocorrelazione"),
  p_value_mod4n = c(format(res_sh_n$p.value, scientific = T, digits = 3), round(res_bp_n$p.value, 3), round(res_dw_n$p.value, 3)),
  p_value_mod4 = res_test_results$p_value
)

kable(res_test_results_n)%>%
  kable_styling(full_width = FALSE, position = 'left')

Test	Caratteristica_testata	p_value_mod4n	p_value_mod4
Shapiro-Wilk	Normalità	8.65e-12	8.11e-13
Breusch-Pagan	Omoschedasticità	0.031	5.79e-18
Durbin-Watson	Autocorrelazione	0.103	0.098

I test per i residui danno gli stessi risultati precedenti, anche se l’autocorrelazione e l’omoschedasticità risultano migliorati in termini di significatività, soprattutto l’omoschedasticità.

Trasformazione box-cox

Si prova ad effettuare una trasformazione sulla variabile y (peso) per vedere se l’analisi dei residui migliora

model_box  <- boxcox(mod4_n)

La lambda ottimale risulta essere 0.4646465.

lambda <- model_box$x[which.max(model_box$y)]
lambda

[1] 0.4646465

Si trasforma la variabile peso con la lambda trovata e si costruisce un nuovo modello di regressione, dipendente dalle stesse variabili selezionate precedentemente, con la variabile peso trasformata.

df_2$peso_trans <- (df_2$Peso^lambda - 1) / lambda

mod_bc <- lm(peso_trans ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2), data = df_2)

summary(mod_bc)


Call:
lm(formula = peso_trans ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2), 
    data = df_2)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.6298  -2.3001  -0.0706   2.1534  19.0174 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -9.344e+01  1.182e+01  -7.903 4.04e-15 ***
N.gravidanze     1.812e-01  5.527e-02   3.279  0.00106 ** 
Gestazione       4.183e+00  8.670e-01   4.825 1.49e-06 ***
Lunghezza        4.039e-02  5.942e-02   0.680  0.49669    
Cranio           1.334e-01  5.509e-03  24.211  < 2e-16 ***
Tipo.partoNat    3.509e-01  1.541e-01   2.278  0.02284 *  
SessoM           9.523e-01  1.431e-01   6.656 3.44e-11 ***
I(Gestazione^2) -4.782e-02  1.136e-02  -4.209 2.66e-05 ***
I(Lunghezza^2)   1.105e-04  6.051e-05   1.826  0.06799 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.493 on 2490 degrees of freedom
Multiple R-squared:  0.7709,    Adjusted R-squared:  0.7702 
F-statistic:  1047 on 8 and 2490 DF,  p-value: < 2.2e-16

Si osserva un aumento dell’R2 aggiustato, tuttavia si perde la significatività di alcune variabili come la lunghezza.

Inoltre, i test sui residui continuano a non passare:

res_sh_bc <- shapiro.test(residuals(mod_bc))
res_bp_bc <- lmtest::bptest(mod_bc) 
res_dw_bc <- lmtest::dwtest(mod_bc) 

res_test_results_bc <-
  
data.frame(
  Test = c("Shapiro-Wilk", "Breusch-Pagan", "Durbin-Watson"),
  Caratteristica_testata = c("Normalità", "Omoschedasticità", "Autocorrelazione"),
  p_value_modbc = c(format(res_sh_bc$p.value, scientific = T, digits = 3), round(res_bp_bc$p.value, 3), round(res_dw_bc$p.value, 3)),
  p_value_mod4n =  res_test_results_n$p_value_mod4n,
  p_value_mod4 = res_test_results$p_value
)

kable(res_test_results_bc)%>%
  kable_styling(full_width = FALSE, position = 'left')

Test	Caratteristica_testata	p_value_modbc	p_value_mod4n	p_value_mod4
Shapiro-Wilk	Normalità	9.43e-10	8.65e-12	8.11e-13
Breusch-Pagan	Omoschedasticità	0.001	0.031	5.79e-18
Durbin-Watson	Autocorrelazione	0.095	0.103	0.098

In base alle analisi fatte si decide di utilizzare il modello 4 ottenuto dopo la rimozione dell’outlier (mod4_n). Nonostante, infatti, i test statistici non passino, dai grafici diagnostici si possono osservare solo alcune deviazioni dalla normalità sulle code della distribuzione, mentre la distribuzione centrale appare conforme. Inoltre i grafici relativi ad omoschedasticità, distribuzione attorno alla media e leverage non evidenziano anomalie. Anche i valori outlier e leverage risultano avere tutti distanza di cook inferiore a 1 e pertanto non influenzano significativamente le stime della regressione.

Il modello finale è dunque il seguente

mod4_n$call

lm(formula = df_2$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2), 
    data = df_2)

Si rimuove dal dataset la variabile che era stata aggiunta con il peso trasformato

df_2$peso_trans <- NULL

3. Previsioni e risultati

Si stima, con il modello scelto, il peso di una neonata considerando una madre alla terza gravidanza (N.gravidanze =2) che partorirà alla 39esima settimana. Per le variabili non citate, ma presenti nel modello, si utilizza il valore medio per quelle quantitative (Lunghezza e Cranio) e la moda per quella qualitativa (Tipo.Parto).

mean_lunghezza <- mean(df_2$Lunghezza)
mean_cranio <- mean(df_2$Cranio)
mode_tipoparto <- names(sort(table(df_2$Tipo.parto), decreasing = TRUE))[1]


newdata <- data.frame(
  N.gravidanze = 2,
  Gestazione = 39,
  Lunghezza = mean_lunghezza,
  Cranio = mean_cranio,
  Tipo.parto = factor(mode_tipoparto, levels = levels(df_2$Tipo.parto)),
  Sesso = factor("F", levels = levels(df_2$Sesso))
)

predict(mod4_n, newdata = newdata)

       1 
3257.494

Il valore di peso predetto è 3257.494 g.

4. Visualizzazioni

Si osserva ora la dipendenza del peso dalla gestazione in base al sesso e in base al tipo di parto.

g1 <-
ggplot(data=df_2)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso),position = "jitter")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso),se=F,method = "lm")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso), col="black", se=F, method=lm)

g2 <-
ggplot(data=df_2)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Tipo.parto),position = "jitter")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Tipo.parto),se=F,method = "lm")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso), col="black", se=F, method=lm)

g1+g2

Si può osservare come in media i neonati maschi abbiano sempre peso maggiore delle neonate femmine.

Relativamente al tipo di parto si osserva invece come esso influenzi il peso del nascituro principalmente sui parti inferiori alle 37 settimane, dove il parto cesareo porta mediamente a neonati con peso maggiore. Invece per un numero di settimane di gestazione superiore a 37 le rette di regressione vanno quasi a coincidere, ad indicare che il tipo di parto perde di influenza.

Si osserva ora la dipendenza del peso dalla lunghezza e dalle dimensioni del cranio del neonato.

g1<-
ggplot(data=df_2)+
  geom_point(aes(x=Lunghezza,
                 y=Peso,
                 col=Sesso),position = "jitter")+
  geom_smooth(aes(x=Lunghezza,
                  y=Peso,
                  col=Sesso),se=F,method = "lm")+
  geom_smooth(aes(x=Lunghezza,
                  y=Peso), col="black", se=F, method=lm)

g2 <-
ggplot(data=df_2)+
  geom_point(aes(x=Cranio,
                 y=Peso,
                 col=Sesso),position = "jitter")+
  geom_smooth(aes(x=Cranio,
                  y=Peso,
                  col=Sesso),se=F,method = "lm")+
  geom_smooth(aes(x=Cranio,
                  y=Peso), col="black", se=F, method=lm)

g1+g2

Si osserva la forte dipendenza del peso dalle caratteristiche antropometriche del neonato.

Si osserva questa dipendenza anche con un grafico 3D.

scatter3d(Peso~Lunghezza + Cranio, point.col = Sesso, groups=Sesso, axis.col = c("black", "black", "black"), residuals=FALSE, surface.alpha=0.3)
rglwidget()

Si osserva anche il grafico 3D con le variabili Gestazione e N.Gravidanze.

scatter3d(Peso~Gestazione + N.gravidanze, point.col = Sesso, groups=Sesso, axis.col = c("black", "black", "black"), surface.alpha=0.3, residuals=FALSE)

rglwidget()

Si può osservare la dipendenza del peso dalla variabile Gestazione, mentre non se ne osserva una rispetto al numero di gravidanze. Tuttavia la visualizzazione grafica non tiene conto del fatto che in realtà il modello include anche tutte le altre variabili che non sono incluse nel grafico, pertanto la significatività del numero di gravidanze potrebbe non essere direttamente visualizzabile con poche variabili.

5. Conclusioni

È stato sviluppato un modello di regressione in grado di spiegare il 74% della variabilità del peso neonatale e di prevedere il peso dei neonati a partire dalle informazioni sulla madre e sulla gravidanza. In particolare, il peso risulta molto dipendente dalle caratteristiche antropometriche del neonato, dal numero di settimane di gestazione e dal numero di gravidanze pregresse della madre. Il modello sviluppato è il seguente:

summary(mod4_n)


Call:
lm(formula = df_2$Peso ~ N.gravidanze + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Sesso + I(Gestazione^2) + I(Lunghezza^2), 
    data = df_2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1168.71  -179.48   -11.03   161.56  1312.24 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -2.879e+03  9.015e+02  -3.194 0.001422 ** 
N.gravidanze     1.465e+01  4.214e+00   3.477 0.000516 ***
Gestazione       2.009e+02  6.611e+01   3.038 0.002403 ** 
Lunghezza       -1.910e+01  4.531e+00  -4.216 2.58e-05 ***
Cranio           1.006e+01  4.200e-01  23.951  < 2e-16 ***
Tipo.partoNat    2.829e+01  1.175e+01   2.409 0.016087 *  
SessoM           7.338e+01  1.091e+01   6.727 2.14e-11 ***
I(Gestazione^2) -2.143e+00  8.664e-01  -2.474 0.013437 *  
I(Lunghezza^2)   3.078e-02  4.614e-03   6.671 3.12e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 266.3 on 2490 degrees of freedom
Multiple R-squared:  0.7432,    Adjusted R-squared:  0.7424 
F-statistic: 900.7 on 8 and 2490 DF,  p-value: < 2.2e-16

Si possono fare le seguenti osservazioni:

Ogni gravidanza in più in media aumenta il peso del neonato di circa 15 g
Ogni mm in più del diametro del cranio aumenta il peso di circa 10 g
Il peso rispetto alla variabile lunghezza ha un andamento quadratico con concavità verso l’alto. Il punto di minimo corrisponde a circa 311 mm, oltre il quale all’aumentare di quest’ultima il peso cresce. Il campione a disposizione ha il minimo di lunghezza a 310 mm e media e mediana a un valore di circa 500 mm. Inoltre, dai grafici in cui è presente la variabile lunghezza si può osservare che in prossimità del valore 311 i dati sono molto pochi, infatti l’andamento che si osserva tra peso e lunghezza è di tipo crescente.
Il peso rispetto alle settimane di gestazione ha un andamento quadratico con concavità verso il basso, il cui massimo risulta a 46 settimane di gestazione. In questo caso questo valore è anche oltre il massimo dei dati a disposizione (43), per cui anche qui si osserva nei grafici e nella relazione tra peso e gestazione un andamento crescente.

Inoltre si è osservato con lo studio del dataset e del modello di regressione che i neonati maschi hanno un peso medio maggiore delle neonate femmine e che, per settimane di gestazione inferiori a circa 37, il parto cesareo in media porta a nascite con peso inferiore rispetto a quello naturale.

Modello Statistico per la Previsione del Peso Neonatale

Raffaella Todino

2025-11-17