Funzioni

gini.index = function(x) {
  ni = table(x)
  fi = ni/length(x)
  fi2 = fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
}

1. Raccolta dei Dati e Struttura del Dataset

Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

Età della madre: Misura dell’età in anni -> variabile quantitativa discreta.
Numero di gravidanze: Quante gravidanze ha avuto la madre -> variabile quantitativa discreta.
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice) -> variabile qualitativa dicotomica su scala nominale.
Durata della gravidanza: Numero di settimane di gestazione -> variabile quantitativa discreta.
Peso del neonato: Peso alla nascita in grammi -> variabile quantitativa continua.
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie -> variabili quantitative continue.
Tipo di parto: Naturale o cesareo -> variabile qualitativa su scala nominale.
Ospedale di nascita: Ospedale 1, 2 o 3 -> variabile qualitativa su scala nominale.
Sesso del neonato: Maschio (M) o femmina (F) -> variabile qualitativa dicotomica su scala nominale.

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.

Iniziamo importando il dataset e mostrando le prime 5 osservazioni:

dati <- read.csv("neonati.csv", stringsAsFactors = T)
attach(dati)
knitr::kable((head(dati,5)), "pipe")

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

2. Analisi e Modellizzazione

Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie. Un focus particolare sarà dato alla relazione tra le variabili materne (età, fumo, gravidanze precedenti) e il peso del neonato.

Analizziamo le variabili singolarmente.

VARIABILE Anni.madre

Usiamo il pacchetto psych per avere una overview più completa dei dati:

knitr::kable(psych::describe(Anni.madre))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	28.164	5.273578	28	28.102	4.4478	0	46	46	0.0427858	0.3777127	0.1054716

Da una prima analisi risulta che alcune osservazioni contengono degli errori di raccolta dati, infatti, sono presenti i valori:
* 0 -> osservazione numero 1380.
* 1 -> osservazione numero 1152.

which(Anni.madre == 0)

## [1] 1380

which(Anni.madre == 1)

## [1] 1152

Dopo esserci assicurati che questi sono gli unici valori errati presenti per questa variabile possiamo proseguire in due modi:
* eliminare le osservazioni
* sostituire Anni.madre[1380] e Anni.madre[1152] con il valore della mediana, ovvero 28 (fonte: https://labdisia.disia.unifi.it/gmm/eco/lez/sett2.pdf).

Proseguiamo con la seconda strada:

Anni.madre[1380] <- median(Anni.madre)
Anni.madre[1152] <- median(Anni.madre)
knitr::kable(psych::describe(Anni.madre))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	28.186	5.215121	28	28.1085	4.4478	13	46	33	0.1511176	-0.1055943	0.1043024

Dai dati mostrati sopra vediamo che la variabile Anni.madre ha un valore minimo di 13 e massimo di 46, con un valor medio di 28.18 e mediana pari a 28. La piccola differenza tra media e mediana è confermata dal piccolo coefficiente di skewness (0.15), che indica una leggera asimmetria positiva (valori bassi più frequenti). L’indice di curtosi è negativo (-0.10), indicando una forma platicurtica. Il coefficiente di variazione è pari al 18.48%.

boxplot(Anni.madre, col = "orange",
     main = "Boxplot Anni.madre")

plot(density(Anni.madre), col = 1,
     main = "Distribuzione Anni.madre",
     xlab = "Anni madre",
     ylab = "Densità di probabilita")
abline(v = mean(Anni.madre), col = 2)
abline(v = median(Anni.madre), col = 3)
abline(h=c(0, 0.02, 0.04, 0.06), col="grey10", lty="dotted") 
abline(v=c(10, 20, 30, 40, 50), col="grey10", lty="dotted")

VARIABILE N.gravidanze

knitr::kable(psych::describe(N.gravidanze))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	0.9812	1.280587	1	0.745	1.4826	0	12	12	2.512746	10.97822	0.0256117

Il summary mostra che il valor medio del numero di gravidanze è 1, con valori minimi di zero (quindi nessuna gravidanza precedente) e un valore massimo di 12. L’indice di skewness è pari a 2.51, indicando una simmetria positiva, il che significa che un numero basso di gravidanze è più frequente. Il coefficiente di curtosi positivo (10.97) indica un andamento leptocurtico, con valori centrali più probabili rispetto alle code. Trattandosi di una variabile discreta, possiamo vedere un grafico a barre che mostra la distribuzione del numero di gravidanze nel dataset, che mostra come il dato più frequente è zero, con 1096 donne che non hanno mai avuto una gravidanza.

my_bar <- barplot(height=table(N.gravidanze),
                  col = "orange",
                  main = "Distribuzione N.gravidanze",
                  xlab = "Numero di gravidanze precedenti",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 1200))
text(my_bar, table(N.gravidanze)+ 45, paste(table(N.gravidanze), sep="") ,cex=1)

VARIABILE Gestazione

knitr::kable(psych::describe(Gestazione))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	38.9804	1.868639	39	39.185	1.4826	25	43	18	-2.064074	8.249145	0.0373728

Il valore medio della variabile gestazione è 38.98 settimane (corrsipondente a 272 giorni, 9 mesi) e il valore mediano 39. Il coefficente di skewness è pari a -2.06, indicando un’asimmetria negativa (valori alti maggiormente frequenti); l’indice di curtosi è positivo (8.24), che indica un andamento leptocurtico.

boxplot(Gestazione, col = "orange",
     main = "Boxplot settimane di Gestazione")

my_bar <- barplot(height=table(Gestazione),
                  col = "orange",
                  main = "Distribuzione Gestazione",
                  xlab = "Settimane di gestazione",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 1200),
                  las = 2)
                  
