Modello di previsione peso neonati
Riccardo Siega
2024-10-21
1. Raccolta dati e struttura dataset
In questo documento verrà creato un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali.
Il dataframe ‘neonati’ contiene informazioni su 2,500 neonati provenienti da tre ospedali. Le variabili raccolte includono:
Età della madre (anni).
Numero di gravidanze (antecedenti l’attuale).
Fumo materno (0=non fumatrice, 1=fumatrice).
Durata della gravidanza (settimane).
Peso neonato alla nascita (grammi).
Lunghezza e diametro del cranio (cm), misurabili anche durante la gravidanza tramite ecografie.
Tipo di parto (Nat, Ces).
Ospedale di nascita (osp1, osp2, osp3).
Sesso del neonato (M, F)
Per prima cosa, dobbiamo importare il dataframe in R. Il working path dovrà essere modificato dall’utente finale per personalizzarlo in base alle proprie esigenze.
path = "/Users/Riccardo/My Drive/Master Data Science/Mod3_Statistica_inferenziale/Cap9_Project"
setwd(path)
neonati <- read.csv("https://drive.google.com/uc?export=download&id=1ChfwftuOSH-WLIto_1AvV-_sQIksGeTq", sep = ",")
attach(neonati)
head(neonati) # visualizza le prime righe## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
## 1 26 0 0 42 3380 490 325 Nat osp3 M
## 2 21 2 0 39 3150 490 345 Nat osp1 F
## 3 34 3 0 38 3640 500 375 Nat osp2 M
## 4 28 1 0 41 3690 515 365 Nat osp2 M
## 5 20 0 0 38 3700 480 335 Nat osp3 F
## 6 32 0 0 40 3200 495 340 Nat osp2 Fstr(neonati) # visualizza la struttura del dataframe## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : chr "Nat" "Nat" "Nat" "Nat" ...
## $ Ospedale : chr "osp3" "osp1" "osp2" "osp2" ...
## $ Sesso : chr "M" "F" "M" "M" ...sum(is.na(neonati)) # verifica la presenza di valori mancanti## [1] 0Importiamo quindi librerie e definiamo le funzioni che ci serviranno.
# Installazione librerie
install.packages("ggplot2")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("moments")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("dplyr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("tidyr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("clipr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("gghalves")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("lubridate")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packagesinstall.packages("RColorBrewer")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpB5m7k8/downloaded_packageslibrary(ggplot2)
library(moments)
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, unionlibrary(tidyr)
library(clipr)## Welcome to clipr. See ?write_clip for advisories on writing to the clipboard in R.library(gghalves)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(RColorBrewer)
# Funzioni
# Coefficiente di variazione
CV <- function(x){
return(sd(x)/mean(x)*100)
}
# Indice di eterogeneità di Gini
gini.index = function(x){
ni=table(x)
fi=ni/length(x)
fi2=fi^2
J=length(table(x))
gini = 1-sum(fi2)
gini.norm = gini/((J-1)/J)
return(gini.norm)
}2a. Analisi preliminare dei dati
i. Peso neonato (grammi)
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 830, 4930, 4100 | |
| Media, Mediana, Moda | 3284.08, 3300, 3300 | Monomodale |
| Q1, Q3, IQR | 2990, 3620, 630 | |
| Lower, Upper Whisker | 2045, 4565 | |
| Outlier | Numerosi, da 830 a 4930, sia sotto il lower whisker che sopra l’upper | 69 valori, numero basso |
| Dev standard, coeff di variazione | 525.03, 15.98% | |
| Indice di Gini | 0.99 | Molto alto, indice di eterogeneità |
| Indice di Fisher | -0.64 | Negativo, distribuzione concentrata verso valori alti |
| Indice di curtosi | 2.03 | Positivo, distribuzione leptocurtica (più allungata della normale) |
- Il ridotto numero di outlier (69, 2.8% delle osservazioni) può essere mantenuto nel dataframe.
- Distribuzione con concentrazione su valori alti. Dal momento che la variabile si esprime su un numero elevato di modalità , per la rappresentazione possiamo fare uso do una funzione densità .
- La distribuzione del peso del neonato è piuttosto distante da quella di una normale. La cosa è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo.
# Minimo, massimo e range
range(Peso)## [1] 830 4930Peso_min = range(Peso)[1]
Peso_max = range(Peso)[2]
# Media (aritmetica)
Peso_avg = mean(Peso)
# Moda
obs_peso = table(Peso)
sorted_obs_peso = sort(obs_peso, decreasing = TRUE)
print(sorted_obs_peso) ## Peso
## 3300 3500 3100 3200 3400 3250 3150 2900 3600 3550 3000 3050 3350 3380 3180 3450 3700 3800 3280 3290 2950 3140 2880 3170 3220
## 56 45 43 38 35 34 33 31 31 30 29 29 27 27 26 25 25 25 24 24 23 23 22 22 22
## 3750 3030 3340 3440 2850 2940 3080 3260 3330 3520 3530 3950 3190 3370 2750 3240 3310 3480 3580 3650 3820 2800 2980 3110 3230
## 22 21 21 21 20 20 20 20 20 20 20 20 19 19 18 18 18 18 18 18 18 17 17 17 17
## 3270 3540 3640 3680 3850 2740 3900 2700 3320 3620 4000 2990 3040 3060 3070 3360 3420 3570 3630 3670 2920 3430 3510 3560 3590
## 17 17 17 17 17 16 16 15 15 15 15 14 14 14 14 14 14 14 14 14 13 13 13 13 13
## 3780 3840 2820 2830 2910 2970 3090 3160 3720 3760 4050 2680 2890 2960 3020 3120 3130 3460 3730 3770 3860 3870 3910 2500 2840
## 13 13 12 12 12 12 12 12 12 12 12 11 11 11 11 11 11 11 11 11 11 11 11 10 10
## 2930 3010 3410 3660 2600 2720 2780 2810 2860 3210 3470 3610 3710 3830 2650 3690 3740 3790 3920 4140 2550 2560 2710 2870 3390
## 10 10 10 10 9 9 9 9 9 9 9 9 9 9 8 8 8 8 8 8 7 7 7 7 7
## 3940 4060 4150 4200 2400 2450 2580 2590 2670 2760 3490 3810 3960 3990 4100 2100 2520 2620 2660 2770 3890 3930 4020 4130 4250
## 7 7 7 7 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 5 5 5 5 5
## 2430 2440 2630 2690 2730 3880 3970 4030 4090 4120 4160 4220 4240 4260 4330 4400 1280 1750 2000 2050 2220 2260 2300 2320 2340
## 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3
## 2410 2510 2640 2790 3980 4010 4070 4080 4110 4180 4290 4310 4320 4350 930 1170 1500 1780 1970 1980 2040 2060 2150 2160 2200
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2
## 2230 2250 2290 2350 2380 2420 2470 2480 2530 2570 2610 4040 4270 4280 4370 4410 4420 4440 4480 4600 4720 830 900 980 990
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1
## 1140 1180 1190 1285 1300 1340 1370 1390 1410 1430 1450 1550 1560 1580 1600 1615 1620 1690 1720 1730 1770 1800 1840 1850 1890
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1900 1950 1960 2090 2120 2180 2210 2270 2280 2310 2330 2370 2390 2460 2490 2540 2862 4170 4190 4230 4300 4340 4470 4520 4540
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 4550 4560 4580 4620 4650 4680 4690 4700 4760 4810 4900 4930
## 1 1 1 1 1 1 1 1 1 1 1 1# Quartili
Peso_Q1 = quantile(Peso, probs = 0.25)[1]
Peso_median = quantile(Peso, probs = 0.5)[1]
Peso_Q3 = quantile(Peso, probs = 0.75)[1]
Peso_IQR = IQR(Peso)
Peso_Top = Peso_Q3 + 1.5*Peso_IQR
Peso_Low = Peso_Q1 - 1.5*Peso_IQR
# Outlier
Peso_outliers = sort(Peso[Peso < Peso_Low | Peso > Peso_Top])
Peso_outliers## [1] 830 900 930 930 980 990 1140 1170 1170 1180 1190 1280 1280 1280 1285 1300 1340 1370 1390 1410 1430 1450 1500 1500
## [25] 1550 1560 1580 1600 1615 1620 1690 1720 1730 1750 1750 1750 1770 1780 1780 1800 1840 1850 1890 1900 1950 1960 1970 1970
## [49] 1980 1980 2000 2000 2000 2040 2040 4580 4600 4600 4620 4650 4680 4690 4700 4720 4720 4760 4810 4900 4930length(Peso_outliers)## [1] 69# Deviazione standard, coefficiente di variazione e di Gini
Peso_sd = sd(Peso)
Peso_CV = CV(Peso)
gini.index(Peso)## [1] 0.9962164# Indici di forma
skewness(Peso) ## [1] -0.6470308kurtosis(Peso)-3 ## [1] 2.031532# Boxplot
boxplot(Peso)
text(1, Peso_median, labels=round(Peso_median, 2), col="blue", pos=4)
text(1, Peso_Q1, labels=round(Peso_Q1, 2), col="red", pos=4)
text(1, Peso_Q3, labels=round(Peso_Q3, 2), col="red", pos=4)
text(1, Peso_Top, labels=round(Peso_Top, 2), col="orange", pos=4)
text(1, Peso_Low, labels=round(Peso_Low, 2), col="orange", pos=4)
# Distribuzione della variabile e confronto con la normale
ggplot(data = neonati, aes(x = Peso)) +
geom_density(fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Peso neonato", x = "Peso (g)",
y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Peso_avg, sd = Peso_sd),
color = "red", size = 1) +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Shapiro Test normalità distribuzione
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16ii. Età della madre
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 13, 46, 33 | Il minimo sarebbe 0, ma si tratta di un errore (come pure 1). Il range sarebbe 46, ma va corretto per lo stesso motivo. |
| Media, Mediana, Moda | 28.16, 28, 30 | Monomodale |
| Q1, Q3, IQR | 25, 32, 7 | |
| Lower, Upper Whisker | 14.5, 42.5 | |
| Outlier | 0,1,13,14,43,44,45,46 | 13 osservazioni, numero molto basso |
| Dev standard, coeff di variazione | 5.27, 18.72% | |
| Indice di Gini | 0.97 | Molto vicino a 1, massima eterogeneità |
| Indice di Fisher | 0.04 | Molto basso, distribuzione virtualmente simmetrica |
| Indice di curtosi | 0.38 | Molto basso, distribuzione virtualmente mesocurtica |
| Coeff. Bravais-Pearson età /peso | -0.022 | Leggera correlazione negativa (maggiore l’età della madre, minore il peso del bambino) |
- Il ridotto numero di outlier (13, 0.5% delle osservazioni) può essere mantenuto nel dataframe. Due outlier presentano però valori di età non plausibili (0,1): vanno corretti e li sostituiremo con il valore mediano.
- Dal momento che la variabile si esprime su un numero discreto di modalità , per la rappresentazione meglio usare istogrammi.
- La distribuzione è vicina a quella di una normale ma la cosa non è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo. Verosimilmente, il test fallisce a causa delle seppur ridotte deviazioni dalla curva normale.
- Nell’analisi della relazione tra età madre e peso del nascituro non si nota un trend particolare dallo scatterplot, cosa confermata dal basso valore del coefficiente di correlazione lineare Bravais-Pearson.
# Minimo, massimo e range
range(Anni.madre)## [1] 0 46Anni_min = range(Anni.madre)[1]
Anni_max = range(Anni.madre)[2]
# Media (aritmetica)
Anni_avg = mean(Anni.madre)
# Moda
obs_anni = table(Anni.madre)
sorted_obs_anni = sort(obs_anni, decreasing = TRUE)
print(sorted_obs_anni) ## Anni.madre
## 30 27 26 25 29 28 32 31 24 23 33 22 34 21 20 35 36 19 37 38 39 18 40 17 16 41 42 15 44 14 43
## 200 197 184 180 174 172 159 147 131 115 110 100 96 74 66 66 64 45 41 38 27 24 19 18 13 13 8 6 4 2 2
## 0 1 13 45 46
## 1 1 1 1 1# Quartili
Anni_Q1 = quantile(Anni.madre, probs = 0.25)[1]
Anni_median = quantile(Anni.madre, probs = 0.5)[1]
Anni_Q3 = quantile(Anni.madre, probs = 0.75)[1]
Anni_IQR = IQR(Anni.madre)
Anni_Top = Anni_Q3 + 1.5*Anni_IQR
Anni_Low = Anni_Q1 - 1.5*Anni_IQR
# Outlier
Anni_outliers = sort(Anni.madre[Anni.madre < Anni_Low | Anni.madre > Anni_Top])
Anni_outliers## [1] 0 1 13 14 14 43 43 44 44 44 44 45 46length(Anni_outliers)## [1] 13# Deviazione standard, coefficiente di variazione e di Gini
Anni_sd = sd(Anni.madre)
Anni_CV = CV(Anni.madre)
gini.index(Anni.madre)## [1] 0.9727146# Indici di forma
skewness(Anni.madre) ## [1] 0.0428115kurtosis(Anni.madre)-3 ## [1] 0.3804165# Boxplot
boxplot(Anni.madre)
text(1, Anni_median, labels=round(Anni_median, 2), col="blue", pos=4)
text(1, Anni_Q1, labels=round(Anni_Q1, 2), col="red", pos=4)
text(1, Anni_Q3, labels=round(Anni_Q3, 2), col="red", pos=4)
text(1, Anni_Top, labels=round(Anni_Top, 2), col="orange", pos=4)
text(1, Anni_Low, labels=round(Anni_Low, 2), col="orange", pos=4)
# Distribuzione variabile e confronto con la normale
ggplot(data = neonati, aes(x = Anni.madre)) +
geom_histogram(aes(y = after_stat(density)), fill = "lightblue", color = "black") +
stat_bin(aes(y = after_stat(density), label = round(after_stat(density), 2)),
geom = "text", color = "black", vjust = -0.5, size = 3) +
scale_y_continuous(limits = c(0, 0.1)) +
labs(title = "Anni Madre", x = "Anni madre", y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Anni_avg, sd = Anni_sd),
color = "red", size = 1) +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Shapiro Test normalità distribuzione
shapiro.test(Anni.madre)##
## Shapiro-Wilk normality test
##
## data: Anni.madre
## W = 0.99308, p-value = 1.639e-09# Analisi relazione età - peso neonato
plot(Anni.madre,Peso, pch=20)
pearson_rho_Anni_Peso = cov(Anni.madre, Peso)/(Anni_sd*Peso_sd)
pearson_rho_Anni_Peso## [1] -0.02247017# Correzione outlier
neonati$Anni.madre[neonati$Anni.madre < 13] = Anni_medianiii. Numero gravidanze
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 0, 12, 12 | |
| Media, Mediana, Moda | 0.98, 1, 0 | Monomodale |
| Q1, Q3, IQR | 0, 1, 1 | |
| Lower, Upper Whisker | -1.5, 2.5 | |
| Outlier | 3,4,5,6,7,8,9,10,11,12 | 246 valori, numero significativo |
| Dev standard, coeff di variazione | 1.28, 130.5% | |
| Indice di Gini | 0.73 | Piuttosto alto, indice di eterogeneità |
| Indice di Fisher | 2.51 | Positivo, distribuzione concentrata verso valori bassi |
| Indice di curtosi | 10.98 | Positivo, distribuzione leptocurtica (più allungata della normale) |
| Coeff. Bravais-Pearson numero gravidanze/peso | 0.002 | Correlazione positiva virtualmente nulla. |
È presente un significativo numero di outlier (246, 9.8% del totale osservazioni). Nonostante ciò, meglio non fare nulla e lasciarli nel dataframe per non perdere informazioni.
Distribuzione concentrata su bassi valori. Dal momento che la variabile si esprime su un numero discreto di modalità , per la rappresentazione meglio usare istogrammi.
La distribuzione del numero di gravidanze è decisamente distante da quella di una normale. La cosa è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo.
Nell’analisi della relazione tra numero delle gravidanze e peso del nascituro non si nota un trend particolare dallo scatterplot, cosa confermata dal basso valore del coefficiente di correlazione lineare Bravais-Pearson.
# Minimo, massimo e range
range(N.gravidanze)## [1] 0 12Grav_min = range(N.gravidanze)[1]
Grav_max = range(N.gravidanze)[2]
# Media (aritmetica)
Grav_avg = mean(N.gravidanze)
# Moda
obs_grav = table(N.gravidanze)
sorted_obs_grav = sort(obs_grav, decreasing = TRUE)
print(sorted_obs_grav) ## N.gravidanze
## 0 1 2 3 4 5 6 8 10 9 7 11 12
## 1096 818 340 150 48 21 11 8 3 2 1 1 1# Quartili
Grav_Q1 = quantile(N.gravidanze, probs = 0.25)[1]
Grav_median = quantile(N.gravidanze, probs = 0.5)[1]
Grav_Q3 = quantile(N.gravidanze, probs = 0.75)[1]
Grav_IQR = IQR(N.gravidanze)
Grav_Top = Grav_Q3 + 1.5*Grav_IQR
Grav_Low = Grav_Q1 - 1.5*Grav_IQR
# Outlier
Grav_outliers = sort(N.gravidanze[N.gravidanze < Grav_Low | N.gravidanze > Grav_Top])
Grav_outliers## [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [41] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [81] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [121] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
## [161] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5
## [201] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 7 8 8 8 8 8 8 8 8 9
## [241] 9 10 10 10 11 12length(Grav_outliers)## [1] 246# Deviazione standard, coefficiente di variazione e di Gini
Grav_sd = sd(N.gravidanze)
Grav_CV = CV(N.gravidanze)
gini.index(N.gravidanze)## [1] 0.7346931# Indici di forma
skewness(N.gravidanze) ## [1] 2.514254kurtosis(N.gravidanze)-3 ## [1] 10.98941# Boxplot
boxplot(N.gravidanze)
text(1, Grav_median, labels=round(Grav_median, 2), col="blue", pos=4)
text(1, Grav_Q1, labels=round(Grav_Q1, 2), col="red", pos=4)
text(1, Grav_Q3, labels=round(Grav_Q3, 2), col="red", pos=4)
text(1, Grav_Top, labels=round(Grav_Top, 2), col="orange", pos=4)
text(1, Grav_Low, labels=round(Grav_Low, 2), col="orange", pos=4)
# Distribuzione della variabile e confronto con la normale
ggplot(data = neonati, aes(x = N.gravidanze)) +
geom_histogram(aes(y = ..density..), binwidth = 1, fill = "lightblue", color = "black") +
stat_bin(aes(y = ..density.., label = round(..density.., 2)), geom = "text",
binwidth = 1, vjust = -0.5, color = "black", size = 3) +
scale_y_continuous(limits = c(0, 0.5)) +
labs(title = "Numero di gravidanze", x = "Gravidanze", y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Grav_avg, sd = Grav_sd),
color = "red", size = 1) +
theme_minimal()## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Shapiro Test normalità distribuzione
shapiro.test(N.gravidanze)##
## Shapiro-Wilk normality test
##
## data: N.gravidanze
## W = 0.72125, p-value < 2.2e-16# Analisi relazione numero gravidanze - peso neonato
plot(N.gravidanze,Peso, pch=20)
pearson_rho_Grav_Peso = cov(N.gravidanze, Peso)/(Grav_sd*Peso_sd)
pearson_rho_Grav_Peso## [1] 0.0024073iv. Fumatrici
| Fumatrice | Ospedale1 | Ospedale2 | Ospedale3 | TOTAL |
|---|---|---|---|---|
| NO | 31.16% | 32.52% | 32.16% | 95.84% |
| SI | 1.48% | 1.44% | 1.24% | 4.16% |
| TOTAL | 32.64% | 33.96% | 33.40% | 100.00% |
La grandissima maggioranza delle gestanti (96%) non fuma. Non sembrano esserci differenze significative tra i gruppi nei vari ospedali.
Per saggiare l’effetto del fumo sul peso del neonato, per variabili categoriali come questa possiamo fare uso del t-test. Il pvalue fornito è 0.30, pertanto non possiano rifiutare l’ipotesi nulla che le medie del peso nei due gruppi siano sostanzialmente uguali (3286.15 per le non fumatrici, 3236.34 per le fumatrici).
Distr_fum = table(Fumatrici, Ospedale)/dim(neonati)[1]*100
Distr_fum_total = addmargins(Distr_fum, margin = c(1,2))
Distr_fum_total## Ospedale
## Fumatrici osp1 osp2 osp3 Sum
## 0 31.16 32.52 32.16 95.84
## 1 1.48 1.44 1.24 4.16
## Sum 32.64 33.96 33.40 100.00# T-test per saggiare uguaglianza della media del peso nei due gruppi
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.346v. Settimane di gestazione
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 25, 43, 18 | |
| Media, Mediana, Moda | 38.98, 39, 40 | Monomodale |
| Q1, Q3, IQR | 38, 40, 2 | |
| Lower, Upper Whisker | 35, 43 | |
| Outlier | 25,26,27,28,29,30,31,32,33,34 | 67 valori, numero basso |
| Deviazione standard | 1.86 | |
| Coefficiente di variazione | 4.79% | |
| Indice di Gini | 0.84 | Piuttosto alto, indice di eterogeneità |
| Indice di Fisher | -2.06 | Negativo, distribuzione concentrata verso valori alti |
| Indice di curtosi | 8.25 | Positivo, distribuzione leptocurtica (più allungata della normale) |
| Coeff. Bravais-Pearson età /peso | 0.59 | Correlazione positiva (maggiore il numero di settimane di gestazione, maggiore il peso del bambino) |
- Il ridotto numero di outlier (67, 2.7% delle osservazioni) può essere mantenuto nel dataframe.
- Distribuzione concentrata su bassi valori. Dal momento che la variabile si esprime su un numero discreto di modalità , per la rappresentazione meglio usare istogrammi.
- La distribuzione delle settimane di gestazione non è affatto vicina a quella di una normale. La cosa è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo.
- Nell’analisi della relazione tra settimane di gestazione e peso del nascituro si nota un trend di correlazione positiva nello scatterplot, confermata dal valore positivo piuttosto elevato del coefficiente di correlazione lineare Bravais-Pearson.
# Minimo, massimo e range
range(Gestazione)## [1] 25 43Gest_min = range(Gestazione)[1]
Gest_max = range(Gestazione)[2]
# Media (aritmetica)
Gest_avg = mean(Gestazione)
# Moda
obs_gest = table(Gestazione)
sorted_obs_gest = sort(obs_gest, decreasing = TRUE)
print(sorted_obs_gest) ## Gestazione
## 40 39 38 41 37 36 42 35 33 34 32 31 30 28 29 27 43 25 26
## 741 581 437 329 192 62 56 33 18 16 9 8 5 4 3 2 2 1 1# Quartili
Gest_Q1 = quantile(Gestazione, probs = 0.25)[1]
Gest_median = quantile(Gestazione, probs = 0.5)[1]
Gest_Q3 = quantile(Gestazione, probs = 0.75)[1]
Gest_IQR = IQR(Gestazione)
Gest_Top = Gest_Q3 + 1.5*Gest_IQR
Gest_Low = Gest_Q1 - 1.5*Gest_IQR
# Outlier
Gest_outliers = sort(Gestazione[Gestazione < Gest_Low | Gestazione > Gest_Top])
Gest_outliers## [1] 25 26 27 27 28 28 28 28 29 29 29 30 30 30 30 30 31 31 31 31 31 31 31 31 32 32 32 32 32 32 32 32 32 33 33 33 33 33 33 33
## [41] 33 33 33 33 33 33 33 33 33 33 33 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34 34length(Gest_outliers)## [1] 67# Deviazione standard, coefficiente di variazione e di Gini
Gest_sd = sd(Gestazione)
Gest_CV = CV(Gestazione)
gini.index(Gestazione)## [1] 0.8475571# Indici di forma
skewness(Gestazione) ## [1] -2.065313kurtosis(Gestazione)-3 ## [1] 8.25815# Boxplot
boxplot(Gestazione)
text(1, Gest_median, labels=round(Gest_median, 2), col="blue", pos=4)
text(1, Gest_Q1, labels=round(Gest_Q1, 2), col="red", pos=4)
text(1, Gest_Q3, labels=round(Gest_Q3, 2), col="red", pos=4)
text(1, Gest_Top, labels=round(Gest_Top, 2), col="orange", pos=4)
text(1, Gest_Low, labels=round(Gest_Low, 2), col="orange", pos=4)
# Distribuzione della variabile e confronto con la normale
ggplot(data = neonati, aes(x = Gestazione)) +
geom_histogram(aes(y = ..density..), binwidth = 1, fill = "lightblue", color = "black") +
stat_bin(aes(y = ..density.., label = round(..density.., 2)), geom = "text",
binwidth = 1, vjust = -0.5, color = "black", size = 3) +
scale_y_continuous(limits = c(0, 0.4)) +
labs(title = "Settimane di gestazione", x = "Settimane", y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Gest_avg, sd = Gest_sd),
color = "red", size = 1) +
theme_minimal()
# Shapiro Test normalità distribuzione
shapiro.test(Gestazione)##
## Shapiro-Wilk normality test
##
## data: Gestazione
## W = 0.83328, p-value < 2.2e-16# Analisi relazione numero gravidanze - peso neonato
plot(Gestazione,Peso, pch=20)
pearson_rho_Gest_Peso = cov(Gestazione, Peso)/(Gest_sd*Peso_sd)
pearson_rho_Gest_Peso## [1] 0.5917687vi. Lunghezza neonato (cm)
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 310, 565, 255 | |
| Media, Mediana, Moda | 494.69, 500, 500 | Monomodale |
| Q1, Q3, IQR | 480, 510, 30 | |
| Lower, Upper Whisker | 435, 555 | |
| Outlier | Numerosi, da 310 a 565, sia sotto il lower whisker che sopra l’upper | 59 valori, numero basso |
| Dev standard, coeff. di variazione | 26.31, 5.32% | |
| Indice di Gini | 0.94 | Molto alto, indice di eterogeneità |
| Indice di Fisher | -1.51 | Negativo, distribuzione concentrata verso valori alti |
| Indice di curtosi | 6.48 | Positivo, distribuzione leptocurtica (più allungata della normale) |
| Coeff. Bravais-Pearson lunghezza/peso | 0.79 | Correlazione positiva (maggiore la lunghezza, maggiore il peso del neonato) |
- Il ridotto numero di outlier (69, 2.3% delle osservazioni) può essere mantenuto nel dataframe.
- Distribuzione con concentrazione su valori alti. Dal momento che la variabile si esprime su un numero elevato di modalità , per la rappresentazione possiamo fare uso do una funzione densità .
- La distribuzione della lunghezza del neonato è decisamente distante da quella di una normale. La cosa è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo.
- Nell’analisi della relazione tra lunghezza e peso del nascituro si nota un trend di correlazione positiva nello scatterplot, confermata dal valore positivo elevato del coefficiente di correlazione lineare Bravais-Pearson.
# Minimo, massimo e range
range(Lunghezza)## [1] 310 565Lun_min = range(Lunghezza)[1]
Lun_max = range(Lunghezza)[2]
# Media (aritmetica)
Lun_avg = mean(Lunghezza)
# Moda
obs_lun = table(Lunghezza)
sorted_obs_lun = sort(obs_lun, decreasing = TRUE)
print(sorted_obs_lun) ## Lunghezza
## 500 510 490 480 520 495 470 485 505 515 530 460 475 525 450 465 540 535 550 455 440 445 545 430 410 420 390 400 370 380 405
## 396 273 270 195 172 150 126 122 117 94 89 88 66 57 43 32 32 24 23 16 13 13 12 10 8 8 5 4 3 3 3
## 555 355 360 435 446 498 560 310 315 320 325 340 345 385 425 473 492 494 497 502 504 514 518 523 565
## 3 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1# Quartili
Lun_Q1 = quantile(Lunghezza, probs = 0.25)[1]
Lun_median = quantile(Lunghezza, probs = 0.5)[1]
Lun_Q3 = quantile(Lunghezza, probs = 0.75)[1]
Lun_IQR = IQR(Lunghezza)
Lun_Top = Lun_Q3 + 1.5*Lun_IQR
Lun_Low = Lun_Q1 - 1.5*Lun_IQR
# Outlier
Lun_outliers = sort(Lunghezza[Lunghezza < Lun_Low | Lunghezza > Lun_Top])
Lun_outliers## [1] 310 315 320 325 340 345 355 355 360 360 370 370 370 380 380 380 385 390 390 390 390 390 400 400 400 400 405 405 405 410
## [31] 410 410 410 410 410 410 410 420 420 420 420 420 420 420 420 425 430 430 430 430 430 430 430 430 430 430 560 560 565length(Lun_outliers)## [1] 59# Deviazione standard, coefficiente di variazione e di Gini
Lun_sd = sd(Lunghezza)
Lun_CV = CV(Lunghezza)
gini.index(Lunghezza)## [1] 0.9404748# Indici di forma
skewness(Lunghezza) ## [1] -1.514699kurtosis(Lunghezza)-3 ## [1] 6.487174# Boxplot
boxplot(Lunghezza)
text(1, Lun_median, labels=round(Lun_median, 2), col="blue", pos=4)
text(1, Lun_Q1, labels=round(Lun_Q1, 2), col="red", pos=4)
text(1, Lun_Q3, labels=round(Lun_Q3, 2), col="red", pos=4)
text(1, Lun_Top, labels=round(Lun_Top, 2), col="orange", pos=4)
text(1, Lun_Low, labels=round(Lun_Low, 2), col="orange", pos=4)
# Distribuzione della variabile e confronto con la normale
ggplot(data = neonati, aes(x = Lunghezza)) +
geom_density(fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Lunghezza neonato", x = "Lunghezza (cm)",
y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Lun_avg, sd = Lun_sd),
color = "red", size = 1) +
theme_minimal()
# Shapiro Test normalità distribuzione
shapiro.test(Lunghezza)##
## Shapiro-Wilk normality test
##
## data: Lunghezza
## W = 0.90941, p-value < 2.2e-16# Analisi relazione lunghezza - peso neonato
plot(Lunghezza,Peso, pch=20)
pearson_rho_Lun_Peso = cov(Lunghezza, Peso)/(Lun_sd*Peso_sd)
pearson_rho_Lun_Peso## [1] 0.7960368vii. Dimensione cranio (cm)
| Indice | Valore | Commenti |
|---|---|---|
| Minimo, Massimo, Range | 235, 390, 155 | |
| Media, Mediana, Moda | 340.02, 340, 340 | Monomodale |
| Q1, Q3, IQR | 330, 350, 20 | |
| Lower, Upper Whisker | 300, 380 | |
| Outlier | Numerosi, da 235 a 390, sia sotto il lower whisker che sopra l’upper | 48 valori, numero basso |
| Dev standard, coeff. di variazione | 16.42, 4.83% | |
| Indice di Gini | 0.97 | Molto alto, indice di eterogeneità |
| Indice di Fisher | -0.78 | Negativo, distribuzione concentrata verso valori alti |
| Indice di curtosi | 2.94 | Positivo, distribuzione leptocurtica (più allungata della normale) |
| Coeff. Bravais-Pearson diametro crano/peso | 0.70 | Correlazione positiva (maggiore il diametro del cranio, maggiore il peso del neonato) |
- Il ridotto numero di outlier (48, 1.9% delle osservazioni) può essere mantenuto nel dataframe.
- Distribuzione con concentrazione su valori alti. Dal momento che la variabile si esprime su un numero elevato di modalità , per la rappresentazione possiamo fare uso do una funzione densità .
- La distribuzione del dimatetro del cranio del neonato è piuttosto distante da quella di una normale. La cosa è confermata dal test Shapiro-Wilk, che presenta un p-value praticamente nullo.
- Nell’analisi della relazione tra diametro cranio e peso del nascituro si nota un trend di correlazione positiva nello scatterplot, confermata dal valore positivo elevato del coefficiente di correlazione lineare Bravais-Pearson.
# Minimo, massimo e range
range(Cranio)## [1] 235 390Cra_min = range(Cranio)[1]
Cra_max = range(Cranio)[2]
# Media (aritmetica)
Cra_avg = mean(Cranio)
# Moda
obs_cra = table(Cranio)
sorted_obs_cra = sort(obs_cra, decreasing = TRUE)
print(sorted_obs_cra) ## Cranio
## 340 350 345 335 330 325 355 360 320 342 343 332 336 346 344 338 347 352 334 365 348 315 333 357 326 370 337 353 310 323 328
## 238 169 159 155 148 101 101 94 68 53 48 46 45 45 43 42 42 41 38 36 34 33 30 30 29 29 28 27 26 26 25
## 354 358 322 356 327 339 362 331 305 349 341 300 324 364 312 368 318 351 363 361 329 375 313 366 367 380 317 321 359 372 373
## 25 25 24 24 19 19 19 17 16 16 15 14 14 12 11 11 10 10 10 9 8 8 7 7 6 6 5 5 5 5 5
## 316 390 290 295 304 307 308 309 319 273 276 277 280 292 298 303 314 369 371 374 376 378 382 235 245 253 254 265 266 267 270
## 4 4 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1
## 272 274 275 278 285 287 289 293 294 297 299 301 306 379 381 383 384 385 386
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1# Quartili
Cra_Q1 = quantile(Cranio, probs = 0.25)[1]
Cra_median = quantile(Cranio, probs = 0.5)[1]
Cra_Q3 = quantile(Cranio, probs = 0.75)[1]
Cra_IQR = IQR(Cranio)
Cra_Top = Cra_Q3 + 1.5*Cra_IQR
Cra_Low = Cra_Q1 - 1.5*Cra_IQR
# Outlier
Cra_outliers = sort(Cranio[Cranio < Cra_Low | Cranio > Cra_Top])
Cra_outliers## [1] 235 245 253 254 265 266 267 270 272 273 273 274 275 276 276 277 277 278 280 280 285 287 289 290 290 290 292 292 293 294
## [31] 295 295 295 297 298 298 299 381 382 382 383 384 385 386 390 390 390 390length(Cra_outliers)## [1] 48# Deviazione standard, coefficiente di variazione e di Gini
Cra_sd = sd(Cranio)
Cra_CV = CV(Cranio)
gini.index(Cranio)## [1] 0.9723873# Indici di forma
skewness(Cranio) ## [1] -0.7850527kurtosis(Cranio)-3 ## [1] 2.946206# Boxplot
boxplot(Cranio)
text(1, Cra_median, labels=round(Cra_median, 2), col="blue", pos=4)
text(1, Cra_Q1, labels=round(Cra_Q1, 2), col="red", pos=4)
text(1, Cra_Q3, labels=round(Cra_Q3, 2), col="red", pos=4)
text(1, Cra_Top, labels=round(Cra_Top, 2), col="orange", pos=4)
text(1, Cra_Low, labels=round(Cra_Low, 2), col="orange", pos=4)
# Distribuzione della variabile e confronto con la normale
ggplot(data = neonati, aes(x = Cranio)) +
geom_density(fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Cranio neonato", x = "Diametro (cm)",
y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Cra_avg, sd = Cra_sd), color = "red", size = 1) +
theme_minimal()
# Shapiro Test normalità distribuzione
shapiro.test(Cranio)##
## Shapiro-Wilk normality test
##
## data: Cranio
## W = 0.96357, p-value < 2.2e-16# Analisi relazione cranio - peso neonato
plot(Cranio,Peso, pch=20)
pearson_rho_Cra_Peso = cov(Cranio, Peso)/(Cra_sd*Peso_sd)
pearson_rho_Cra_Peso## [1] 0.7048015viii. Tipo parto
| Tipo parto | Ospedale1 | Ospedale2 | Ospedale3 | TOTAL |
|---|---|---|---|---|
| Cesareo | 9.68% | 10.16% | 9.28% | 29.12% |
| Naturale | 22.96% | 23.80% | 24.12% | 70.88% |
| TOTAL | 32.64% | 33.96% | 33.40% | 100.00% |
- La gran parte delle gestanti (71%) partorisce naturalmente. L’incidenza del parto naturale aumenta progressivamente da Ospedale 1 a 2 e infine 3.
- Per saggiare l’effetto del tipo di parto sul peso del neonato, per variabili categoriali come questa possiamo fare uso del t-test. Il pvalue fornito è 0.89, pertanto non possiano rifiutare l’ipotesi nulla che le medie del peso nei due gruppi siano sostanzialmente uguali (3282.04 per i parti cesarei, 3284.91 per quelli naturali).
Distr_parto = table(Tipo.parto, Ospedale)/dim(neonati)[1]*100
Distr_parto_total = addmargins(Distr_parto, margin = c(1,2))
Distr_parto_total## Ospedale
## Tipo.parto osp1 osp2 osp3 Sum
## Ces 9.68 10.16 9.28 29.12
## Nat 22.96 23.80 24.12 70.88
## Sum 32.64 33.96 33.40 100.00# T-test per saggiare uguaglianza della media del peso nei due gruppi
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.916ix. Ospedale
| Ospedale 1 | Ospedale 2 | Ospedale 3 | TOTAL |
|---|---|---|---|
| 816 | 849 | 835 | 2500 |
| 32.64% | 33.96% | 33.40% | 100.00% |
La ripartizione delle gestanti è quasi equidistribuita sui tre ospedali.
Per saggiare l’effetto dell’ospedale sul peso del neonato, per variabili categoriali a più di due valori come questa possiamo fare uso del test ANOVA. Il pvalue fornito è elevato (0.18), pertanto non rifiutiamo l’ipotesi nulla che le medie del peso nei tre gruppi siano uguali. L’ospedale quindi non ha effetto significativo sul peso del neonato.
table(Ospedale)## Ospedale
## osp1 osp2 osp3
## 816 849 835# ANOVA test per saggiare uguaglianza della media del peso nei tre gruppi
Osp_test = aov(Peso ~ Ospedale, data = neonati)
summary(Osp_test)## Df Sum Sq Mean Sq F value Pr(>F)
## Ospedale 2 936237 468118 1.699 0.183
## Residuals 2497 687952305 275512x. Sesso nascituro
| Sesso | Ospedale1 | Ospedale2 | Ospedale3 | TOTAL |
|---|---|---|---|---|
| F | 16.36% | 17.40% | 16.48% | 50.24% |
| M | 16.28% | 16.56% | 16.92% | 49.76% |
| TOTAL | 32.64% | 33.96% | 33.40% | 100.00% |
La ripartizione F/M è praticamente 50/50. Ospedale 3 presenta in effetti più maschi che femmine, diversamente da quanto accade nelle altre due strutture.
Per saggiare l’effetto del sesso sul peso del neonato, per variabili categoriali come questa possiamo fare uso del t-test. Il pvalue fornito è praticamente nullo, pertanto rifiutiamo l’ipotesi nulla che le medie del peso nei due gruppi siano sostanzialmente uguali (3161.13 per le femmine, 3408.21 per i maschi).
Distr_sesso = table(Sesso, Ospedale)/dim(neonati)[1]*100
Distr_sesso_total = addmargins(Distr_sesso, margin = c(1,2))
Distr_sesso_total## Ospedale
## Sesso osp1 osp2 osp3 Sum
## F 16.36 17.40 16.48 50.24
## M 16.28 16.56 16.92 49.76
## Sum 32.64 33.96 33.40 100.00# T-test per saggiare uguaglianza della media del peso nei due gruppi
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.2152b. Creazione del modello di regressione
Step1 - Normalità della variabile risposta
Come visto in precedenza, la variabile è leptocurtica (+2.03) e asimmetrica negativa (-0.64), quindi non normale, come conferma lo Shapiro-Wilk test). Ciononostante, la differenza tra la distribuzione e quella normale non appare così forte e procederemo ugualmente alla creazione del modello di regressione lineare.
ggplot(data = neonati, aes(x = Peso)) +
geom_density(fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Peso neonato", x = "Peso (g)", y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = Peso_avg, sd = Peso_sd),
color = "red", size = 1) +
theme_minimal()
skewness(Peso)## [1] -0.6470308kurtosis(Peso)-3## [1] 2.031532shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Step2 - Matrice di correlazione
Creaimo la matrice per verificare quali siano le variabili numeriche più significative.
neonati_numeric = neonati[sapply(neonati, is.numeric)]
matrix_corr = round(cor(neonati_numeric, use = "complete.obs"),2)
print(matrix_corr)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## Anni.madre 1.00 0.38 0.01 -0.13 -0.02 -0.06 0.02
## N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06 0.04
## Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02 -0.01
## Gestazione -0.13 -0.10 0.03 1.00 0.59 0.62 0.46
## Peso -0.02 0.00 -0.02 0.59 1.00 0.80 0.70
## Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00 0.60
## Cranio 0.02 0.04 -0.01 0.46 0.70 0.60 1.00Correlazione negativa moderata con l’età della madre (-0.02).
Correlazione quasi nulla con il numero delle gravidanze (0.00).
Correlazione negativa moderata con il fumo (-0.02).
Correlazione positiva significativa con le settimane di gestazione (+0.59).
Correlazione positiva significativa con la lunghezza (+0.80).
Correlazione positiva significativa con il diametro cranio (+0.70).
Per le variabili qualitative bisogna fare uso di strumenti diversi, come i boxplot. Tale analisi è già stata effettuata in precedenza e ci porta a dire che:
- Fumatrici - Il peso medio dei neonati non è significativamente differente nei due gruppi. Dal punto di vista clinico, il fumo della madre è però riconosciuto come un fattore che può influenzare negativamente la salute del neonato, pertanto è bene includerlo in prima battuta nel modello di regressione per rispettare la teoria scientifica.
- Tipo Parto - Il peso medio dei neonati non è significativamente differente nei due gruppi. Dal punto di vista clinico, la modalità del parto è spesso considerata un fattore potenzialmente rilevante per il benessere del neonato, quindi è bene includerla in prima battuta nel modello di regressione.
- Ospedale - Il peso medio dei neonati non è significativamente differenti nei tre gruppi. Anche se il test non rileva differenze significative, potrebbe comunque esserci un piccolo effetto che diventa più evidente in un modello di regressione multipla, quindi è bene includerlo in prima battuta.
- Sesso - Come anticipato, il peso medio dei neonati è significativamente differente nei due gruppi e la variabile va assolutamente inclusa nel modello di regressione.
Step3 - Modello con tutte le variabili
mod1 = lm(Peso ~ . , data = neonati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = neonati)
##
## 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-16Considerazioni:
Le variabili età della madre e fumo non sono statisticamente significative.
Le variabili settimane di gestazione, lunghezza, cranio e sesso sono statisticamente significative.
Le variabili numero di gravidanze, tipo di parto e ospedale sono statisticamente abbastanza significative.
Adjusted R-squared: 72.78%
Step4 - Aggiornamenti successivi
Cominciamo con il rimuovere le variabili non statisticamente significative, quali età della madre e fumo. Adjusted R-squared resta pari 72.78%. Si tratta però di variabili cui onestamente è difficile rinunciare: la prima è una classica variabile di controllo, la seconda è sospetta di avere un ruolo nel peso del neonato.
mod2 = update(mod1, ~. -Anni.madre -Fumatrici)
summary(mod2) # Adjusted R-squared 72.78% ##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati)
##
## 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-16car::vif(mod2) # tutti minori di 5, no multicollinearità ## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.025244 1 1.012544
## Gestazione 1.670237 1 1.292376
## Lunghezza 2.081921 1 1.442886
## Cranio 1.626507 1 1.275346
## Tipo.parto 1.003872 1 1.001934
## Ospedale 1.003038 2 1.000759
## Sesso 1.040337 1 1.019969Proviamo a verificare un’eventuale presenza di interazioni tra le variabili, partendo da mod1.
mod3 = update(mod1, ~. +N.gravidanze*Gestazione)
summary(mod3)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + N.gravidanze:Gestazione,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1122.68 -181.66 -15.26 160.38 2611.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6647.2095 163.7446 -40.595 < 2e-16 ***
## Anni.madre 0.7909 1.1463 0.690 0.4903
## N.gravidanze -62.9547 69.9781 -0.900 0.3684
## Fumatrici -30.6939 27.5435 -1.114 0.2652
## Gestazione 30.2415 4.3793 6.906 6.32e-12 ***
## Lunghezza 10.2943 0.3007 34.235 < 2e-16 ***
## Cranio 10.4784 0.4261 24.593 < 2e-16 ***
## Tipo.partoNat 29.4362 12.0846 2.436 0.0149 *
## Ospedaleosp2 -11.2707 13.4385 -0.839 0.4017
## Ospedaleosp3 27.7924 13.4999 2.059 0.0396 *
## SessoM 77.4792 11.1778 6.932 5.28e-12 ***
## N.gravidanze:Gestazione 1.9172 1.8001 1.065 0.2869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2488 degrees of freedom
## Multiple R-squared: 0.729, Adjusted R-squared: 0.7278
## F-statistic: 608.4 on 11 and 2488 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod4 = update(mod1, ~. +Sesso*Gestazione)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Gestazione:Sesso,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1117.93 -180.40 -14.42 161.62 2610.43
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6611.7541 171.1164 -38.639 < 2e-16 ***
## Anni.madre 0.8085 1.1462 0.705 0.4806
## N.gravidanze 11.3550 4.6661 2.434 0.0150 *
## Fumatrici -31.1488 27.5469 -1.131 0.2583
## Gestazione 29.2584 4.5922 6.371 2.23e-10 ***
## Lunghezza 10.2994 0.3007 34.254 < 2e-16 ***
## Cranio 10.4733 0.4260 24.585 < 2e-16 ***
## Tipo.partoNat 29.8886 12.0870 2.473 0.0135 *
## Ospedaleosp2 -10.8430 13.4403 -0.807 0.4199
## Ospedaleosp3 28.6393 13.5021 2.121 0.0340 *
## SessoM -222.2860 234.4661 -0.948 0.3432
## Gestazione:SessoM 7.6834 6.0015 1.280 0.2006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2488 degrees of freedom
## Multiple R-squared: 0.7291, Adjusted R-squared: 0.7279
## F-statistic: 608.6 on 11 and 2488 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod5 = update(mod1, ~. +Tipo.parto*Gestazione)
summary(mod5)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Gestazione:Tipo.parto,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1128.60 -180.42 -14.09 162.16 2613.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6496.8679 246.9721 -26.306 < 2e-16 ***
## Anni.madre 0.7404 1.1473 0.645 0.5188
## N.gravidanze 11.4204 4.6661 2.448 0.0145 *
## Fumatrici -30.6589 27.5408 -1.113 0.2657
## Gestazione 26.6399 6.2925 4.234 2.38e-05 ***
## Lunghezza 10.2853 0.3008 34.195 < 2e-16 ***
## Cranio 10.4660 0.4261 24.565 < 2e-16 ***
## Tipo.partoNat -281.4365 264.4989 -1.064 0.2874
## Ospedaleosp2 -11.3113 13.4379 -0.842 0.4000
## Ospedaleosp3 27.6444 13.5017 2.047 0.0407 *
## SessoM 78.0996 11.1869 6.981 3.74e-12 ***
## Gestazione:Tipo.partoNat 7.9711 6.7735 1.177 0.2394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2488 degrees of freedom
## Multiple R-squared: 0.729, Adjusted R-squared: 0.7278
## F-statistic: 608.5 on 11 and 2488 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod6 = update(mod1, ~. +Ospedale*Gestazione)
summary(mod6)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Gestazione:Ospedale,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1122.20 -181.07 -14.15 161.10 2611.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6724.6815 208.0044 -32.330 < 2e-16 ***
## Anni.madre 0.7947 1.1472 0.693 0.4885
## N.gravidanze 11.4442 4.6730 2.449 0.0144 *
## Fumatrici -30.2271 27.5536 -1.097 0.2727
## Gestazione 32.2616 5.6125 5.748 1.01e-08 ***
## Lunghezza 10.2951 0.3008 34.220 < 2e-16 ***
## Cranio 10.4720 0.4265 24.552 < 2e-16 ***
## Tipo.partoNat 29.5555 12.0954 2.444 0.0146 *
## Ospedaleosp2 -63.7453 285.9050 -0.223 0.8236
## Ospedaleosp3 42.5179 273.8772 0.155 0.8766
## SessoM 77.5125 11.1867 6.929 5.38e-12 ***
## Gestazione:Ospedaleosp2 1.3480 7.3303 0.184 0.8541
## Gestazione:Ospedaleosp3 -0.3690 7.0168 -0.053 0.9581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7276
## F-statistic: 557.2 on 12 and 2487 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod7 = update(mod2, ~. +Ospedale*Tipo.parto)
summary(mod7)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + Tipo.parto:Ospedale, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1105.96 -182.96 -15.09 160.04 2622.86
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6718.7724 136.9141 -49.073 < 2e-16 ***
## N.gravidanze 12.4288 4.3340 2.868 0.00417 **
## Gestazione 32.1198 3.7921 8.470 < 2e-16 ***
## Lunghezza 10.2964 0.3007 34.246 < 2e-16 ***
## Cranio 10.5110 0.4259 24.680 < 2e-16 ***
## Tipo.partoNat 37.6461 21.0090 1.792 0.07327 .
## Ospedaleosp2 -12.1290 24.6136 -0.493 0.62221
## Ospedaleosp3 48.3611 25.1916 1.920 0.05501 .
## SessoM 77.5161 11.1775 6.935 5.16e-12 ***
## Tipo.partoNat:Ospedaleosp2 1.6347 29.3735 0.056 0.95562
## Tipo.partoNat:Ospedaleosp3 -27.5427 29.8567 -0.922 0.35636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2489 degrees of freedom
## Multiple R-squared: 0.7288, Adjusted R-squared: 0.7277
## F-statistic: 669 on 10 and 2489 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod8 = update(mod1, ~. +Ospedale*Sesso)
summary(mod8)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Ospedale:Sesso,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1120.12 -180.63 -15.29 161.37 2614.15
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6736.3053 141.4894 -47.610 < 2e-16 ***
## Anni.madre 0.8013 1.1468 0.699 0.4848
## N.gravidanze 11.3792 4.6709 2.436 0.0149 *
## Fumatrici -30.5203 27.5648 -1.107 0.2683
## Gestazione 32.5328 3.8202 8.516 < 2e-16 ***
## Lunghezza 10.2932 0.3009 34.208 < 2e-16 ***
## Cranio 10.4740 0.4265 24.557 < 2e-16 ***
## Tipo.partoNat 29.4735 12.0899 2.438 0.0148 *
## Ospedaleosp2 -6.9260 18.8871 -0.367 0.7139
## Ospedaleosp3 27.3555 19.1589 1.428 0.1535
## SessoM 80.0552 19.3724 4.132 3.71e-05 ***
## Ospedaleosp2:SessoM -8.7550 26.8895 -0.326 0.7448
## Ospedaleosp3:SessoM 1.4224 27.0314 0.053 0.9580
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7276
## F-statistic: 557.2 on 12 and 2487 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa.
mod9 = update(mod1, ~. +Gestazione*Lunghezza)
summary(mod9)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + Gestazione:Lunghezza,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1140.11 -179.52 -13.24 161.90 2625.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.129e+03 9.192e+02 -2.316 0.020617 *
## Anni.madre 9.052e-01 1.141e+00 0.793 0.427594
## N.gravidanze 1.179e+01 4.644e+00 2.539 0.011166 *
## Fumatrici -2.710e+01 2.741e+01 -0.989 0.322973
## Gestazione -9.157e+01 2.476e+01 -3.698 0.000222 ***
## Lunghezza 1.430e-01 2.024e+00 0.071 0.943674
## Cranio 1.069e+01 4.261e-01 25.085 < 2e-16 ***
## Tipo.partoNat 2.849e+01 1.203e+01 2.369 0.017907 *
## Ospedaleosp2 -9.527e+00 1.338e+01 -0.712 0.476382
## Ospedaleosp3 2.845e+01 1.343e+01 2.118 0.034258 *
## SessoM 7.193e+01 1.118e+01 6.435 1.48e-10 ***
## Gestazione:Lunghezza 2.682e-01 5.288e-02 5.071 4.25e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.6 on 2488 degrees of freedom
## Multiple R-squared: 0.7317, Adjusted R-squared: 0.7305
## F-statistic: 616.7 on 11 and 2488 DF, p-value: < 2.2e-16car::vif(mod9)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.190614 1 1.091153
## N.gravidanze 1.189563 1 1.090671
## Fumatrici 1.007898 1 1.003941
## Gestazione 72.024973 1 8.486753
## Lunghezza 95.458996 1 9.770312
## Cranio 1.647392 1 1.283508
## Tipo.parto 1.004531 1 1.002263
## Ospedale 1.004939 2 1.001232
## Sesso 1.050972 1 1.025169
## Gestazione:Lunghezza 263.237961 1 16.224610# Interazione statisticamente significativa, ma genera multicollinearità .Proviamo ora a valutare le relazioni non lineari.
mod10 = update(mod1, ~. +I(Gestazione^2))
summary(mod10)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(Gestazione^2),
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1148.24 -179.95 -11.61 163.03 2634.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4830.8240 897.2543 -5.384 7.97e-08 ***
## Anni.madre 0.8782 1.1461 0.766 0.4436
## N.gravidanze 11.3755 4.6631 2.439 0.0148 *
## Fumatrici -28.9924 27.5249 -1.053 0.2923
## Gestazione -73.9035 49.6668 -1.488 0.1369
## Lunghezza 10.3928 0.3039 34.198 < 2e-16 ***
## Cranio 10.5623 0.4278 24.690 < 2e-16 ***
## Tipo.partoNat 29.0581 12.0778 2.406 0.0162 *
## Ospedaleosp2 -10.4856 13.4334 -0.781 0.4351
## Ospedaleosp3 27.7499 13.4884 2.057 0.0398 *
## SessoM 75.4810 11.2111 6.733 2.06e-11 ***
## I(Gestazione^2) 1.4212 0.6612 2.149 0.0317 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.7 on 2488 degrees of freedom
## Multiple R-squared: 0.7294, Adjusted R-squared: 0.7282
## F-statistic: 609.6 on 11 and 2488 DF, p-value: < 2.2e-16car::vif(mod10)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.191462 1 1.091541
## N.gravidanze 1.189267 1 1.090535
## Fumatrici 1.007800 1 1.003892
## Gestazione 287.271649 1 16.949090
## Lunghezza 2.133550 1 1.460668
## Cranio 1.646639 1 1.283215
## Tipo.parto 1.004550 1 1.002272
## Ospedale 1.005731 2 1.001430
## Sesso 1.048359 1 1.023894
## I(Gestazione^2) 280.671253 1 16.753246# Interazione statisticamente abbastanza significativa, ma genera multicollineraritÃ
mod11 = update(mod1, ~. +I(N.gravidanze^2))
summary(mod11) # # Adjusted R-squared 72.81% ##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(N.gravidanze^2),
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1128.37 -179.65 -13.48 162.42 2610.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6734.4377 141.3292 -47.651 < 2e-16 ***
## Anni.madre 0.5107 1.1562 0.442 0.65873
## N.gravidanze 24.6024 8.5097 2.891 0.00387 **
## Fumatrici -32.9140 27.5663 -1.194 0.23259
## Gestazione 33.0335 3.8258 8.634 < 2e-16 ***
## Lunghezza 10.2818 0.3006 34.200 < 2e-16 ***
## Cranio 10.4340 0.4263 24.473 < 2e-16 ***
## Tipo.partoNat 30.3493 12.0875 2.511 0.01211 *
## Ospedaleosp2 -10.5898 13.4365 -0.788 0.43069
## Ospedaleosp3 27.5909 13.4934 2.045 0.04098 *
## SessoM 77.8055 11.1733 6.963 4.23e-12 ***
## I(N.gravidanze^2) -2.4371 1.3151 -1.853 0.06397 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.8 on 2488 degrees of freedom
## Multiple R-squared: 0.7293, Adjusted R-squared: 0.7281
## F-statistic: 609.2 on 11 and 2488 DF, p-value: < 2.2e-16car::vif(mod11) # tutti sotto 5## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.212039 1 1.100926
## N.gravidanze 3.958708 1 1.989650
## Fumatrici 1.010353 1 1.005163
## Gestazione 1.703712 1 1.305263
## Lunghezza 2.086956 1 1.444630
## Cranio 1.634792 1 1.278590
## Tipo.parto 1.005692 1 1.002842
## Ospedale 1.006370 2 1.001589
## Sesso 1.040812 1 1.020202
## I(N.gravidanze^2) 3.571580 1 1.889862# Interazione statisticamente significativa, non genera multicollinearitÃ
mod12 = update(mod11, ~. +I(Lunghezza^2))
summary(mod12)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(N.gravidanze^2) +
## I(Lunghezza^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.31 -180.34 -11.09 158.10 1757.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.82295 724.01608 0.232 0.816717
## Anni.madre 0.33392 1.13527 0.294 0.768682
## N.gravidanze 28.90724 8.36621 3.455 0.000559 ***
## Fumatrici -26.78442 27.07063 -0.989 0.322550
## Gestazione 43.21382 3.89947 11.082 < 2e-16 ***
## Lunghezza -20.24395 3.15649 -6.413 1.70e-10 ***
## Cranio 10.54331 0.41872 25.180 < 2e-16 ***
## Tipo.partoNat 28.08677 11.86923 2.366 0.018041 *
## Ospedaleosp2 -9.02436 13.19231 -0.684 0.494000
## Ospedaleosp3 28.31455 13.24736 2.137 0.032665 *
## SessoM 69.86018 10.99991 6.351 2.54e-10 ***
## I(N.gravidanze^2) -2.91161 1.29200 -2.254 0.024310 *
## I(Lunghezza^2) 0.03166 0.00326 9.713 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.8 on 2487 degrees of freedom
## Multiple R-squared: 0.7391, Adjusted R-squared: 0.7379
## F-statistic: 587.3 on 12 and 2487 DF, p-value: < 2.2e-16car::vif(mod12)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.212350 1 1.101068
## N.gravidanze 3.969848 1 1.992448
## Fumatrici 1.010902 1 1.005436
## Gestazione 1.836371 1 1.355128
## Lunghezza 238.690769 1 15.449620
## Cranio 1.635973 1 1.279052
## Tipo.parto 1.006079 1 1.003035
## Ospedale 1.006521 2 1.001626
## Sesso 1.046599 1 1.023034
## I(N.gravidanze^2) 3.576693 1 1.891215
## I(Lunghezza^2) 230.640954 1 15.186868# Interazione statisticamente significativa, ma genera multicollinearitÃ
mod13 = update(mod11, ~. +I(Cranio^2))
summary(mod13)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(N.gravidanze^2) +
## I(Cranio^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1119.33 -177.95 -13.81 161.61 2593.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.69932 1151.52362 0.081 0.93515
## Anni.madre 0.41263 1.14837 0.359 0.71939
## N.gravidanze 27.16592 8.46192 3.210 0.00134 **
## Fumatrici -29.33366 27.38266 -1.071 0.28416
## Gestazione 39.64599 3.95734 10.018 < 2e-16 ***
## Lunghezza 10.52303 0.30128 34.928 < 2e-16 ***
## Cranio -32.27614 7.16163 -4.507 6.89e-06 ***
## Tipo.partoNat 28.65944 12.00745 2.387 0.01707 *
## Ospedaleosp2 -8.97624 13.34654 -0.673 0.50130
## Ospedaleosp3 29.14488 13.40280 2.175 0.02976 *
## SessoM 72.93803 11.12612 6.556 6.72e-11 ***
## I(N.gravidanze^2) -2.84922 1.30782 -2.179 0.02945 *
## I(Cranio^2) 0.06317 0.01057 5.974 2.64e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 271.9 on 2487 degrees of freedom
## Multiple R-squared: 0.7331, Adjusted R-squared: 0.7318
## F-statistic: 569.2 on 12 and 2487 DF, p-value: < 2.2e-16car::vif(mod13)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.212286 1 1.101039
## N.gravidanze 3.968914 1 1.992213
## Fumatrici 1.010837 1 1.005404
## Gestazione 1.848303 1 1.359523
## Lunghezza 2.125129 1 1.457782
## Cranio 467.700642 1 21.626388
## Tipo.parto 1.006250 1 1.003120
## Ospedale 1.006896 2 1.001720
## Sesso 1.046423 1 1.022948
## I(N.gravidanze^2) 3.581543 1 1.892497
## I(Cranio^2) 455.057600 1 21.332079# Interazione statisticamente significativa, ma genera multicollinearitÃ
mod14 = update(mod11, ~. +Tipo.parto*I(Gestazione^2))
summary(mod14)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(N.gravidanze^2) +
## I(Gestazione^2) + Tipo.parto:I(Gestazione^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1137.50 -180.25 -12.18 161.13 2636.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.423e+03 9.195e+02 -4.810 1.60e-06 ***
## Anni.madre 5.054e-01 1.156e+00 0.437 0.66205
## N.gravidanze 2.568e+01 8.513e+00 3.016 0.00258 **
## Fumatrici -3.243e+01 2.754e+01 -1.177 0.23916
## Gestazione -8.713e+01 4.995e+01 -1.744 0.08122 .
## Lunghezza 1.038e+01 3.037e-01 34.179 < 2e-16 ***
## Cranio 1.052e+01 4.278e-01 24.597 < 2e-16 ***
## Tipo.partoNat -1.740e+02 1.391e+02 -1.252 0.21084
## Ospedaleosp2 -9.817e+00 1.343e+01 -0.731 0.46475
## Ospedaleosp3 2.659e+01 1.349e+01 1.972 0.04873 *
## SessoM 7.619e+01 1.121e+01 6.797 1.33e-11 ***
## I(N.gravidanze^2) -2.645e+00 1.317e+00 -2.009 0.04467 *
## I(Gestazione^2) 1.507e+00 6.622e-01 2.275 0.02297 *
## Tipo.partoNat:I(Gestazione^2) 1.338e-01 9.084e-02 1.472 0.14102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.5 on 2486 degrees of freedom
## Multiple R-squared: 0.7301, Adjusted R-squared: 0.7286
## F-statistic: 517.2 on 13 and 2486 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa
mod15 = update(mod11, ~. +I(Anni.madre^2))
summary(mod15)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + I(N.gravidanze^2) +
## I(Anni.madre^2), data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1094.5 -180.2 -13.5 160.5 2612.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6894.8293 177.1640 -38.918 < 2e-16 ***
## Anni.madre 13.1766 8.5189 1.547 0.12205
## N.gravidanze 24.6311 8.5076 2.895 0.00382 **
## Fumatrici -32.4640 27.5610 -1.178 0.23895
## Gestazione 32.7698 3.8289 8.559 < 2e-16 ***
## Lunghezza 10.2718 0.3006 34.167 < 2e-16 ***
## Cranio 10.4379 0.4263 24.488 < 2e-16 ***
## Tipo.partoNat 29.8502 12.0890 2.469 0.01361 *
## Ospedaleosp2 -10.5726 13.4331 -0.787 0.43133
## Ospedaleosp3 27.2508 13.4919 2.020 0.04351 *
## SessoM 78.3825 11.1771 7.013 3e-12 ***
## I(N.gravidanze^2) -2.3473 1.3161 -1.784 0.07462 .
## I(Anni.madre^2) -0.2225 0.1483 -1.501 0.13357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.7 on 2487 degrees of freedom
## Multiple R-squared: 0.7295, Adjusted R-squared: 0.7282
## F-statistic: 558.9 on 12 and 2487 DF, p-value: < 2.2e-16# Interazione non statisticamente significativa2c. Selezione del modello
Non resta che analizzare la bontà dei modelli ottenuti tramite i criteri di informazione Akaike e Bayes. I due criteri portano entrambi alla scelta di mod2.
AIC(mod1,mod2,mod11)## df AIC
## mod1 12 35172.09
## mod2 10 35169.79
## mod11 13 35170.64# Vince mod2
BIC(mod1,mod2,mod11)## df BIC
## mod1 12 35241.97
## mod2 10 35228.03
## mod11 13 35246.35# Vince mod2
car::vif(mod1) # assenti multicollinearità significative## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.190207 1 1.090966
## N.gravidanze 1.189252 1 1.090528
## Fumatrici 1.007410 1 1.003698
## Gestazione 1.695001 1 1.301922
## Lunghezza 2.085773 1 1.444220
## Cranio 1.630914 1 1.277073
## Tipo.parto 1.004255 1 1.002125
## Ospedale 1.004235 2 1.001057
## Sesso 1.040650 1 1.020123car::vif(mod2) # assenti multicollinearità significative## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.025244 1 1.012544
## Gestazione 1.670237 1 1.292376
## Lunghezza 2.081921 1 1.442886
## Cranio 1.626507 1 1.275346
## Tipo.parto 1.003872 1 1.001934
## Ospedale 1.003038 2 1.000759
## Sesso 1.040337 1 1.019969car::vif(mod11) # assenti multicollinearità significative, numero gravidanze con VIF=3.6## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.212039 1 1.100926
## N.gravidanze 3.958708 1 1.989650
## Fumatrici 1.010353 1 1.005163
## Gestazione 1.703712 1 1.305263
## Lunghezza 2.086956 1 1.444630
## Cranio 1.634792 1 1.278590
## Tipo.parto 1.005692 1 1.002842
## Ospedale 1.006370 2 1.001589
## Sesso 1.040812 1 1.020202
## I(N.gravidanze^2) 3.571580 1 1.8898622d. Analisi qualità del modello
Per la valutazione del parametro Adjusted R-squared possiamo fare uso della funzione summary che restituisce un valore del 72.78%, non molto alto.
Per la valutazione del parametro MSE meglio riferirsi alla sua radice quadrata, RMSE, da confrontare con la varianza della variabile risposta. RMSE (273)<<Peso_sd (525), segno che il modello ha una buona capacità predittiva.
# Adj R-squared
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati)
##
## 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# RMSE
Peso_est = predict(mod2, neonati)
RMSE = sqrt(mean((Peso - Peso_est)^2))
RMSE## [1] 273.4227Peso_sd## [1] 525.0387Passiamo ora ad un’analisi dei residui, partendo da una valutazione grafica:
I residui si dispongono generalmente attorno alla media nulla, ma non sulle code.
La maggior parte dei residui si allinea al QQplot della normale, ma non sulle code.
L’andamento della varianza dei residui è abbastanza costante, ma non sulle code.
Il test sulla distribuzione normale dei residui fallisce, ma considerato che la stessa variabile peso non è per sua natura normale lo riterrei accettabile. Fallisce il test sulla omoschedasticità , come si vede dal relativo grafico. Fortunatamente passa il test sulla non correlazione.
par(mfrow=c(2,2))
plot(mod2)
# Test normalità residui
shapiro.test(residuals(mod2))##
## Shapiro-Wilk normality test
##
## data: residuals(mod2)
## W = 0.97408, p-value < 2.2e-16# pvalue nullo, si rifiuta ipotesi di normalitÃ
# Test omoschedasticità residui
lmtest::bptest(mod2)##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 91.768, df = 8, p-value < 2.2e-16# pvalue nullo, si rifiuta ipotesi di omoschedasticitÃ
# Test autocorrelazione
lmtest::dwtest(mod2)##
## Durbin-Watson test
##
## data: mod2
## DW = 1.9527, p-value = 0.1184
## alternative hypothesis: true autocorrelation is greater than 0# pvalue alto, non si rifiuta l'ipotesi di non correlazionePassiamo ora alla valutazione di leverage e outlier:
Le osservazioni presentano tutte distanza di Cook bassa, ad eccezione di O1551 per cui è piuttosto alta (0.56). In totale ci sono 96 osservazioni ad alto leverage nello spazio dei regressori.
Sono presenti 4 outlier significativi della variabile risposta (O1551, O155, O1306).
# La distanza di Cook misura l'influenza combinata di leverage (distanza dei valori indipendenti) e del residuo (distanza della predizione dal valore reale). Un'osservazione con alta distanza di Cook è considerata influente perché cambia i coefficienti della regressione se viene rimossa.
cook = cooks.distance(mod2)
plot(cook)
max(cook) ## [1] 0.5557527which.max(cook) # O1551 con distanza critica (0.56)## 1551
## 1551which(cook > 0.5)## 1551
## 1551neonati[which(cook > 0.5), ] # Identificazione obs con distanza Cook > 0.5## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
## 1551 35 1 0 38 4370 315 374 Nat osp3 F# Leverage misura quanto un punto è lontano dalla media delle variabili indipendenti. Un'osservazione con alto leverage ha il potenziale per essere influente, ma non necessariamente lo è perché potrebbe avere un piccolo residuo. Ci limiteremo quindi all'analisi della distanza di Cook, il codice di seguito è riportato solo per completezza.
par(mfrow=c(1,1))
lev = hatvalues(mod2)
plot(lev)
p = sum(lev)
soglia = 2*p/nrow(neonati)
abline(h=soglia,col=2) 
sum(lev > soglia) # 96 osservazioni ad alto leverage nello spazio dei regressori## [1] 96lev[lev > soglia]## 13 15 34 89 101 106 131 134 151 155
## 0.007361639 0.008775066 0.007643421 0.014420225 0.008426883 0.016173741 0.008130949 0.008471856 0.012661857 0.008046220
## 161 204 206 220 294 310 312 378 442 445
## 0.021627904 0.015436038 0.010436735 0.008270600 0.007586623 0.029729279 0.014137567 0.016839143 0.008861521 0.008412092
## 492 516 582 587 592 638 684 697 748 750
## 0.009242952 0.013908498 0.012606967 0.009287608 0.008383030 0.007581023 0.010870924 0.007611658 0.009450875 0.007855118
## 757 805 828 913 928 946 947 956 985 1014
## 0.009114441 0.015287466 0.008231836 0.007343981 0.023513670 0.007869751 0.009313671 0.008619578 0.008572090 0.009409396
## 1067 1091 1130 1188 1219 1248 1273 1291 1311 1321
## 0.009421177 0.009900661 0.032772693 0.008116936 0.032134700 0.015556185 0.008781482 0.007228173 0.010581446 0.010293353
## 1357 1385 1411 1428 1429 1450 1505 1551 1553 1556
## 0.007914331 0.013700942 0.009986260 0.009158583 0.022532944 0.015916785 0.015080350 0.049396095 0.009500259 0.007664618
## 1610 1619 1686 1701 1712 1718 1727 1780 1781 1809
## 0.010346542 0.017207090 0.011378179 0.011661008 0.008125639 0.008051861 0.014494385 0.026525402 0.018109714 0.010507291
## 1977 2040 2086 2089 2114 2115 2120 2140 2148 2149
## 0.008671741 0.013591813 0.014459304 0.008203459 0.014249699 0.012646408 0.019714874 0.008049584 0.009815107 0.014447413
## 2175 2200 2216 2221 2224 2244 2307 2317 2359 2408
## 0.033456296 0.012545760 0.009073495 0.022730816 0.007812558 0.007810288 0.014903596 0.008640738 0.010984743 0.010657742
## 2422 2437 2452 2458 2471 2478
## 0.022330175 0.024704776 0.026001843 0.009562118 0.022832989 0.007416139neonati[lev > soglia, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
## 13 36 5 0 38 3060 455 325 Ces osp1 F
## 15 33 3 0 34 2400 470 298 Ces osp3 M
## 34 27 0 0 39 3150 480 382 Nat osp1 F
## 89 36 8 0 39 3610 500 351 Ces osp1 M
## 101 31 0 0 34 1370 390 287 Nat osp2 F
## 106 29 4 0 30 1340 400 273 Ces osp1 M
## 131 30 0 0 34 2290 450 285 Nat osp2 M
## 134 38 6 0 37 3950 500 350 Nat osp3 M
## 151 20 0 0 41 2280 450 280 Ces osp3 M
## 155 30 0 0 36 3610 410 330 Nat osp1 M
## 161 35 9 0 42 3760 540 348 Nat osp2 F
## 204 30 8 0 40 3850 518 340 Nat osp3 F
## 206 39 1 0 31 1500 405 295 Nat osp3 M
## 220 23 1 0 40 3520 445 363 Nat osp1 F
## 294 34 5 0 40 3100 500 330 Ces osp2 F
## 310 40 3 0 28 1560 420 379 Nat osp3 F
## 312 26 1 0 32 1280 360 276 Nat osp2 M
## 378 27 0 0 28 1285 400 274 Nat osp1 F
## 442 35 6 1 38 2430 460 324 Nat osp2 F
## 445 27 0 0 32 1550 410 289 Nat osp1 F
## 492 34 2 0 33 1410 380 295 Nat osp2 F
## 516 40 8 0 38 3520 470 341 Nat osp3 M
## 582 30 7 0 35 2220 470 316 Nat osp3 M
## 587 16 1 0 31 1900 440 300 Nat osp2 F
## 592 30 1 0 32 2260 440 322 Ces osp3 F
## 638 25 0 0 33 1720 420 300 Nat osp1 M
## 684 30 1 0 39 3000 475 390 Ces osp2 F
## 697 30 0 0 39 2820 510 300 Ces osp3 F
## 748 35 0 0 33 1390 390 277 Nat osp1 F
## 750 24 0 0 35 1450 405 280 Nat osp1 F
## 757 30 6 0 35 2680 450 322 Nat osp1 F
## 805 30 2 0 29 1190 360 272 Nat osp2 F
## 828 32 6 0 40 4200 510 350 Nat osp1 M
## 913 34 5 0 40 3060 485 335 Ces osp2 F
## 928 25 0 0 28 830 310 254 Nat osp1 F
## 946 36 5 0 35 2900 490 340 Nat osp3 M
## 947 34 3 0 32 1615 390 297 Nat osp3 F
## 956 25 0 0 41 2210 430 310 Nat osp3 F
## 985 24 5 0 42 3600 510 335 Ces osp3 M
## 1014 17 0 0 37 2050 390 295 Nat osp2 F
## 1067 26 3 0 31 1960 420 300 Nat osp2 F
## 1091 30 1 0 33 1770 410 275 Nat osp3 M
## 1130 33 11 0 43 3400 475 360 Nat osp1 M
## 1188 21 0 0 40 4140 550 320 Ces osp1 M
## 1219 38 12 0 39 3350 490 344 Nat osp2 M
## 1248 26 1 0 30 1170 370 266 Nat osp2 M
## 1273 32 1 0 33 2040 480 307 Ces osp1 F
## 1291 39 5 0 37 3360 490 320 Nat osp2 M
## 1311 40 6 0 34 1840 430 305 Nat osp1 F
## 1321 36 6 0 40 4340 515 383 Nat osp1 M
## 1357 22 0 0 32 2340 445 304 Nat osp1 F
## 1385 33 0 0 29 1620 410 292 Nat osp3 F
## 1411 32 6 0 39 3290 480 360 Ces osp2 M
## 1428 30 1 0 36 1280 385 292 Nat osp2 F
## 1429 24 4 0 29 1280 390 355 Nat osp1 F
## 1450 36 8 0 41 3730 480 335 Nat osp3 M
## 1505 30 8 0 39 2860 490 337 Ces osp2 F
## 1551 35 1 0 38 4370 315 374 Nat osp3 F
## 1553 30 4 0 35 4520 520 360 Nat osp2 F
## 1556 37 0 0 41 2420 490 300 Ces osp1 M
## 1610 37 3 0 33 2000 470 293 Ces osp1 F
## 1619 31 0 0 31 990 340 278 Ces osp2 F
## 1686 27 0 0 31 1800 430 308 Ces osp3 M
## 1701 22 0 0 32 1430 380 301 Nat osp1 M
## 1712 28 0 0 39 3800 520 300 Nat osp3 F
## 1718 34 4 0 42 2660 500 320 Nat osp2 F
## 1727 36 8 0 36 2860 460 334 Nat osp2 F
## 1780 25 2 0 25 900 325 253 Nat osp3 F
## 1781 35 9 0 37 3150 490 335 Nat osp2 M
## 1809 35 0 0 32 1780 420 277 Ces osp1 F
## 1977 39 4 0 34 2970 480 350 Ces osp2 F
## 2040 27 1 0 38 3240 410 359 Ces osp1 F
## 2086 26 8 0 40 3250 500 355 Nat osp2 M
## 2089 32 1 1 33 1780 400 305 Ces osp1 F
## 2114 36 0 0 31 1180 355 270 Nat osp3 F
## 2115 35 1 0 32 1890 500 309 Nat osp2 F
## 2120 32 0 0 27 1140 370 267 Nat osp3 F
## 2140 30 2 0 33 1600 410 290 Ces osp1 F
## 2148 33 6 0 37 2900 465 350 Ces osp2 F
## 2149 39 3 0 30 1300 380 276 Nat osp1 M
## 2175 37 8 0 28 930 355 235 Nat osp1 F
## 2200 33 0 0 30 1750 410 294 Nat osp2 M
## 2216 22 0 0 32 2580 470 330 Nat osp1 F
## 2221 35 10 0 39 2950 495 335 Nat osp1 F
## 2224 41 1 0 33 2000 425 312 Ces osp3 M
## 2244 35 6 0 39 3300 500 350 Nat osp3 M
## 2307 26 1 0 30 1170 370 273 Nat osp3 M
## 2317 25 6 0 38 2530 460 340 Nat osp1 F
## 2359 25 6 0 33 2230 430 313 Nat osp3 F
## 2408 37 2 0 31 1690 405 290 Nat osp2 M
## 2422 33 10 0 40 3090 485 353 Nat osp3 M
## 2437 28 1 0 27 980 320 265 Nat osp1 M
## 2452 28 0 0 26 930 345 245 Ces osp3 F
## 2458 31 0 0 31 1730 430 300 Nat osp3 F
## 2471 34 10 0 38 2880 470 345 Ces osp2 M
## 2478 32 1 0 33 2740 475 324 Ces osp2 F# Outlier
plot(rstudent(mod2))
abline(h=c(-2,2), col=2)
res_stud = rstudent(mod2)
obs_outlier = which(res_stud < -2 | res_stud > 2)
num_obs_outlier = length(obs_outlier) # 104 osservazioni outlier nello spazio della risposta
num_obs_outlier## [1] 104neonati[obs_outlier, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
## 5 20 0 0 38 3700 480 335 Nat osp3 F
## 90 29 1 0 39 3800 490 333 Ces osp3 M
## 119 31 0 0 40 3410 550 372 Nat osp2 M
## 130 30 2 0 39 4240 485 352 Nat osp2 M
## 140 29 1 1 41 4420 530 362 Ces osp2 F
## 146 24 1 0 40 3820 500 320 Ces osp2 F
## 155 30 0 0 36 3610 410 330 Nat osp1 M
## 262 28 1 0 39 2600 480 350 Nat osp1 F
## 295 18 0 0 40 1850 460 305 Nat osp3 F
## 310 40 3 0 28 1560 420 379 Nat osp3 F
## 318 21 0 0 39 3580 485 320 Nat osp3 M
## 329 25 1 0 40 4560 540 340 Nat osp1 M
## 361 25 1 0 39 2720 495 345 Nat osp3 F
## 364 22 2 0 38 2410 470 336 Nat osp2 M
## 375 33 4 0 38 4270 510 353 Nat osp2 M
## 377 31 1 0 39 3400 460 325 Nat osp1 M
## 390 38 0 0 40 3700 470 320 Nat osp1 M
## 403 31 0 0 41 4110 500 350 Nat osp2 M
## 418 26 0 0 40 3580 480 326 Nat osp1 F
## 455 24 0 0 38 2920 505 351 Ces osp3 M
## 472 25 2 0 41 3990 495 335 Nat osp3 M
## 582 30 7 0 35 2220 470 316 Nat osp3 M
## 616 29 0 0 42 3540 540 368 Nat osp1 M
## 623 34 0 0 41 3940 500 340 Nat osp1 F
## 632 21 0 0 37 2750 510 333 Nat osp1 M
## 633 30 0 0 40 3860 480 338 Nat osp1 M
## 648 35 1 0 38 2850 500 360 Nat osp1 F
## 653 21 0 0 41 2700 500 336 Nat osp2 F
## 709 29 1 0 38 4130 520 349 Nat osp1 F
## 762 28 0 0 39 4120 520 343 Nat osp1 F
## 791 30 1 0 41 4440 510 335 Nat osp3 M
## 850 28 1 0 39 3920 485 345 Nat osp2 F
## 890 33 2 0 39 3660 490 318 Nat osp2 F
## 901 24 0 0 39 2930 500 350 Nat osp3 M
## 908 35 2 0 39 3300 460 320 Ces osp1 M
## 928 25 0 0 28 830 310 254 Nat osp1 F
## 950 25 0 0 41 2980 500 356 Nat osp3 M
## 1036 26 0 0 40 4330 500 355 Nat osp3 F
## 1132 21 0 0 39 3950 500 330 Nat osp1 F
## 1137 26 0 0 40 3000 520 340 Ces osp3 M
## 1192 24 2 0 40 2560 490 323 Nat osp1 M
## 1194 41 1 0 39 2160 450 325 Ces osp3 M
## 1207 23 1 0 39 4100 505 350 Nat osp1 F
## 1211 34 1 0 40 2580 495 325 Nat osp3 F
## 1230 29 1 0 41 4010 470 353 Nat osp2 M
## 1268 28 1 0 40 3790 460 332 Nat osp2 F
## 1287 35 3 0 40 2660 485 335 Ces osp3 M
## 1293 30 3 0 38 4600 485 380 Nat osp1 M
## 1297 29 2 0 38 3870 500 343 Nat osp2 F
## 1306 23 0 0 41 4900 510 352 Nat osp2 F
## 1341 32 3 0 36 3780 480 349 Nat osp3 M
## 1360 31 1 0 41 2940 515 337 Nat osp1 F
## 1395 30 2 0 39 3790 505 304 Ces osp3 M
## 1399 42 2 0 38 2560 525 349 Ces osp2 M
## 1429 24 4 0 29 1280 390 355 Nat osp1 F
## 1433 31 4 0 41 4810 530 364 Ces osp3 M
## 1499 30 1 0 39 2380 495 325 Nat osp1 F
## 1541 30 0 0 38 4540 530 343 Ces osp3 M
## 1551 35 1 0 38 4370 315 374 Nat osp3 F
## 1553 30 4 0 35 4520 520 360 Nat osp2 F
## 1585 21 0 0 40 4050 500 340 Nat osp3 M
## 1588 22 0 0 42 3100 510 361 Nat osp2 F
## 1593 41 3 0 35 1500 420 304 Nat osp1 M
## 1635 32 2 0 39 3430 445 322 Ces osp1 F
## 1639 39 3 0 40 4760 550 365 Nat osp2 F
## 1694 23 1 0 36 3850 460 334 Ces osp3 F
## 1712 28 0 0 39 3800 520 300 Nat osp3 F
## 1718 34 4 0 42 2660 500 320 Nat osp2 F
## 1780 25 2 0 25 900 325 253 Nat osp3 F
## 1837 29 1 0 41 3940 500 340 Nat osp1 F
## 1838 34 3 0 40 4580 515 360 Nat osp2 M
## 1856 33 0 0 36 3110 465 300 Ces osp3 M
## 1868 29 0 0 40 3470 525 390 Nat osp1 M
## 1915 27 1 0 38 2480 500 315 Ces osp2 M
## 1920 26 0 1 39 4930 550 350 Ces osp2 F
## 1937 36 1 0 41 2940 500 358 Nat osp3 M
## 1944 31 1 0 38 3920 490 345 Nat osp2 M
## 1962 38 4 0 38 4370 530 340 Nat osp3 F
## 1963 27 0 0 42 4700 540 362 Nat osp3 F
## 1993 31 1 0 40 4080 510 343 Nat osp2 F
## 2012 24 0 0 38 3920 505 333 Nat osp3 M
## 2023 27 1 0 39 4650 510 354 Nat osp2 M
## 2040 27 1 0 38 3240 410 359 Ces osp1 F
## 2076 23 0 0 40 4720 540 360 Nat osp3 F
## 2115 35 1 0 32 1890 500 309 Nat osp2 F
## 2123 32 0 0 39 3950 520 330 Ces osp1 F
## 2135 33 3 0 39 3950 480 345 Nat osp3 M
## 2151 34 0 0 39 2690 495 335 Nat osp1 M
## 2179 25 2 0 41 4010 500 345 Nat osp3 F
## 2185 25 0 0 40 2990 515 350 Ces osp1 F
## 2195 38 1 0 40 3980 480 335 Nat osp1 F
## 2204 33 0 0 39 4000 490 360 Ces osp3 F
## 2219 37 1 0 39 2500 490 352 Ces osp1 M
## 2225 27 0 0 35 3140 465 290 Nat osp2 F
## 2287 33 3 0 40 3640 470 337 Nat osp3 F
## 2315 24 0 0 42 2800 520 340 Ces osp2 M
## 2343 23 1 0 39 4250 520 350 Nat osp1 M
## 2370 27 1 0 40 3890 495 330 Nat osp2 M
## 2392 28 1 0 40 4720 540 355 Nat osp3 M
## 2398 26 0 0 40 2660 500 343 Nat osp1 F
## [ reached 'max' / getOption("max.print") -- omitted 4 rows ]car::outlierTest(mod2) # test di significatività outlier## rstudent unadjusted p-value Bonferroni p
## 1551 10.004278 3.9642e-23 9.9104e-20
## 155 5.046640 4.8210e-07 1.2053e-03
## 1306 4.872419 1.1716e-06 2.9291e-03Proviamo quindi a rimuovere le osservazioni critiche.
obs_removal = c(155, 1306, 1551)
neonati_cleaned <- neonati[-obs_removal, ]
mod2_cleaned = update(mod2, data = neonati_cleaned)
summary(mod2_cleaned)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati_cleaned)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1132.9 -180.1 -13.4 160.4 1132.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6720.0967 132.0938 -50.874 < 2e-16 ***
## N.gravidanze 13.8472 4.2076 3.291 0.00101 **
## Gestazione 28.9007 3.6889 7.835 6.91e-15 ***
## Lunghezza 11.0647 0.2991 36.992 < 2e-16 ***
## Cranio 9.7786 0.4177 23.408 < 2e-16 ***
## Tipo.partoNat 27.8811 11.7292 2.377 0.01753 *
## Ospedaleosp2 -11.6119 13.0500 -0.890 0.37366
## Ospedaleosp3 27.1593 13.1012 2.073 0.03827 *
## SessoM 76.9682 10.8565 7.090 1.75e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.9 on 2488 degrees of freedom
## Multiple R-squared: 0.7432, Adjusted R-squared: 0.7424
## F-statistic: 900.1 on 8 and 2488 DF, p-value: < 2.2e-16Ripetiamo a questo punto il test sui residui per il nuovo modello, mod2_cleaned:
Graficamente non si notano grandi variazioni rispetto allo scenario precedente.
Il grosso miglioramento è che ora i residui passano il test di omoschedasticità .
par(mfrow=c(2,2))
plot(mod2_cleaned)
# Test normalità residui
shapiro.test(residuals(mod2_cleaned))##
## Shapiro-Wilk normality test
##
## data: residuals(mod2_cleaned)
## W = 0.99221, p-value = 2.522e-10# pvalue nullo, si rifiuta ipotesi di normalitÃ
# Test omoschedasticità residui
lmtest::bptest(mod2_cleaned)##
## studentized Breusch-Pagan test
##
## data: mod2_cleaned
## BP = 12.413, df = 8, p-value = 0.1337# pvalue alto, non si rifiuta ipotesi di omoschedasticitÃ
# Test autocorrelazione
lmtest::dwtest(mod2_cleaned)##
## Durbin-Watson test
##
## data: mod2_cleaned
## DW = 1.9539, p-value = 0.1247
## alternative hypothesis: true autocorrelation is greater than 0# pvalue alto, non si rifiuta l'ipotesi di non correlazioneRivalutiamo a questo punto la bontà del modello, mod2_cleaned, selezionato:
Adjusted R-squared restituisce un valore del 74.24%, non altissimo ma il migliore sin qui trovato.
RMSE (265) < Peso_sd (524), segno che il modello mantiene una buona capacità predittiva.
# Adj R-squared
summary(mod2_cleaned)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati_cleaned)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1132.9 -180.1 -13.4 160.4 1132.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6720.0967 132.0938 -50.874 < 2e-16 ***
## N.gravidanze 13.8472 4.2076 3.291 0.00101 **
## Gestazione 28.9007 3.6889 7.835 6.91e-15 ***
## Lunghezza 11.0647 0.2991 36.992 < 2e-16 ***
## Cranio 9.7786 0.4177 23.408 < 2e-16 ***
## Tipo.partoNat 27.8811 11.7292 2.377 0.01753 *
## Ospedaleosp2 -11.6119 13.0500 -0.890 0.37366
## Ospedaleosp3 27.1593 13.1012 2.073 0.03827 *
## SessoM 76.9682 10.8565 7.090 1.75e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.9 on 2488 degrees of freedom
## Multiple R-squared: 0.7432, Adjusted R-squared: 0.7424
## F-statistic: 900.1 on 8 and 2488 DF, p-value: < 2.2e-16# RMSE
Peso_est = predict(mod2_cleaned, neonati_cleaned)
RMSE = sqrt(mean((neonati_cleaned$Peso - Peso_est)^2))
RMSE## [1] 265.4121sd(neonati_cleaned$Peso)## [1] 523.86493. Previsioni e risultati
Vogliamo ora stimare il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice, che partorirà alla 39^ settimana. Per le variabili regressori non specificate, faremo uso di valori mediani o modali a seconda del tipo di variabile.
Il modello ci fornisce un valore di peso atteso pari a 3308g.
neonata = data.frame(N.gravidanze = 2, Fumatrici = 0, Gestazione = 39,
Anni.madre = median(neonati_cleaned$Anni.madre),
Lunghezza = median(neonati_cleaned$Lunghezza),
Cranio = median(neonati_cleaned$Cranio),
Tipo.parto = "Nat", Ospedale = "osp2", Sesso = "F")
Peso_est_neonata = predict(mod2_cleaned, neonata)
Peso_est_neonata## 1
## 3308.0664. Visualizzazioni
Mostriamo infine alcune visualizzazioni relative ai risultati del modello.
i. Gestazione vs Peso, scatterplot
Visualizza le osservazioni del peso vs settimane di gestazione e sovrappone la linea di regressione del modello, per mostrare come si adatta ai dati.
neonati_cleaned$Peso_est = predict(mod2_cleaned, neonati_cleaned)
ggplot(neonati_cleaned, aes(x = Gestazione)) +
geom_point(aes(y = Peso), alpha = 0.6, color = "blue") +
geom_smooth(aes(y = Peso_est), method = "lm", col = "red", se = FALSE) +
labs(title = "Gestazione vs Peso",
x = "Settimane di Gestazione", y = "Peso (g)") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
ii. Confronto tra valori osservati e predetti, scatterplot
Mostra i valori predetti dal modello vs quelli osservati. Se il modello è buono, i punti dovrebbero allinearsi lungo la linea diagonale (y = x).
ggplot(neonati_cleaned, aes(x = Peso, y = Peso_est)) +
geom_point(color = "blue", alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
labs(title = "Valori Predetti vs Osservati",
x = "Peso Osservato", y = "Peso Predetto") +
theme_minimal()
iii. Distribuzione dei residui, istogramma
Un modello buono dovrebbe avere dei residui distribuiti simmetricamente attorno allo zero.
neonati_cleaned$res = neonati_cleaned$Peso - neonati_cleaned$Peso_est
ggplot(neonati_cleaned, aes(x = res)) +
geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
labs(title = "Distribuzione dei residui",
x = "Residuo", y = "Frequenza") +
theme_minimal()
i. Sesso vs Peso, boxplot
Confrontano i pesi predetti dal modello in base a variabili categoriche, come il sesso.
ggplot(neonati_cleaned, aes(x = Sesso, y = Peso_est, fill = Sesso)) +
geom_boxplot() +
stat_summary(fun = median, geom = "text", aes(label = round(..y.., 0)),
vjust = -0.5, color = "black", size = 3) +
stat_summary(fun = function(x) quantile(x, 0.25),
geom = "text", aes(label = round(..y.., 0)),
vjust = 1.5, color = "blue", size = 3) +
stat_summary(fun = function(x) quantile(x, 0.75),
geom = "text", aes(label = round(..y.., 0)),
vjust = -1.5, color = "blue", size = 3) +
labs(title = "Previsioni Peso per Sesso", x = "Sesso", y = "Peso Predetto (g)") +
theme_minimal()