Angelo Caramia - Modello Statistico per la previsione peso di neonati
Angelo Caramia
14/12/2024
R Markdown
dati<-read.csv("C:/Users/angel/Desktop/Data Science/R statistica inferenziale/neonati.csv")
summary(dati)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
## #Analisi Preliminare
Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie. Inoltre si saggeranno le seguenti ipotesi con i test adatti: in alcuni ospedali si fanno più parti cesarei La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione Le misure antropometriche sono significativamente diverse tra i due sessi
Calcolo di statistiche descrittive
library(psych)## Warning: il pacchetto 'psych' è stato creato con R versione 4.1.3describe(dati) ## vars n mean sd median trimmed mad min max range
## Anni.madre 1 2500 28.16 5.27 28 28.10 4.45 0 46 46
## N.gravidanze 2 2500 0.98 1.28 1 0.74 1.48 0 12 12
## Fumatrici 3 2500 0.04 0.20 0 0.00 0.00 0 1 1
## Gestazione 4 2500 38.98 1.87 39 39.19 1.48 25 43 18
## Peso 5 2500 3284.08 525.04 3300 3302.90 459.61 830 4930 4100
## Lunghezza 6 2500 494.69 26.32 500 496.45 22.24 310 565 255
## Cranio 7 2500 340.03 16.43 340 340.68 14.83 235 390 155
## Tipo.parto* 8 2500 1.71 0.45 2 1.76 0.00 1 2 1
## Ospedale* 9 2500 2.01 0.81 2 2.01 1.48 1 3 2
## Sesso* 10 2500 1.50 0.50 1 1.50 0.00 1 2 1
## skew kurtosis se
## Anni.madre 0.04 0.38 0.11
## N.gravidanze 2.51 10.98 0.03
## Fumatrici 4.59 19.06 0.00
## Gestazione -2.06 8.25 0.04
## Peso -0.65 2.03 10.50
## Lunghezza -1.51 6.48 0.53
## Cranio -0.78 2.94 0.33
## Tipo.parto* -0.92 -1.16 0.01
## Ospedale* -0.01 -1.49 0.02
## Sesso* 0.01 -2.00 0.01Statistiche per tutte le variabili è mostrata la media, la deviazione standard, la mediana, minimo e massimo, il numero di elementi, range, curtosi e asimmetria, errore standard e le altre descrittive per le variabili in esame. notiamo il massimo della gestazione di 43 settimane è un valore limite secondo la letteratura ma ancora accettabile sicuramente non vale per il valore minimo di età 0 della madre
Individua eventuiali valori mancanti
colSums(is.na(dati)) ## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0Analisi grafica delle distribuzioni Istogrammi e grafici di tutte le variabili ad eccezione di Tipo parto, fumatrici e Ospedale che verranno visualizzate in tabella a doppia entrata
par(mfrow = c(3, 2))
hist(dati$Peso, main = "Distribuzione Peso Neonatale", xlab = "Peso (grammi)", col = "lightblue")
boxplot(dati$Peso, main = "Boxplot Peso Neonatale", col = "lightblue", horizontal = TRUE)
hist(dati$Gestazione, main = "Distribuzione Settimane di Gestazione", xlab = "Settimane", col = "lightgreen")
boxplot(dati$Gestazione, main = "Boxplot Settimane di Gestazione", col = "lightgreen", horizontal = TRUE)
hist(dati$Lunghezza, main = "Distribuzione Lunghezza Neonati", xlab = "Lunghezza (mm)", col = "lightpink")
boxplot(dati$Lunghezza, main = "Boxplot Lunghezza Neonati", col = "lightpink", horizontal = TRUE)
Gli istogrammi presentano una forma simil gaussiana. Andremo a verificare che la media non sia diversa da quella della popolazione. I box plot espricitano la mediana e come si distribuiscono i valori andando ad evidenziare outlier. Per settimane e lunghezza notiamo molti outlier nei valori inferiori e quasi nessuno per valori superiori. Per il peso invece ci sono anche outlier verso l’estremo superiore anche se si notano anche qui molti più elementi con valori più bassi.
par(mfrow = c(2, 2))
hist(dati$Anni.madre, main = "Distribuzione Anni madre", xlab = "Anni", col = "lightblue")
boxplot(dati$Anni.madre, main = "Boxplot Anni madre", col = "lightblue", horizontal = TRUE)
hist(dati$Cranio, main = "Distribuzione Cranio Neonati", xlab = "Circonferenza (mm)", col = "lightpink")
boxplot(dati$Cranio, main = "Boxplot Cranio Neonati", col = "lightpink", horizontal = TRUE)
Anche qui si nota una forma simil gaussiana oltre a degli outlier anomali per l’età della madre essendo molto bassi 2 valori
Tabella che divide fumatrici (1) da non fumatrici (0)
table(dati$Fumatrici)##
## 0 1
## 2396 104Identificazione di outlier per Peso, gestazione, Lunghezza, Cranio e Anni Madre
par(mfrow = c(1, 1))
peso_outlier <- boxplot.stats(dati$Peso)$out
peso_outlier # valori anomali del peso## [1] 1370 1340 4680 1500 1850 1560 1280 1750 4600 1285 1550 1410 1900 1720 1980
## [16] 1390 1450 1970 1190 2000 830 1615 1960 1770 1750 1170 2040 1980 4600 4900
## [31] 1840 1620 1280 1280 4810 2040 4620 1500 2000 990 4760 1800 1430 1950 1970
## [46] 900 1780 4580 4930 4700 4650 4720 1780 1180 1890 1140 1600 1300 930 1750
## [61] 2000 4690 1580 1170 4720 1690 980 930 1730gestazione_outlier <- boxplot.stats(dati$Gestazione)$out
gestazione_outlier ## [1] 34 33 34 30 34 34 31 34 33 28 32 28 34 34 32 33 33 31 32 34 33 33 33 33 29
## [26] 34 28 32 31 33 34 34 33 30 33 34 32 29 34 29 33 31 31 32 25 32 33 34 34 33
## [51] 31 32 27 33 30 28 30 32 33 34 30 33 31 27 26 31 33Lunghezza_outlier <- boxplot.stats(dati$Lunghezza)$out
Lunghezza_outlier## [1] 390 400 410 405 420 420 360 430 400 410 560 380 420 390 405 360 430 310 390
## [20] 430 390 420 410 420 370 430 430 410 385 390 560 315 420 340 410 430 430 380
## [39] 325 420 420 410 400 430 355 370 410 380 355 410 425 400 370 430 565 405 320
## [58] 345 430Cranio_outlier <- boxplot.stats(dati$Cranio)$out
Cranio_outlier## [1] 298 382 287 273 285 280 390 384 295 276 382 274 289 295 386 390 277 280 272
## [20] 254 297 385 295 275 266 383 390 292 292 293 278 253 277 390 381 270 267 290
## [39] 276 235 294 299 298 290 273 290 265 245AnniM_outlier <- boxplot.stats(dati$Anni.madre)$out
AnniM_outlier## [1] 13 45 43 44 44 43 14 46 1 0 14 44 44Sostituisco il valore anomalo di 0 e 1 anni con la media delle età (calcolata escludendo il valore zero e 1)
media_anni_madre <- mean(dati$Anni.madre[!(dati$Anni.madre %in% c(0, 1))], na.rm = TRUE)
dati$Anni.madre[dati$Anni.madre %in% c(0, 1)] <- media_anni_madreControllo il dataset dopo la sostituzione
summary(dati) ## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. :13.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.19 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
## Distribuzione per Sesso
library(ggplot2)## Warning: il pacchetto 'ggplot2' è stato creato con R versione 4.1.2##
## Caricamento pacchetto: 'ggplot2'## I seguenti oggetti sono mascherati da 'package:psych':
##
## %+%, alphaggplot(dati, aes(x = Sesso, fill = Sesso)) +
geom_bar() +
labs(title = "Distribuzione per Sesso", x = "Sesso", y = "Frequenza") +
scale_fill_manual(values = c("lightblue", "pink")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)) +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size = 5)
Il grafico mostra una quasi equa distribuzione nei due generi. Questo da valore alle conclusioni riguardo le differenze di genere
Distribuzione per Tipo di Parto indipendentemente dagli ospedali
ggplot(dati, aes(x = Tipo.parto, fill = Tipo.parto)) +
geom_bar() +
labs(title = "Distribuzione per Tipo di Parto", x = "Tipo di Parto", y = "Frequenza") +
scale_fill_brewer(palette = "Set2") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)) +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size = 5)
si notano 728 osservazioni per i parti cesarei e 1772 per i parti naturali come mostra il grafico.
Distribuzione per Ospedale
ggplot(dati, aes(x = Ospedale, fill = Ospedale)) +
geom_bar() +
labs(title = "Distribuzione per Ospedale", x = "Ospedale", y = "Frequenza") +
scale_fill_viridis_d() +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)) +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.15, size = 5)
Si nota una quasi pari distribuzione di osservazioni tra i tre ospedali
Distribuzione per Fumatrici
ggplot(dati, aes(x = Fumatrici, fill = Fumatrici)) +
geom_bar() +
labs(title = "Distribuzione per Fumo Materno", x = "Fumatrici", y = "Frequenza") +
scale_fill_manual(values = c("lightgreen", "lightcoral")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)) +
geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size = 5)
Come mostrato dai grafici, il campione dispone di pochi dati riguardo le madri fumatrici
Relazione tra Peso e Settimane di Gestazione
plot(dati$Gestazione, dati$Peso, main = "Peso vs Settimane di Gestazione",
xlab = "Settimane di Gestazione", ylab = "Peso (grammi)", col = "darkblue", pch = 20)
Si può notare un andamento crescente tra settimane di gestazione e peso.
Peso Neonatale rispetto al Fumo Materno o meno
boxplot(dati$Peso ~ dati$Fumatrici, main = "Peso Neonatale e Fumo Materno",
xlab = "Fumo Materno (0 = No, 1 = Si)", ylab = "Peso (grammi)", col = c("lightblue", "lightcoral"))
Nonostante le mediane appaiano molto simili, le madri che non fumano risultano avere bambini con pesi che si disperdono maggiormente oltre i limiti del boxplot (per cui si notano molti outlier). Per cui si nota una maggiore variabilità in questo sottogruppo piuttosto che nelle non fumatrici (anche perchè sono di numerosità molto inferiore).
tabella_parti <- table(dati$Ospedale, dati$Tipo.parto)
df_parti <- as.data.frame(tabella_parti)
colnames(df_parti) <- c("Ospedale", "TipoParto", "Frequenza")
# Creazione del grafico a barre
library(ggplot2)
ggplot(df_parti, aes(x = TipoParto, y = Frequenza, fill = Ospedale)) +
geom_bar(stat = "identity", position = position_dodge()) + # Barre affiancate
geom_text(
aes(label = Frequenza),
position = position_dodge(width = 0.9),
vjust = -0.5,
size = 3
) + # Etichette sopra le barre
scale_fill_brewer(palette = "Set2") + # Tavolozza di colori personalizzata
labs(
title = "Distribuzione dei Tipi di Parto per Ospedale",
x = "Tipo di Parto",
y = "Frequenza",
fill = "Ospedale"
)
La distribuzione appare abbastanza bilanciata tra gli ospedali nel tipo di parto, ma eventuali differenze vanno verificate con un test statistico
valutazioni sul campione in esame rispetto alla popolazione
chisq_test <- chisq.test(tabella_parti)
chisq_test##
## Pearson's Chi-squared test
##
## data: tabella_parti
## X-squared = 1.0972, df = 2, p-value = 0.5778NON Si può rifiutare l’ipotesi nulla di ugualianza. il tipo di parto non è legato all’ospedale
t_test_peso <- t.test(dati$Peso, mu = 3300, alternative = "two.sided")
t_test_peso ##
## One Sample t-test
##
## data: dati$Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081la media del campione non è significativamente diversa da quella della popolazione
t_test_lunghezza <- t.test(dati$Lunghezza, mu = 500, alternative = "two.sided")
t_test_lunghezza ##
## One Sample t-test
##
## data: dati$Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692la media del campione è significativamente diversa da quella della popolazione per la lunghezza anche se quasi sempre entro circa 7mm poichè l’estremo inferiore è 493.6
#test condizionati dal Sesso
t_test_peso_sesso <- t.test(Peso ~ Sesso, data = dati)
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.215t_test_lunghezza_sesso <- t.test(Lunghezza ~ Sesso, data = dati)
t_test_lunghezza_sesso##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.582, df = 2459.3, 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:
## -11.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672t_test_cranio_sesso <- t.test(Cranio ~ Sesso, data = dati)
t_test_cranio_sesso##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4486Sono tutti significativi in quanto p < 0.01 indicando differenze tra i due generi per peso lunghezza e cranio, come previsto dalla letteratura a vantaggio del genere maschile come mostrato dagli output
shapiro.test(dati$Peso)##
## Shapiro-Wilk normality test
##
## data: dati$Peso
## W = 0.97066, p-value < 2.2e-16Poichè il p-value è inferiore a 0.05, rifiutiamo l’ipotesi nulla di normalità i dati relativi al Peso non seguono una distribuzione normale
Conversione in fattori
dati$Ospedale <- as.factor(dati$Ospedale)
dati$Tipo.parto <- as.factor(dati$Tipo.parto)
dati$Sesso <- as.factor(dati$Sesso)
str(dati)## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : num 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 : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
## $ Ospedale : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
## $ Sesso : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...Creazione del Modello di Regressione
Correlazioni
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- (cor(x, y))
txt <- format(c(r, 1), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 1.5)
}
pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)
Tenendo presente lo scopo reale del modello di previsione farei delle considerazioni preliminari: sembra che con il peso incidano le variabili GESTAZIONE LUNGHEZZA, CRANIO, SESSO poichè i coefficienti corr sono rispettivamente r= 0.59 ; 0.8 ; 0.7 e 0.24.
Studi in letteratura indicano che il fumo può incidere negativamente sul peso, lunghezza e parto prematura a causa del minore apporto di ossigeno (fonte: salute.gov) inoltre è¨ noto in letteratura come una madre fumatrice possa innescare o favorire altre patologie al bambino. Per queste considerazioni anche se in questo specifico campione non risulta incisivo l’apporto del fumo si dovrebbero fare considerazioni sul modello. Ad esempio se inserire la variabile può dare maggiore generalizzabilità al modello
Oltre il campione specifico, altra motivazione della non incidenza del fumo potrebbe essere la modalità in cui questa variabile è stata raccolta ovvero che la variabile è stata caratterizzata solo in due livelli. Valutare i mesi o anni di fumo precedenti oppure la quanitità di sigarette medie potrebbe essere maggiormente utile per future ricerche per maggiore prevedibilità del peso Visto che lo scopo è utilizzare al meglio le TIV (terapie intensive) e le risorse umane terrei gli anni della madre in considerazione anche se non sembra esserci correlazione in questo campione
Inoltre, visivamente la relazione PESO - LUNGHEZZA E PESO-CRANIO sembrano lineari, la relazione PESO-GESTAZIONE appare curvilinea per questo si potrebbe provare a valutare tale relazione per una maggiore prevedibilità  (in letteratura è noto che nascere con molte settimane di anticipo compromette in modo maggiore rispetto a differenze più vicine al completamento della gestazione per cui suggerisce una relazione non lineare)
t.test(dati$Peso~dati$Sesso)##
## Welch Two Sample t-test
##
## data: dati$Peso by dati$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.215t.test(dati$Peso~dati$Tipo.parto)##
## Welch Two Sample t-test
##
## data: dati$Peso by dati$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.916il genere come prevedibile incide sul peso di circa 150g mentre il tipo di parto no è significativo per cui può essere escluso dal successivo modello
Anova per valutare eventuali differenze tra tipo di parto e gli Ospedali
anova_result <- aov(Peso ~ Ospedale, data = dati)
summary(anova_result)## Df Sum Sq Mean Sq F value Pr(>F)
## Ospedale 2 936237 468118 1.699 0.183
## Residuals 2497 687952305 275512boxplot(Peso ~ Ospedale, data = dati, pch=20,
main = "Peso per Ospedale",
xlab = "Ospedale", ylab = "Peso",
col = "lightblue")
come si evince dal p value dell’anova 0.18 maggiore di 0.05 e dai boxplot non ci sono associazioni significative tra Peso e Ospedale
Modello di regressione lineare multipla inserendo solo inizialmente tutte le variabili
modello_peso <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
Tipo.parto + Ospedale + Sesso, data = dati)
summary(modello_peso) ##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
##
## 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.1400 141.3974 -47.633 < 2e-16 ***
## Anni.madre 0.7975 1.1463 0.696 0.4867
## N.gravidanze 11.4130 4.6665 2.446 0.0145 *
## Fumatrici -30.1567 27.5396 -1.095 0.2736
## Gestazione 32.5262 3.8179 8.519 < 2e-16 ***
## Lunghezza 10.2951 0.3007 34.237 < 2e-16 ***
## Cranio 10.4725 0.4261 24.580 < 2e-16 ***
## Tipo.partoNat 29.5025 12.0848 2.441 0.0147 *
## Ospedaleosp2 -11.2217 13.4388 -0.835 0.4038
## Ospedaleosp3 28.0985 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-16appare come tutte le variabili tranne anni madre ed essere fumatrici o meno, sono significative
Selezione del modello
library(MASS)## Warning: il pacchetto 'MASS' è stato creato con R versione 4.1.3modello_aic <- stepAIC(modello_peso, direction = "both") ## 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 36321 186809099 28074
## - Fumatrici 1 89979 186862757 28075
## <none> 186772779 28075
## - Tipo.parto 1 447229 187220007 28079
## - N.gravidanze 1 448861 187221640 28079
## - Ospedale 2 686401 187459180 28081
## - Sesso 1 3611594 190384372 28121
## - Gestazione 1 5446472 192219251 28145
## - Cranio 1 45338176 232110954 28617
## - Lunghezza 1 87959836 274732615 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 36321 186772779 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 37239 186862757 28075
## - Tipo.parto 1 440684 187340680 28077
## - Ospedale 2 701680 187601677 28079
## - 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 29039Calcolo in automatico il modello migliore partendo da quello completo
summary(modello_aic)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = dati)
##
## 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-16con il metodo AIC si ottiene il modello composto N.gravidanze Gestazione Lunghezza Cranio Tipo.partoNat Ospedale Sesso come variabili esplicative del peso R quadro aggiustato risualta 0.7265 Gli anni madre e fumatrici sono state elimiante dal modello migliore come ci aspettavamo dalle considerazioni sui p value
n <- nrow(dati)
modello_bic <- stepAIC(modello_peso, direction = "both", k = log(n))## Start: AIC=28139.46
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36321 186809099 28132
## - Fumatrici 1 89979 186862757 28133
## - Ospedale 2 686401 187459180 28133
## - Tipo.parto 1 447229 187220007 28138
## - N.gravidanze 1 448861 187221640 28138
## <none> 186772779 28140
## - Sesso 1 3611594 190384372 28180
## - Gestazione 1 5446472 192219251 28204
## - Cranio 1 45338176 232110954 28675
## - Lunghezza 1 87959836 274732615 29096
##
## Step: AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28126
## - Ospedale 2 692738 187501837 28126
## - Tipo.parto 1 448222 187257321 28130
## <none> 186809099 28132
## - N.gravidanze 1 633756 187442855 28133
## + Anni.madre 1 36321 186772779 28140
## - Sesso 1 3618736 190427835 28172
## - Gestazione 1 5412879 192221978 28196
## - Cranio 1 45588236 232397335 28670
## - Lunghezza 1 87950050 274759149 29089
##
## Step: AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 701680 187601677 28119
## - Tipo.parto 1 440684 187340680 28124
## <none> 186899996 28126
## - N.gravidanze 1 610840 187510837 28126
## + Fumatrici 1 90897 186809099 28132
## + Anni.madre 1 37239 186862757 28133
## - Sesso 1 3602797 190502794 28165
## - Gestazione 1 5346781 192246777 28188
## - Cranio 1 45632149 232532146 28664
## - Lunghezza 1 88355030 275255027 29086
##
## Step: AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Ospedale 2 701680 186899996 28126
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 43769 187557908 28127
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto 1 463870 187601677 28119
## + Ospedale 2 724866 187340680 28124
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 44972 188020574 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(modello_bic)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16con il metodo BIC si ottiene il modello composto N.gravidanze Gestazione Lunghezza Cranio Sesso come variabili esplicative del peso R quadro aggiustato risualta 0.7265 (spiega il 72.65% della varianza) identico al metodo AIC anche come p value, dalle considerazioni fatte precedentemente questo metodo sembra dare il risultato più simile a ciò che ci si aspetta dalla letteratura , tranne per il fumo che è stato nuovamente eliminato. Forse come spiegato sopra per il campione o per il metodo dicotomico di valutarlo
Rispetto al modello preferito dal metodo AIC è stato eliminato anche l’ospedale. In effetti logicamente questa variabile non dovrebbe contribuire al peso. la significatività statistca trovata per l’ospedale 3 (p =0.034) può essere una seppur anomala, fluttiazione statistica. Nel contesto reale potrebbe accadere che alcuni ospedali siano preferiti ad altri nei casi in cui si rischia un parto anticipato e ciò può indirizzare la scelta senza significare che abbiano una incidenza sul peso
Riprendo tale modello indagando possibili interazioni
modello_interazioni <- lm(Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio)^2, data = dati)
summary(modello_interazioni)##
## Call:
## lm(formula = Peso ~ (Anni.madre + N.gravidanze + Fumatrici +
## Gestazione + Lunghezza + Cranio)^2, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1021.04 -177.20 -16.37 164.88 2317.73
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3034.43888 1417.96167 2.140 0.032452 *
## Anni.madre -56.29135 26.75357 -2.104 0.035473 *
## N.gravidanze -126.28224 95.68228 -1.320 0.187021
## Fumatrici 102.22367 916.32621 0.112 0.911183
## Gestazione -282.09085 66.66166 -4.232 2.40e-05 ***
## Lunghezza 19.03932 5.44024 3.500 0.000474 ***
## Cranio -20.69010 8.21599 -2.518 0.011856 *
## Anni.madre:N.gravidanze -0.17443 0.85893 -0.203 0.839092
## Anni.madre:Fumatrici 2.16084 5.75098 0.376 0.707147
## Anni.madre:Gestazione 2.75425 0.75999 3.624 0.000296 ***
## Anni.madre:Lunghezza -0.24395 0.06026 -4.048 5.32e-05 ***
## Anni.madre:Cranio 0.20844 0.08524 2.445 0.014543 *
## N.gravidanze:Fumatrici -9.63511 24.17947 -0.398 0.690309
## N.gravidanze:Gestazione -4.92195 2.99569 -1.643 0.100508
## N.gravidanze:Lunghezza 0.39762 0.24595 1.617 0.106075
## N.gravidanze:Cranio 0.41698 0.37188 1.121 0.262284
## Fumatrici:Gestazione -53.44547 22.36459 -2.390 0.016935 *
## Fumatrici:Lunghezza 5.72437 1.71406 3.340 0.000851 ***
## Fumatrici:Cranio -2.64128 2.41281 -1.095 0.273758
## Gestazione:Lunghezza 0.00801 0.11039 0.073 0.942158
## Gestazione:Cranio 0.73650 0.21549 3.418 0.000642 ***
## Lunghezza:Cranio -0.00661 0.01405 -0.471 0.637968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.3 on 2478 degrees of freedom
## Multiple R-squared: 0.7314, Adjusted R-squared: 0.7291
## F-statistic: 321.3 on 21 and 2478 DF, p-value: < 2.2e-16in combinazione con la letteratura, delle interazioni significative potremmo considerare Anni.madre:Gestazione (2.7) Anni.madre:Lunghezza (-0.24) Fumatrici:Gestazione (-53.45) Fumatrici:Lunghezza (5.7) Gestazione:Cranio (+0.74) però è possibile notare che seppur significativi i contributi sono di ordine del grammo o decimo di grammo tranne per Fumatrici:gestazione per cui trascurabili l’essere fumatrice modifica significativamente l’effetto della gestazione sul peso del bambino con -53.45 g p <0.05 quindi questa è forse l’unica interazione da considerare
Modello con trasformazioni non lineari per la gestazione come suggerito dal grafico precedente
modello_nonlineare <- lm(Peso ~ poly(Gestazione, 2) + Lunghezza + Anni.madre + Fumatrici, data = dati)
summary(modello_nonlineare)##
## Call:
## lm(formula = Peso ~ poly(Gestazione, 2) + Lunghezza + Anni.madre +
## Fumatrici, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1177.9 -204.2 -19.3 194.5 3598.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3710.7084 155.3852 -23.881 < 2e-16 ***
## poly(Gestazione, 2)1 4388.5893 401.6457 10.927 < 2e-16 ***
## poly(Gestazione, 2)2 152.0372 317.2673 0.479 0.631832
## Lunghezza 13.8878 0.3072 45.204 < 2e-16 ***
## Anni.madre 4.4601 1.2022 3.710 0.000212 ***
## Fumatrici -26.4350 31.1363 -0.849 0.395958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 310.2 on 2494 degrees of freedom
## Multiple R-squared: 0.6516, Adjusted R-squared: 0.6509
## F-statistic: 933 on 5 and 2494 DF, p-value: < 2.2e-16Il termine quadratico dela Gestazione non è significativo e potrebbe essere rimosso per semplificare il modello.
Confronto tra modelli - modello_peso, modello_peso_ridotto, modello_peso_int
modello_peso_ridotto <- lm(Peso ~ N.gravidanze + Gestazione +Lunghezza + Cranio + Sesso, data = dati)
summary(modello_peso_ridotto) ##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16modello_peso_int <- lm(Peso ~ N.gravidanze + Gestazione +Lunghezza + Cranio + Sesso + Fumatrici:Gestazione, data = dati)
summary(modello_peso_int) ##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Fumatrici:Gestazione, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.33 -181.33 -15.97 162.95 2636.22
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6682.3544 135.7187 -49.237 < 2e-16 ***
## N.gravidanze 12.7287 4.3450 2.929 0.00343 **
## Gestazione 32.6216 3.8062 8.571 < 2e-16 ***
## Lunghezza 10.2334 0.3009 34.008 < 2e-16 ***
## Cranio 10.5355 0.4262 24.717 < 2e-16 ***
## SessoM 78.1998 11.2029 6.980 3.76e-12 ***
## Gestazione:Fumatrici -0.8028 0.7024 -1.143 0.25316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265
## F-statistic: 1107 on 6 and 2493 DF, p-value: < 2.2e-16inserito nel modello ridotto Gestazione:Fumatrici risulta non significativo per cui eliminabile
AIC(modello_peso, modello_peso_ridotto, modello_peso_int) ## df AIC
## modello_peso 12 35172.09
## modello_peso_ridotto 7 35179.33
## modello_peso_int 8 35180.02BIC(modello_peso, modello_peso_ridotto, modello_peso_int)## df BIC
## modello_peso 12 35241.97
## modello_peso_ridotto 7 35220.10
## modello_peso_int 8 35226.61il metodo BIC che predilige modelli piu semplici conferma che che il modello migliore è modello_peso_ridotto ci sarebbe da decidere se includere nel modello di previsione il fumo (e/o sue interazioni) o meno per dare maggiore generalizzabilità  anche se questi dati non suggeriscono che vada inserita tale variabile. Questa considerazione dovrebbe essere fatta con un confronto adeguato con esperti del settore per valutazioni cliniche. Per quanto mostrato nei dati non dovrebbe essere inserito
per cui il miglior modello appare quello che associa al Peso: N.gravidanze, Gestazione,Lunghezza ,Cranio,Sesso
scelto il modello valutiamo le altre info sui residui
analisi residui
par(mfrow=c(2,2))
plot(modello_peso_ridotto)
I grafici mostrano un buon adattamento tranne verso i valori più estremi. Viene visualizzata la distanza di cook e sembra indicare che il valore 1551 peggiori il modello lo stesso valore come è evidente dai residui e dalla distribuzione dei quartili appare anomalo
verifichiamo con i test
lmtest::bptest(modello_peso_ridotto) ##
## studentized Breusch-Pagan test
##
## data: modello_peso_ridotto
## BP = 90.253, df = 5, p-value < 2.2e-16Poichè il p-value è < 0.01, rifiutiamo l’ipotesi nulla di omoschedasticitÃ
lmtest::dwtest(modello_peso_ridotto) ##
## Durbin-Watson test
##
## data: modello_peso_ridotto
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0I residui non presentano una correlazione significativa tra di loro e soddisfano l’assunzione di indipendenza. Non possiamo rifiutare l’ipotesi nulla di autocorrelazione dei residui
shapiro.test(modello_peso_ridotto$residuals) ##
## Shapiro-Wilk normality test
##
## data: modello_peso_ridotto$residuals
## W = 0.97408, p-value < 2.2e-16rifiutiamo l’ipotesi nulla di normalità dei residui. I residui non seguono una distribuzione normale, il che potrebbe indicare che il modello lineare non è perfettamente appropriato (ed esempio effetti non lineari) o che nei dati vi siano outlier. plot(density(residuals(modello_peso_ridotto))) Ci aspettiamo che l’osservazione 1551 sia oltre la soglia dello 0.5 dell’indice di Cook
Verifica dell’indice di Cook, Leverage e Outlier
Leverage: quanto un’osservazione è distante dallo spazio medio delle variabili indipendenti. Un’osservazione con alta leverage può influenzare notevolmente i coefficienti del modello.
lev <- hatvalues(modello_peso_ridotto)
plot(lev)
p <- sum(lev)
n <- length(lev)
soglia <- 2 * p / n
abline(h = soglia, col = 2) 
lev[lev > soglia]## 13 15 34 67 89 96
## 0.005630918 0.007050422 0.006744143 0.005891590 0.012816470 0.005351586
## 101 106 131 134 151 155
## 0.007526751 0.014479871 0.007229339 0.007552819 0.010883937 0.007207682
## 161 189 190 204 205 206
## 0.020335483 0.004893297 0.005366557 0.014490604 0.005351634 0.009476652
## 220 294 305 310 312 315
## 0.007393997 0.005912765 0.005442061 0.028812123 0.013169272 0.005385800
## 378 440 442 445 486 492
## 0.015934061 0.005404736 0.007723662 0.007509382 0.005164446 0.008274018
## 497 516 582 587 592 614
## 0.005166306 0.013079851 0.011665555 0.008412325 0.006384116 0.005299262
## 638 656 657 684 697 702
## 0.006688287 0.005927777 0.005322685 0.008818987 0.005863826 0.005202259
## 729 748 750 757 765 805
## 0.005023115 0.008565543 0.006942097 0.008145491 0.006070298 0.014356657
## 828 893 895 913 928 946
## 0.007179817 0.005075205 0.005295896 0.005571144 0.022742332 0.006909044
## 947 956 985 1008 1014 1049
## 0.008409465 0.007784123 0.007039416 0.005343037 0.008470133 0.004956169
## 1067 1091 1106 1130 1166 1181
## 0.008465430 0.008933360 0.005967317 0.031728597 0.005513559 0.005677676
## 1188 1200 1219 1238 1248 1273
## 0.006477203 0.005492370 0.030694311 0.005908078 0.014622914 0.007085831
## 1291 1293 1311 1321 1325 1356
## 0.006117497 0.006073639 0.009625908 0.009293111 0.004857169 0.005303442
## 1357 1385 1395 1400 1402 1411
## 0.006965051 0.012636943 0.005126697 0.005925069 0.004811441 0.008048184
## 1420 1428 1429 1450 1505 1551
## 0.005155654 0.008192811 0.021757172 0.015104831 0.013330439 0.048769569
## 1553 1556 1573 1593 1606 1610
## 0.008504889 0.005919673 0.005047204 0.005623758 0.005001812 0.008722184
## 1617 1619 1628 1686 1693 1701
## 0.004866796 0.015067498 0.005069731 0.009349313 0.005077858 0.010842957
## 1712 1718 1727 1735 1780 1781
## 0.006992084 0.006958857 0.013300523 0.004884846 0.025538678 0.016832361
## 1809 1827 1868 1892 1962 1967
## 0.008707504 0.006065698 0.005205637 0.005332812 0.005540442 0.005337356
## 1977 2037 2040 2046 2086 2089
## 0.006927281 0.004889127 0.011494872 0.005471670 0.013193090 0.006293550
## 2098 2114 2115 2120 2140 2146
## 0.005094455 0.013316875 0.011772090 0.018659995 0.006244232 0.005802168
## 2148 2149 2157 2175 2200 2215
## 0.007926839 0.013583436 0.005907225 0.032527273 0.011670024 0.004892265
## 2216 2220 2221 2224 2225 2244
## 0.008117864 0.005414040 0.021628717 0.005838076 0.005591261 0.006929217
## 2257 2307 2317 2318 2337 2359
## 0.006170254 0.013965608 0.007673614 0.004831118 0.005230450 0.010067364
## 2408 2422 2436 2437 2452 2458
## 0.009696691 0.021532808 0.004986522 0.023943328 0.023838489 0.008506087
## 2471 2478
## 0.020903740 0.005775173La soglia individua molti dati oltre; tuttavia, leverage alto non implica necessariamente che l’osservazione sia problematica.
Analisi degli outlier nei residui standardizzati
Gli outlier nei residui standardizzati mostrano osservazioni che deviano significativamente dal modello.
plot(rstudent(modello_peso_ridotto))
abline(h = c(-2, 2))
car::outlierTest(modello_peso_ridotto)## rstudent unadjusted p-value Bonferroni p
## 1551 10.051908 2.4906e-23 6.2265e-20
## 155 5.027798 5.3138e-07 1.3285e-03
## 1306 4.827238 1.4681e-06 3.6702e-03Valutazione della distanza di Cook Combina gli effetti di leverage e outlier per valutare se il dato è anomalo rispetto al modello.
cook <- cooks.distance(modello_peso_ridotto)
plot(cook, ylim = c(0, 1))
abline(h = 0.5, col = "red") # Solo un'osservazione è oltre 0.5.
osservazioni_influenti <- which(cook > 0.5)
print(osservazioni_influenti) # Mostra gli indici delle osservazioni influenti.## 1551
## 1551Rimozione dell’osservazione influente 1551 Verifica che il database abbia il numero corretto di dati dopo la modifica.
dati_modificati <- dati[-osservazioni_influenti, ]
nrow(dati) - nrow(dati_modificati)## [1] 1Creazione del nuovo modello con dati modificati (eliminando l’osservazione 1551) Valutazione delle prestazioni del modello aggiornato.
modello_peso_ridotto_mod <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso, data = dati_modificati)
summary(modello_peso_ridotto_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati_modificati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1165.74 -179.59 -12.74 162.89 1410.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6683.4142 133.0802 -50.221 < 2e-16 ***
## N.gravidanze 13.1652 4.2557 3.094 0.002 **
## Gestazione 29.5891 3.7340 7.924 3.43e-15 ***
## Lunghezza 10.8927 0.3017 36.109 < 2e-16 ***
## Cranio 9.9187 0.4225 23.476 < 2e-16 ***
## SessoM 78.1348 10.9840 7.114 1.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7372, Adjusted R-squared: 0.7367
## F-statistic: 1399 on 5 and 2493 DF, p-value: < 2.2e-16AIC(modello_peso_ridotto, modello_peso_ridotto_mod)## Warning in AIC.default(modello_peso_ridotto, modello_peso_ridotto_mod): models
## are not all fitted to the same number of observations## df AIC
## modello_peso_ridotto 7 35179.33
## modello_peso_ridotto_mod 7 35066.98BIC(modello_peso_ridotto, modello_peso_ridotto_mod)## Warning in BIC.default(modello_peso_ridotto, modello_peso_ridotto_mod): models
## are not all fitted to the same number of observations## df BIC
## modello_peso_ridotto 7 35220.10
## modello_peso_ridotto_mod 7 35107.74Ripetizione dei test diagnostici sul nuovo modello Verifica di omoschedasticità , autocorrelazione e normalità dei residui.
par(mfrow = c(2, 2))
plot(modello_peso_ridotto_mod)
lmtest::bptest(modello_peso_ridotto_mod)##
## studentized Breusch-Pagan test
##
## data: modello_peso_ridotto_mod
## BP = 11.393, df = 5, p-value = 0.04411lmtest::dwtest(modello_peso_ridotto_mod)##
## Durbin-Watson test
##
## data: modello_peso_ridotto_mod
## DW = 1.954, p-value = 0.1251
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(modello_peso_ridotto_mod$residuals)##
## Shapiro-Wilk normality test
##
## data: modello_peso_ridotto_mod$residuals
## W = 0.98886, p-value = 4.764e-13Breusch-Pagan Test: p-value=0.04411. Il p-value inferiore a 0.05 suggerisce eteroschedasticità , cioè una variabilità non costante dei residui. Questo viola l’assunzione di omoschedasticità del modello di regressione Anche se dopo l’eliminazione del valore outlier non possiamo rifiutare l’ipotesi nulla al 0.01, migliorando le carattetistiche del modello Durbin-Watson Test: Il p-value maggiore di 0.05 indica che non c’è evidenza significativa di autocorrelazione positiva tra i residui del modello. Shapiro-Wilk Test: Il p-value molto basso indica che i residui non seguono una distribuzione normale, un’altra violazione delle assunzioni classiche del modello.
Calcolo di R^2 e RMSE per il nuovo modello
pred <- predict(modello_peso_ridotto_mod, newdata = dati_modificati)
actual <- dati_modificati$Peso
R2 <- round(cor(pred, actual)^2, 3)
RMSE <- round(sqrt(mean((pred - actual)^2)), 2)
R2## [1] 0.737RMSE## [1] 268.93Il modello spiega il 73.7% della varianza dei dati. 1% in più rispetto al modello precedente (72.65%)
Previsione del peso per un nuovo neonato Peso previsto sulla base delle caratteristiche specificate.
nuovo_neonato <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = 500,
Cranio = 340,
Sesso = "F"
)
peso_previsto <- predict(modello_peso_ridotto_mod, newdata = nuovo_neonato)Analisi dei coefficienti del modello Interpretazione dei coefficienti utili al modello di predizione.
coeff <- summary(modello_peso_ridotto_mod)$coefficients
round(coeff["N.gravidanze", ], 2) ## Estimate Std. Error t value Pr(>|t|)
## 13.17 4.26 3.09 0.00round(coeff["Gestazione", ], 2) ## Estimate Std. Error t value Pr(>|t|)
## 29.59 3.73 7.92 0.00round(coeff["Lunghezza", ], 2) ## Estimate Std. Error t value Pr(>|t|)
## 10.89 0.30 36.11 0.00round(coeff["Cranio", ], 2) ## Estimate Std. Error t value Pr(>|t|)
## 9.92 0.42 23.48 0.00round(coeff["SessoM", ], 2) ## Estimate Std. Error t value Pr(>|t|)
## 78.13 10.98 7.11 0.00L’aumento unitario del numero di gravidanze tende ad aumentare di 13.17 grammi il peso del neonato L’aumento unitario delle settimane di gestazione tende ad aumentare di 29.59 grammi il peso del neonato L’aumento unitario della lunghezza del bambino tende ad aumentare di 10.89 grammi il peso del neonato L’aumento unitario della circonferenza del cranio del bambino tende ad aumentare di 9.92 grammi il peso del neonato Il sesso maschile del bambino tende ad aumentare di 78.13 grammi il peso del neonato
Visualizzazione: scatterplot con linea di regressione Relazione tra settimane di gestazione e peso neonatale.
ggplot(dati, aes(x = Gestazione, y = Peso, color = factor(Sesso))) +
geom_point(position = position_jitter(width = 0.2, height = 0)) +
geom_smooth(method = "lm", color = "blue") +
labs(
title = "Relazione tra Settimane di Gestazione e Peso Neonatale",
x = "Settimane di Gestazione",
y = "Peso Neonatale (g)",
color = "Sesso"
) +
theme_minimal()## `geom_smooth()` using formula 'y ~ x'
Visualizzazione: grafico a violino per gli ospedali Distribuzione del peso neonatale nei tre ospedali.
ggplot(dati, aes(x = Ospedale, y = Peso, fill = Ospedale)) +
geom_violin(alpha = 0.7) +
labs(
title = "Distribuzione del Peso Neonatale nei Tre Ospedali",
x = "Ospedale",
y = "Peso Neonatale (g)"
) +
theme_minimal()
Si possono notare le similarità tra le frequenze e le distribuzioni del Peso dei bambini nei tre Ospedali
Visualizzazione: scatterplot per cranio e peso Relazione tra circonferenza cranio e peso.
ggplot(data = dati) +
geom_point(aes(x = Cranio, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Cranio, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Circonferenza Cranio e Peso (g)",
x = "Circonferenza Cranio",
y = "Peso",
fill = "Peso (g)"
)## `geom_smooth()` using formula 'y ~ x'
La circonferenza dei maschi appare mediamente sempre maggiore rispetto a quella delle femmine come indicato dalle rette.
Visualizzazione: scatterplot per lunghezza e peso Relazione tra lunghezza e peso.
ggplot(data = dati) +
geom_point(aes(x = Lunghezza, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Lunghezza, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Lunghezza (cm) e Peso (g)",
x = "Lunghezza (cm)",
y = "Peso",
fill = "Peso (g)"
)## `geom_smooth()` using formula 'y ~ x'
Notiamo che le due rette hanno un punto di intersezione e il genere femminile man mano che si avvicina il momento del parto tendono ad avere una lunghezza inferiore ai maschi.
Per concludere possiamo dire che il modello (modello_peso_ridotto) trovato con i predittori N.gravidanze, Gestazione Lunghezza, Cranio, Sesso possa considerarsi sufficientemente buono per predire il peso del neonato (R aggiustato = 0.736).