text(my_bar, table(Gestazione)+ 45, paste(table(Gestazione), sep="") ,cex=1)

VARIABILE Peso

knitr::kable(psych::describe(Peso))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	3284.081	525.0387	3300	3302.901	459.606	830	4930	4100	-0.6466427	2.027507	10.50077

Il valore medio è 3284 grammi, il minimo 830 e il massimo 4930. La skewness è 10.50 (asimmetria positiva, con valori inferiori più probabili) e curtosi positiva (2.02, andamento leptocurtico). La deviazione standard è 525 g, che corrisponde a un coefficiente di variazione del 15%.

plot(density(Peso), col = 1,
     main = "Distribuzione Peso",
     xlab = "Peso [g]",
     ylab = "Densità di probabilita")
abline(v = mean(Peso), col = 2)
abline(v = median(Peso), col = 3)
abline(h=c(0, 4e-4, 8e-4), col="grey10", lty="dotted") 
abline(v=c(1000, 2000, 3000, 4000, 5000), col="grey10", lty="dotted")

boxplot(Peso, col = "orange",
     main = "Boxplot Peso")

VARIABILE Lunghezza

knitr::kable(psych::describe(Lunghezza))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	494.692	26.31864	500	496.45	22.239	310	565	255	-1.51379	6.479586	0.5263729

Il valore medio della lunghezza è 494 mm con una mediana di 500 mm. Il valori minimo e massimo sono rispettivamente 310 e 565 mm. La skewness è negativa (-1.5, valori maggiori più probabili) e la forma leptocurtica (indice di curtosi uguale a 6.48). la deviazione standard è 26.31 mm, che corrisponde a un coefficiente di variazione del 6%.

plot(density(Lunghezza), col = 1,
     main = "Distribuzione Lunghezza",
     xlab = "Lunghezza [mm]",
     ylab = "Densità di probabilita")
abline(v = mean(Lunghezza), col = 2)
abline(v = median(Lunghezza), col = 3)
abline(h=c(0.005, 0.01, 0.015, 0.02), col="grey10", lty="dotted") 
abline(v=c(300, 350, 400, 450, 500, 550), col="grey10", lty="dotted")

VARIABILE Cranio

knitr::kable(psych::describe(Cranio))

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
X1	1	2500	340.0292	16.42533	340	340.6845	14.826	235	390	155	-0.7845817	2.94145	0.3285066

La variabile cranio ha un valore medio uguale alla mediana, 340 mm. I valori minimo e massimo sono rispettivamente 235 e 390 mm. Ha un’asimmetria leggermente negativa (skewness pari a -0.78) e an andamento leptocurtico (indice di curtosi uguale a 2.94). La deviazione standard è 16.42, che corrisponde a un coefficiente di variazione del 4.7 %.

plot(density(Cranio), col = 1,
     main = "Distribuzione Diametro cranio",
     xlab = "Diametro Cranio [mm]",
     ylab = "Densità di probabilita")
abline(v = mean(Cranio), col = 2)
abline(v = median(Cranio), col = 3)
abline(h=c(0, 0.01, 0.02), col="grey10", lty="dotted") 
abline(v=c(250, 300, 350, 400), col="grey10", lty="dotted")

VARIABILE Fumatrici

Fumatrici_ni = table(Fumatrici)
Fumatrici_fi = Fumatrici_ni/length(Fumatrici)

Fumatrici_distr = cbind(Fumatrici_ni, Fumatrici_fi)
knitr::kable(Fumatrici_distr, "pipe")

	Fumatrici_ni	Fumatrici_fi
0	2396	0.9584
1	104	0.0416

print(gini.index(Fumatrici))

## [1] 0.1594778

my_bar <- barplot(height=table(Fumatrici),
                  col = "orange",
                  main = "Frequenze assolute Fumatrici",
                  xlab = "",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 3000))
                  
text(my_bar, table(Fumatrici)+ 100, paste(table(Fumatrici), sep=" ") ,cex=1)

Il numero di non fumatrici è 23 volte più alto rispetto alla fumatrici, che rappresentano solo il 4.1 % del campione. L’indice di Gini è pari a 0.16 (scarsa eterogeneità).

VARIABILE Tipo.parto

Tipo.parto_ni = table(Tipo.parto)
Tipo.parto_fi = Tipo.parto_ni/length(Tipo.parto)

Tipo.parto_distr = cbind(Tipo.parto_ni, Tipo.parto_fi)
knitr::kable(Tipo.parto_distr, "pipe")

	Tipo.parto_ni	Tipo.parto_fi
Ces	728	0.2912
Nat	1772	0.7088

print(gini.index(Tipo.parto))

## [1] 0.8256102

my_bar <- barplot(height=table(Tipo.parto),
                  col = "orange",
                  main = "Frequenze assolute Tipo.parto",
                  xlab = "Tipo parto (cesareo/naturale)",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 2500))
                  
text(my_bar, table(Tipo.parto)+ 100, paste(table(Tipo.parto), sep=" ") ,cex=1)

Il numero di parti naturali è 2.4 volte il numero ci parti cesarei, che rappresentano il 29% del campione. L’indice di Gini è 0.82.

VARIABILE Ospedale

Ospedale_ni = table(Ospedale)
Ospedale_fi = Ospedale_ni/length(Ospedale)

Ospedale_distr = cbind(Ospedale_ni, Ospedale_fi)
knitr::kable(Ospedale_distr, "pipe")

	Ospedale_ni	Ospedale_fi
osp1	816	0.3264
osp2	849	0.3396
osp3	835	0.3340

print(gini.index(Ospedale))

## [1] 0.9998683

my_bar <- barplot(height=table(Ospedale),
                  col = "orange",
                  main = "Frequenze assolute per ogni ospedale",
                  xlab = "Ospedale",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 1000))
                  
text(my_bar, table(Ospedale)+ 100, paste(table(Ospedale), sep=" ") ,cex=1)

Le osservazioni sono distribuite in maniera molto equa tra le 3 modalità, con un indice di Gini di 0.9998, indicando massima eterogeneità.

VARIABILE Sesso

Sesso_ni = table(Sesso)
Sesso_fi = Sesso_ni/length(Sesso)

Sesso_distr = cbind(Sesso_ni, Sesso_fi)
knitr::kable(Sesso_distr, "pipe")

	Sesso_ni	Sesso_fi
F	1256	0.5024
M	1244	0.4976

print(gini.index(Sesso))

## [1] 0.999977

my_bar <- barplot(height=table(Sesso),
                  col = "orange",
                  main = "Frequenze assolute Sesso",
                  xlab = "Sesso (M/F)",
                  ylab = "Frequenza assoluta",
                  ylim=c(0, 1500))
                  
text(my_bar, table(Sesso)+ 100, paste(table(Sesso), sep=" ") ,cex=1)

La variabile sesso è distribuita equamente tra M e F, con un indice di Gini di 1 (massima eterogeità).

Variabile Peso in relazione alle altre variabili materne (età, fumo, gravidanze precedenti).

Facciamo un’analisi multidimensionale della variabile Peso in relazione alle altre variabili materne (età, fumo, gravidanze precedenti). Mostriamo i boxplot, gli scatterplot e calcoliamo dei summary con il pacchetto dplyr.

Peso~Anni.madre

dati$Anni.madre_CL <- cut(Anni.madre, breaks = c (12, 24, 35, 46))
knitr::kable(table(dati$Anni.madre_CL), "pipe")

Var1	Freq
(12,24]	595
(24,35]	1687
(35,46]	218

boxplot(Peso~dati$Anni.madre_CL, ,
                  col = "orange",
                  main = "Boxplot Peso~Anni.madre",
                  xlab = "Classi di età Anni.madre")

library(dplyr)

## 
## Attaching package: 'dplyr'

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

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

dati %>%
  group_by(Anni.madre_CL) %>%
  summarize(media = mean(Peso), devstd = sd(Peso))

## # A tibble: 3 × 3
##   Anni.madre_CL media devstd
##   <fct>         <dbl>  <dbl>
## 1 (12,24]       3277.   469.
## 2 (24,35]       3294.   526.
## 3 (35,46]       3223.   645.

Dal boxplot e dal summary si vede che il valor medio della variabile peso per la fascia di età (35,46] subisce un leggero calo (di 54 g e 71 g rispetto alle classi (12,24] e (24,35]); l’aspetto di maggior rilievo è tuttavia l’aumento della deviazione standard; passiamo infatti da un coefficiente di variazione del 14% per la classe (12,24], al 15.9% della classe di età (24,35], per finire al 19.98% per ultima classe di età.

plot(Peso~jitter(Anni.madre,3), pch = 1, cex = 0.5, col = "purple")

cor(Peso, Anni.madre)

## [1] -0.02377346

Dallo scatterplot si vede un andamento orizzontale, confermato dal basso coefficiente di correlazione -0.023.

Peso~Fumatrici

boxplot(Peso~Fumatrici,
        col = "orange",
        main = "Boxplot Peso~Fumatrici",
        xlab = "Fumatrici (SI/NO)")

dati %>%
  group_by(Fumatrici) %>%
  summarize(media = mean(Peso), devstd = sd(Peso))

## # A tibble: 2 × 3
##   Fumatrici media devstd
##       <int> <dbl>  <dbl>
## 1         0 3286.   527.
## 2         1 3236.   479.

Nella classificazione della variabile peso in base alla variabile dicotomica Fumatrici, possiamo notare una diminuzione del peso medio da 3286 g a 3236, corrispondente a una diminuzione percentuale dell’ 1.5 % passando da “non-fumatrici” a “fumatrici” e una diminuzione del coefficiente di variazione che passa dal 16% al 14.7%.

Peso~N.gravidanze

boxplot(Peso~dati$N.gravidanze, ,
                  col = "orange",
                  main = "Boxplot Peso~N.gravidanze",
                  xlab = "Numero di gravidanze precedenti")

Peso_N.gravidanze <- dati %>%
  group_by(N.gravidanze) %>%
  summarize(media = mean(Peso), devstd = sd(Peso))
Peso_N.gravidanze

## # A tibble: 13 × 3
##    N.gravidanze media devstd
##           <int> <dbl>  <dbl>
##  1            0 3279.   519.
##  2            1 3275.   510.
##  3            2 3319.   504.
##  4            3 3289.   595.
##  5            4 3398.   685.
##  6            5 3316.   446.
##  7            6 3063.   829.
##  8            7 2220     NA 
##  9            8 3076.   944.
## 10            9 3455    431.
## 11           10 2973.   107.
## 12           11 3400     NA 
## 13           12 3350     NA

y1 <- as.data.frame(t(Peso_N.gravidanze[2]))
x1 <- seq(0,12,1)
plot(x1, y1, ylim=c(1000, 3500))
lines(x1, y1, col = "red")

Considerato che i dati del campione per i quali il numero di gravidanze è superiore a 3 è molto ridotto, non è possibile individuare in questo momento un trend chiaro che leghi la variabile peso al numero di gravidanze precedenti.

plot(Peso~jitter(N.gravidanze,3), pch = 1, cex = 0.5, col = "purple")

cor(Peso, N.gravidanze)

## [1] 0.0024073

La nuvola di punti è molto addensata nella regione 0-1, quella in cui ci sono la maggior parte delle nascite; tuttavia graficamente non è possibile individuare un trend chiaro, confermato dal coefficiente di correlazione (0.002).

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. Ad esempio, ci aspettiamo che il fumo materno e le gravidanze multiple abbiano effetti negativi sul peso alla nascita, mentre una maggiore durata della gestazione potrebbe aumentare il peso del neonato.

Matrice delle correlazioni

Con la funzione pairs, mostriamo la matrice delle correlazioni. Creiamo un nuovo dataframe per vedere la matrice delle correlazioni tra la variabile Peso e le altre variabili quantitative (Anni.madre, N,gravidanze, Gestazione, Lunghezza, Cranio).

?pairs

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

library(dplyr)

dati_quant <- data.frame(Peso, Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio)
pairs(dati_quant,lower.panel=panel.cor, upper.panel=panel.smooth)

Dalla matrice delle correlazioni possiamo fare queste considerazioni:
- le variabili che mostrano il coefficiente di correlazione maggiore con la variabile Peso sono:
1) Gestazione con un coefficente di 0.59.
2) Lunghezza con correlazione positiva e coefficiente 0.80.
3) Diametro del cranio con coefficiente di correlazione 0.70.
- bisognerà verificare problemi di multicollinearità calcolando i VIF, considerati i coefficenti di correlzione tra:
1) Anni.Madre e N.Gravidanze.
2) Lunghezza e Gestazione.
3) Lunghezza e Cranio.
4) Gestazione e Cranio.

Variabili qualitative

Per quanto riguarda le variabili qualitative (Fumatrici, Tipo.parto, Ospedale e Sesso) vediamo i boxplot condizionati di Peso e della variabile qualitativa.

Peso~Fumatrici

Partendo dalla variabile Fumatrici:

par(mfrow = c(1,2))
boxplot(Peso, main = "Peso", col = "orange")
boxplot(Peso~Fumatrici, main = "Peso~Fumatrici", col = "orange")

Facciamo un t.test per saggiare l’ipotesi di uguaglianza delle medie tra gruppi indipendenti:

t.test(Peso~Fumatrici)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

Il p_value è 0.30, quindi siamo ampiamente nella zona di non rifiuto; pertanto non possiamo rifiutare l’ipotesi nulla H0 secondo la quale la differenza delle medie tra peso dei neonati da madri fumatrici e non fumatrici non è significativamente diversa da zero.

plot(density(rt(100000, 114.1)), main = "t test di differenza tra medie di Peso~Fumatrici")
abline(v = qt(0.025, 114), col = "red")
abline(v = qt(0.975, 114), col = "red")
points(x = 1.034, y = 0, col = "blue", pch = 20, cex = 3)

Inspettatamente la variabile Fumatrici non ha un’influenza significativa sul peso del neonato; tuttavia questo risultato potrebbe essere anche dovuto al fatto che il numero di madri fumatrici è pari a 104, ovvero al 4.16% del campione. In letteratura ci sono studi che smentiscono il t.test effettuato sui dati del campione (https://bmcpregnancychildbirth.biomedcentral.com/articles/10.1186/s12884-018-1694-4#:~:text=Compared%20with%20infants%20born%20to%20nonsmoking%20mothers%2C%20mean%20birth%20weight,cigarettes%20per%20day%20during%20pregnancy.).

count(dati, Fumatrici)

##   Fumatrici    n
## 1         0 2396
## 2         1  104

Peso~Tipo.parto

par(mfrow = c(1,2))
boxplot(Peso, main = "Peso", col = "orange")
boxplot(Peso~Tipo.parto, main = "Tipo.parto", col = "orange")

Ripetiamo anche in questo caso il t.test:

t.test(Peso~Tipo.parto)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -46.27992  40.54037
## sample estimates:
## mean in group Ces mean in group Nat 
##          3282.047          3284.916

Anche in questo caso il p_value è molto alto e non possiamo rifiutare l’ipotesi nulla.

plot(density(rt(100000, 1493)), main = "t test di differenza tra medie di Peso~Tipo.parto")
abline(v = qt(0.025, 1493), col = "red")
abline(v = qt(0.975, 1493), col = "red")
points(x = -0.12968, y = 0, col = "blue", pch = 20, cex = 3)

Peso~Ospedale

par(mfrow = c(1,2))
boxplot(Peso, main = "Peso", col = "orange")
boxplot(Peso~Ospedale, main = "Ospedale", col = "orange")

In questi caso le modalità della variabile qualitative ospedale sono 3, quindi utilizzeremo un t.test per saggiare l’ipotesi di differenze tra medie di gruppi indipendenti della variabile peso utilizzando la correzione di Bonferroni:

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

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Peso and Ospedale 
## 
##      osp1 osp2
## osp2 0.99 -   
## osp3 0.33 0.33
## 
## P value adjustment method: holm

In tutti e 3 i confronti incrociati, i p_value hanno valori molto alti e pertando non possiamo rifiutare l’ipotesi nulla H0: le differenze tra le medie della variabile peso non sono significativamente diverse da zero al variare dell’ospedale.

Peso~Sesso

par(mfrow = c(1,2))
boxplot(Peso, main = "Peso", col = "orange")
boxplot(Peso~Sesso, main = "Sesso", col = "orange")

Facciamo anche in questo caso il t.test:

t.test(Peso~Sesso)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

Il p_value è minimo (<2.2e-16), pertanto possiamo rifiutare l’ipotesi nulla H0 e dire che con la media del Peso del neonato è significativamente diversa a seconda del sesso.

plot(density(rt(100000, 2490.7)), main = "t test di differenza tra medie di Peso~Sesso", xlim = c(-20,4))
abline(v = qt(0.025, 2490.7), col = "red")
abline(v = qt(0.975, 2490.7), col = "red")
points(x = -12.106, y = 0, col = "blue", pch = 20, cex = 3)

Modello di regressione

Per la creazione del modello inseriamo inizialmente tutte le variabili, poi faremo un’analisi per semplificarlo:

mod1 <- lm(Peso~Anni.madre+N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Tipo.parto+Ospedale+Sesso)
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1123.3  -181.2   -14.6   160.7  2612.6 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.1677   141.3977 -47.633  < 2e-16 ***
## Anni.madre        0.7983     1.1463   0.696   0.4862    
## N.gravidanze     11.4118     4.6665   2.445   0.0145 *  
## Fumatrici       -30.1567    27.5396  -1.095   0.2736    
## Gestazione       32.5265     3.8179   8.520  < 2e-16 ***
## Lunghezza        10.2951     0.3007  34.237  < 2e-16 ***
## Cranio           10.4725     0.4261  24.580  < 2e-16 ***
## Tipo.partoNat    29.5027    12.0848   2.441   0.0147 *  
## Ospedaleosp2    -11.2216    13.4388  -0.835   0.4038    
## Ospedaleosp3     28.0984    13.4972   2.082   0.0375 *  
## SessoM           77.5473    11.1779   6.938 5.07e-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.1 on 10 and 2489 DF,  p-value: < 2.2e-16

Gestazione, Lunghezza e Cranio sono le variabili quantitative con i coefficienti di regressione che hanno un peso maggiore all’interno del modello e un p_value molto piccolo. Per quanto riguarda il tipo di parto, si evidenzia che il parto naturale ha un coefficente di regressione beta positivo (29.50) e un beta molto significativo (0.0147). La variabile Sesso ha un p_value molto piccolo e assegna un coefficiente beta di 77.55 nel caso il neonato sia di sesso maschile.
La variabile Anni.madre ha un p_value molto alto che suggerisce una scarsa influenza sul modello, così come le variabili Fumatrici e Ospedale.

Il valore di R^2 adjusted è pari al 72%, ovvero un modello che spiega il 72% della varianza del campione.

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.

In questa fase proveremo a ottimizzare il modello, provando a semplificarlo, riducendo il numero di variabili e provando a massimizzare il coefficiente R^2.

Come primo step proviamo ad eliminare dal modello le variabili che riteniamo non significative.

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

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Sesso)
## 
## 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

Al modello 2 sono state sottratte 3 variabili (Anni.madre, Fumatrici e Ospedale); il nuovo modello ha un coefficiente R^2 simile al modello1 (72%) ma con 3 variabili in meno; per il criterio del rasoio di Occam questo modello è da considerarsi migliore rispetto al modello 1. I p_value sono cambiati significativamente (la variazione maggiore si osserva su N.gravidanze il cui p_value è passato da 0.014 del modello 1 a 0.0033 del modello 2).

Proviamo a semplificare ulteriormente eliminado la variabile N.gravidanze.

mod3 <- update(mod2,~.-N.gravidanze)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1118.43  -185.56   -16.07   161.53  2626.18 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6675.8084   135.7769 -49.167  < 2e-16 ***
## Gestazione       31.1917     3.7821   8.247 2.59e-16 ***
## Lunghezza        10.2412     0.3008  34.049  < 2e-16 ***
## Cranio           10.6397     0.4242  25.080  < 2e-16 ***
## Tipo.partoNat    29.1062    12.1113   2.403   0.0163 *  
## SessoM           79.0670    11.2010   7.059 2.17e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2494 degrees of freedom
## Multiple R-squared:  0.7267, Adjusted R-squared:  0.7262 
## F-statistic:  1326 on 5 and 2494 DF,  p-value: < 2.2e-16

Il modello 3 ha un coefficiente R^2 praticamente uguale al modello 2 ma con una variabile in meno. I p_value non sono cambiati in modo significativo rispetto al caso precedente; pertanto a parità di prestazioni, il modello 3 è attualmente da preferire.

Osservando in dettaglio gli scatterplot, si osserva una flessione della nuvola di punti nella curva Peso~Gestazione, che potrebbe essere sinonimo di una non-linearità. Proviamo ad inserire nel modello una dipendenza di tipo lineare e quadratica di Peso da Gestazione:

mod4 <- update(mod3, ~.+I(Gestazione^2))
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso + I(Gestazione^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1120.27  -184.61   -15.21   164.47  2648.27 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4704.6435   899.0411  -5.233 1.81e-07 ***
## Gestazione        -78.8252    49.7474  -1.585   0.1132    
## Lunghezza          10.3420     0.3040  34.024  < 2e-16 ***
## Cranio             10.7344     0.4260  25.195  < 2e-16 ***
## Tipo.partoNat      28.6370    12.1037   2.366   0.0181 *  
## SessoM             76.9383    11.2333   6.849 9.32e-12 ***
## I(Gestazione^2)     1.4686     0.6622   2.218   0.0267 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2493 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7266 
## F-statistic:  1108 on 6 and 2493 DF,  p-value: < 2.2e-16

Il modello 4 ha lo stesso R^2 del modello 3; sembra effettivamente esserci una non linearità (p_value = 0.0267), ma l’introduzione della nuova variabile ha aumentato il p_value di gestazione da 2.59e-16 a 0.1132; quindi il modello 3 è da preferire.

Per quanto riguarda le ipotesi di modello congiunto, possiamo provare a vedere le relazione congiunta tra Gestazione/Cranio e Gestazione/Lunghezza.

mod5 <- lm(Peso~Tipo.parto+Sesso+Gestazione*Lunghezza+Gestazione*Cranio)
summary(mod5)

## 
## Call:
## lm(formula = Peso ~ Tipo.parto + Sesso + Gestazione * Lunghezza + 
##     Gestazione * Cranio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1105.8  -183.6   -14.7   162.4  2677.3 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -328.28121 1108.04975  -0.296   0.7670    
## Tipo.partoNat          27.80805   12.03760   2.310   0.0210 *  
## SessoM                 73.14955   11.18564   6.540 7.46e-11 ***
## Gestazione           -138.20109   29.55732  -4.676 3.09e-06 ***
## Lunghezza               9.17478    3.76073   2.440   0.0148 *  
## Cranio                 -7.45568    6.43949  -1.158   0.2471    
## Gestazione:Lunghezza    0.03328    0.09749   0.341   0.7328    
## Gestazione:Cranio       0.47458    0.16648   2.851   0.0044 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273 on 2492 degrees of freedom
## Multiple R-squared:  0.7304, Adjusted R-squared:  0.7296 
## F-statistic: 964.3 on 7 and 2492 DF,  p-value: < 2.2e-16

Provando ad eliminare la relazione congiunta Gestazione/Lunghezza:

mod6 <- lm(Peso~Tipo.parto+Sesso+Lunghezza+Gestazione*Cranio)
summary(mod6)

## 
## Call:
## lm(formula = Peso ~ Tipo.parto + Sesso + Lunghezza + Gestazione * 
##     Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1106.68  -183.65   -14.15   162.89  2680.73 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -319.73011 1107.57033  -0.289  0.77285    
## Tipo.partoNat       27.81860   12.03542   2.311  0.02089 *  
## SessoM              73.31415   11.17327   6.562 6.45e-11 ***
## Lunghezza           10.45458    0.30111  34.720  < 2e-16 ***
## Gestazione        -138.28043   29.55117  -4.679 3.03e-06 ***
## Cranio              -9.30632    3.47545  -2.678  0.00746 ** 
## Gestazione:Cranio    0.52232    0.09034   5.782 8.31e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273 on 2493 degrees of freedom
## Multiple R-squared:  0.7303, Adjusted R-squared:  0.7297 
## F-statistic:  1125 on 6 and 2493 DF,  p-value: < 2.2e-16

Il modello 6 ha un R^2 simile agli altri modelli, ma comunque più alto, pertanto al momento da preferire.
Se proviamo a rimuovere Tipo.parto:

mod7 <- lm(Peso~Sesso+Lunghezza+Gestazione*Cranio)
summary(mod7)

## 
## Call:
## lm(formula = Peso ~ Sesso + Lunghezza + Gestazione * Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1125.48  -181.04   -13.99   163.92  2682.19 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -249.1238  1108.1125  -0.225  0.82214    
## SessoM              73.3078    11.1830   6.555 6.72e-11 ***
## Lunghezza           10.4220     0.3010  34.620  < 2e-16 ***
## Gestazione        -139.4557    29.5725  -4.716 2.54e-06 ***
## Cranio              -9.4246     3.4781  -2.710  0.00678 ** 
## Gestazione:Cranio    0.5262     0.0904   5.821 6.62e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7292 
## F-statistic:  1347 on 5 and 2494 DF,  p-value: < 2.2e-16

Il modello 7 ha un R^2 adjusted di 0.7292, mentre il modello 6 di 0.7297, ma il modello 7 ha una variabile in meno, pertanto preferibile. Usiamo il criterio di valutazione AIC e BIC.

AIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7)

##      df      AIC
## mod1 12 35172.09
## mod2  8 35175.16
## mod3  7 35181.82
## mod4  8 35178.89
## mod5  9 35152.40
## mod6  8 35150.52
## mod7  7 35153.87

BIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7)

##      df      BIC
## mod1 12 35241.97
## mod2  8 35221.75
## mod3  7 35222.59
## mod4  8 35225.48
## mod5  9 35204.82
## mod6  8 35197.11
## mod7  7 35194.64

Il criterio AIC predilige il modello 6 più complesso, mentre il criterio BIC il modello 7 più semplice. Diamo priorità al criterio BIC e scegliamo il modello 7.

Possiamo provare anche ad utilizzare il metodo stepwise per vedere il modello migliore offerto dall’algoritmo:

stepwise.model <- step(mod1, direction="both",k=2)

## Start:  AIC=28075.39
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36392 186809099 28074
## - Fumatrici     1     89979 186862686 28075
## <none>                      186772707 28075
## - Tipo.parto    1    447233 187219939 28079
## - N.gravidanze  1    448762 187221469 28079
## - Ospedale      2    686397 187459103 28081
## - Sesso         1   3611588 190384294 28121
## - Gestazione    1   5446558 192219264 28145
## - Cranio        1  45338051 232110758 28617
## - Lunghezza     1  87959790 274732497 29038
## 
## Step:  AIC=28073.88
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28073
## <none>                      186809099 28074
## + Anni.madre    1     36392 186772707 28075
## - Tipo.parto    1    448222 187257321 28078
## - Ospedale      2    692738 187501837 28079
## - N.gravidanze  1    633756 187442855 28080
## - Sesso         1   3618736 190427835 28120
## - Gestazione    1   5412879 192221978 28143
## - Cranio        1  45588236 232397335 28618
## - Lunghezza     1  87950050 274759149 29036
## 
## Step:  AIC=28073.1
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      186899996 28073
## + Fumatrici     1     90897 186809099 28074
## + Anni.madre    1     37311 186862686 28075
## - Tipo.parto    1    440684 187340680 28077
## - Ospedale      2    701680 187601677 28078
## - N.gravidanze  1    610840 187510837 28079
## - Sesso         1   3602797 190502794 28119
## - Gestazione    1   5346781 192246777 28142
## - Cranio        1  45632149 232532146 28617
## - Lunghezza     1  88355030 275255027 29039

summary(stepwise.model)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## N.gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16

AIC(stepwise.model, mod7)

##                df      AIC
## stepwise.model 10 35169.79
## mod7            7 35153.87

BIC(stepwise.model, mod7)

##                df      BIC
## stepwise.model 10 35228.03
## mod7            7 35194.64

Preferiamo il modello 7 rispetto al modello creato con il metodo stepwise, confermato dai criteri AIC e BIC.
Vediamo se ci sono problemi di multicollinearità:

car::vif(mod7, "predictor")

## GVIFs computed for predictors

##                GVIF Df GVIF^(1/(2*Df)) Interacts With
## Sesso      1.047119  1        1.023288           --  
## Lunghezza  2.101626  1        1.449698           --  
## Gestazione 2.095429  3        1.131216         Cranio
## Cranio     2.095429  3        1.131216     Gestazione
##                         Other Predictors
## Sesso      Lunghezza, Gestazione, Cranio
## Lunghezza      Sesso, Gestazione, Cranio
## Gestazione              Sesso, Lunghezza
## Cranio                  Sesso, Lunghezza

Tutti i VIF sono inferiori a 5 e pertanto possiamo escludere problemi di multicollinearità.

Analisi della Qualità del Modello

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

Procediamo facendo l’analisi dei residui del modello, che deve rispettare i criteri di normalità, omoschedasticità e indipendenza tra i residui.

par(mfrow = c (2,2))
plot(mod7, pch = 1, cex = 0.5, col = "orange")

Dal primo grafico si vede che i residui sono distribuiti attorno allo zero, suggerendo che la maggior parte delle osservazioni è stata “captata” dal modello.
Per quanto riguarda il secondo grafico, idealmente i punti dovrebbero stare tutti sulla bisettrice; questo è valido per tutti i valori centrali ma si vede una distorsione all’inizio e alla fine della coda. Lo Shapiro test potrà confermare/smentire l’ipotesi di normalità dei residui.
Nel terzo grafico i punti sembrano assumere una forma che suggerisce omoschedasticità.
Il quarto grafico mostra i valori di leva; l’osservazione 1551 supera il valore di soglia 0.5 ma non il valore 1.

Procediamo con il test di normalità:

shapiro.test(residuals(mod7))

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

Il p_value è molto basso e pertando dobbiamo rifiutare l’ipotesi di normalità:

plot(density(residuals(mod7)))

Vediamo come vanno gli altri test:

lmtest::bptest(mod7)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod7
## BP = 84.129, df = 5, p-value < 2.2e-16

Il test di Breusch-Pagan restituisce un valore molto basso e quindi dobbiamo rifiutare l’ipotesi di omoschedasticità.

lmtest::dwtest(mod7)

## 
##  Durbin-Watson test
## 
## data:  mod7
## DW = 1.9598, p-value = 0.1574
## alternative hypothesis: true autocorrelation is greater than 0

Il Durbin-Watson test resitutisce un p_value di 0.1574 e quindi non possiamo rifiutare H0, pertanto i residui tra loro risultano indipendenti.

Facciamo l’analisi dei valori di leva:

lev <- hatvalues(mod7)
plot(lev, pch = 1, cex = 0.5, col = "orange")
p = sum(lev)
n = length(lev)
soglia = 2*p/n
abline(h = soglia, col = 2)

Si mostrano molti valori di leva nello spazio di regressione.

lev[lev>soglia]

##          15          34          36          67          96         101 
## 0.006662141 0.006279801 0.007212227 0.005426227 0.005630577 0.008131732 
##         106         117         131         151         155         190 
## 0.026268199 0.004869023 0.008412940 0.014212651 0.007352563 0.005411964 
##         205         206         220         249         305         310 
## 0.005463922 0.011641172 0.007464761 0.005754283 0.005582291 0.069150154 
##         312         315         340         378         445         486 
## 0.018074792 0.007390396 0.004924712 0.037966967 0.010029108 0.004864763 
##         492         532         556         565         587         592 
## 0.008962035 0.005163885 0.004906628 0.005749134 0.010296057 0.006400975 
##         615         638         684         697         726         748 
## 0.005875917 0.006669026 0.008813382 0.005940044 0.005207166 0.011875031 
##         750         765         805         895         928         947 
## 0.007484866 0.006658296 0.032331609 0.007245454 0.064330044 0.009033009 
##         956         991        1014        1067        1091        1130 
## 0.008522457 0.005361778 0.008313800 0.009558420 0.012876407 0.008435566 
##        1134        1157        1181        1188        1200        1238 
## 0.006131577 0.004849573 0.007613774 0.006437339 0.005577303 0.005377858 
##        1248        1273        1283        1293        1294        1356 
## 0.031703595 0.007523786 0.005045700 0.005792517 0.004828409 0.006613590 
##        1357        1385        1395        1400        1420        1428 
## 0.007322943 0.018458854 0.004860999 0.005565973 0.005162430 0.008194394 
##        1429        1551        1553        1556        1593        1606 
## 0.034041960 0.049713601 0.008491138 0.007022197 0.004856215 0.004941190 
##        1610        1619        1634        1686        1693        1701 
## 0.009675088 0.021318459 0.004912072 0.009107394 0.005632537 0.010675601 
##        1712        1735        1780        1806        1809        1827 
## 0.007046813 0.005052455 0.104978777 0.004875435 0.014807386 0.006025296 
##        1858        1868        1959        1977        2016        2040 
## 0.004976747 0.005821904 0.004853281 0.007349682 0.008643615 0.011855707 
##        2062        2087        2089        2114        2115        2120 
## 0.004993987 0.005004512 0.006437323 0.023588493 0.012412038 0.053672475 
##        2140        2149        2175        2180        2200        2215 
## 0.007748714 0.024089390 0.097179618 0.005077950 0.014684667 0.005553188 
##        2216        2224        2225        2257        2307        2318 
## 0.008016388 0.005781860 0.006053966 0.005966523 0.026422811 0.005451828 
##        2337        2382        2391        2408        2433        2437 
## 0.005605511 0.004882716 0.006088849 0.013137382 0.005000669 0.058284181 
##        2452        2458        2478 
## 0.109201379 0.009805826 0.005823114

Per quanto riguarda gli outliers:

plot(rstudent(mod7), pch = 1, cex = 0.5, col = "orange")
abline(h = c(-2,2), col = 2)

car::outlierTest(mod4)

##       rstudent unadjusted p-value Bonferroni p
## 1551 10.094831         1.6380e-23   4.0949e-20
## 155   4.999623         6.1434e-07   1.5359e-03
## 1306  4.709843         2.6151e-06   6.5378e-03

Ci sono 3 osservazioni outliers (1551, 155 e 1306). Vediamo se superano la distanza di Cook:

cooks <- cooks.distance(mod7)
plot(cooks, pch = 1, cex = 0.75, col = "orange")

Come già visto precedentemente, l’osservazione 1551 supera il valore di soglia 0.5, ma non quello di allarme 1, pertanto lo possiamo mantenere nel modello.

Calcolo RMSE:

sqrt(sum(mean((Peso-predict(mod7))^2)))

## [1] 272.8809

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, non fumatrice, che partorirà alla 39esima settimana.

predict(mod7, newdata=data.frame(Sesso = 'F', Lunghezza = mean(Lunghezza), Gestazione = 39, Cranio = mean(Cranio)))

##        1 
## 3240.937

mean(Peso)

## [1] 3284.081

La funzione predict restituisce il valore 3240.081 grammi per il peso della neonata.

4. Visualizzazioni

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello.

Mostriamo l’impatto del numero di settimane di gestazione sul peso del neonato. Creiamo un grafico che mostra il peso del neonato in funzione del numero di settimane di gestazione a seconda che sia maschio o femmina.

library(ggplot2)
ggplot(data=dati)+
  geom_point(aes(x=Gestazione, y=Peso, col=Sesso),position='jitter')+
  geom_smooth(aes(x=Gestazione, y=Peso, col=Sesso), se=F, method = "lm")+
    labs(title='Peso in funzione del numero di settimane di gestazione e del sesso')

## `geom_smooth()` using formula = 'y ~ x'

Allo stesso modo possiamo mostrare il peso per le variabili lunghezza e cranio dipendentemente dal sesso:

ggplot(data=dati)+
  geom_point(aes(x=Lunghezza, y=Peso, col=Sesso),position='jitter')+
  geom_smooth(aes(x=Lunghezza, y=Peso, col=Sesso), se=F, method = "lm")+
    labs(title='Peso in funzione della lunghezza e del sesso')

## `geom_smooth()` using formula = 'y ~ x'

ggplot(data=dati)+
  geom_point(aes(x=Cranio, y=Peso, col=Sesso),position='jitter')+
  geom_smooth(aes(x=Cranio, y=Peso, col=Sesso), se=F, method = "lm")+
    labs(title='Peso in funzione del diametro del cranio e del sesso')

## `geom_smooth()` using formula = 'y ~ x'

Conclusioni

Il modello di regressione lineare multivariabile scelto utilizza le variabili sesso, lunghezza, diametro del cranio e settimane di gestazione con un valore R^2 adjuster pari a 0.7292.
Le variabili età della madre e tabagismo non sono state incluse nel modello perchè non apportavano un significativo cambiamento (seppure le informazioni provenienti dalla letteratura sono differenti).

Progetto_03_StatisticaInferenziale

Giuseppe Frascolla

2024-11-25

Funzioni

1. Raccolta dei Dati e Struttura del Dataset

2. Analisi e Modellizzazione

Analisi Preliminare

VARIABILE Anni.madre

VARIABILE N.gravidanze

VARIABILE Gestazione

VARIABILE Peso

VARIABILE Lunghezza

VARIABILE Cranio

VARIABILE Fumatrici

VARIABILE Tipo.parto

VARIABILE Ospedale

VARIABILE Sesso

Variabile Peso in relazione alle altre variabili materne (età, fumo, gravidanze precedenti).

Peso~Anni.madre

Peso~Fumatrici

Peso~N.gravidanze

Creazione del Modello di Regressione

Matrice delle correlazioni

Variabili qualitative

Peso~Fumatrici

Peso~Tipo.parto

Peso~Ospedale

Peso~Sesso

Modello di regressione

Selezione del Modello Ottimale

Analisi della Qualità del Modello

3. Previsioni e Risultati

4. Visualizzazioni

Conclusioni