Progetto_finale_Neonati
Marco Sablone
2024-12-03
1) ANALISI PRELIMINARE
dati = read.csv("neonati.csv", stringsAsFactors = T)
attach(dati)
head(dati)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 Fdim(dati)## [1] 2500 10n = nrow(dati)Il dataset presenta 10 variabili con 2500 osservazioni. In relazione alle consegne della presente ricerca, assumeremo la variabile “Peso” come variabile dipendente. Eseguiamo una prima analisi su questa variabile
plot(density(Peso),
main = "Distribuzione di densità del peso dei bambini alla nascita",
xlab = "Peso del bambino in grammi")
summary(Peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930boxplot(Peso)
La variabile è una variabile di tipo quantitatva continua espressa in
grammi il cui range di variazione è determinato da un minimo di 830 g ed
un massimo di 4930 g. In media il campione assume un peso di 3284 g. Da
una prima analisi grafica, il campione della variaible Peso sembra
distribuirsi secondo una distribuzione standard denotando un’asimmetria
leggermente negativa ed una “campana” molto accentuata. inoltre è
possibile notare come la funzione di densità presenti un notevole
strascico della coda sinistra: comparando la funzione con il relativo
boxplot è possibile evidenziare la presenza di molti outliers sopratutto
al disotto del baffo inferiore.
Cerchiamo di trovare conferma di quanto appena descritto con i dati
moments::skewness(Peso)## [1] -0.6470308moments::kurtosis(Peso) - 3## [1] 2.031532L’indice skewness denota una dstribuzione con asimmetria negativa, facendo si che la distribuzione si concentri sulle osservazioni con valore più alto. Inoltre l’indice di Curtosi definisce una curva leptocurtica, ovvero molto accentuata, come è possibile verificare dal grafico della distribuzione di densità.
library(ggplot2)
ggplot(data=dati)+
geom_histogram(aes(x=Peso),
fill = "green",
col ="blue") +
labs(title = "Peso del bambino alla nascita",
x = "Peso del bambino in grammi") +
scale_x_continuous(breaks = seq(500,5000,500))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Ad ulteriore conferma di quanto esposto in precedenza, la distribuzione presenta quel che potremo definire casi particolari, con neonati il cui peso è inferiore alla media per motivi che dovremo indagare.
VARIABILI INDIPENDENTI Eseguiamo ora una rapida disamina deiregressori presenti nel dataset
Anni.madre
summary(Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00ggplot(data = dati)+
geom_bar(aes(x=Anni.madre),
fill = "pink",
col ="red") +
labs(title = "Età delle madri partorienti",
x = "Età delle madri") +
scale_x_continuous(breaks = seq(0,50,5))
Variabile quantitativa continua che esprime l’età delle madri al momento del parto. L’istogramma descrive che la maggior parte delle partorienti rientrano in un’età compresa tra i 25 e 30 anni. Si notano anche dei valori anomali in quanto il valore minimo risulta essere pari a zero
table(Anni.madre)## Anni.madre
## 0 1 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 1 1 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 172 174 200
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 147 159 110 96 66 64 41 38 27 19 13 8 2 4 1 1Dalla distribuzione di frequenza assoluta, risultano due madri con età di 0 e 1 anno. Evidentemente sono errori, pertanto imputiamo i valori alla mediana della distribuzione
Anni.madre[Anni.madre<=1] = median(Anni.madre)table(Anni.madre)## Anni.madre
## 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
## 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 174 174 200 147 159
## 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 110 96 66 64 41 38 27 19 13 8 2 4 1 1N.gravidanze Variabile quantitativa discreta indicante il numero di gravidanze precedenti alla presente analisi
summary(N.gravidanze)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.9812 1.0000 12.0000hist(N.gravidanze,
col = "red",
main = "Numero gravidanze precedenti",
xlab = "Numero gravidanze")
(table(N.gravidanze) / length(N.gravidanze))*100## N.gravidanze
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 43.84 32.72 13.60 6.00 1.92 0.84 0.44 0.04 0.32 0.08 0.12 0.04 0.04table(N.gravidanze)## N.gravidanze
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1096 818 340 150 48 21 11 1 8 2 3 1 1Possiamo verificare come quasi la metà del campione considerato sia alla prima esperienza neonatale, mentre circa il 45% del campione ha avuto già uno o ue parti. Il resto del campione presenta delle frequenze molto basse e relative a più di due parti, con una punta massima di 12 verficatosi in un solo caso.
Fumatrici Variabile qualitativa nominale dicotomica che divide il campione in due tra madri non fumatrici (valore = 0) e fumatrici (valore = 1)
table(Fumatrici)## Fumatrici
## 0 1
## 2396 104table(Fumatrici) / length(Fumatrici) *100## Fumatrici
## 0 1
## 95.84 4.16ggplot(data = dati, aes(x=Fumatrici))+
geom_bar(col = "black",
fill = c("lightgreen", "red")) +
labs(title = "Classificazione madri fumatrici e non",
x = "Non fumatrici (VERDE) / Fumatrici (ROSSO)")
Il capmione vede la presenza del 95% di madri non fumatrici
Gestazione
Variabile quantitativa continua indicante il numero di settimane di gestazione
table(Gestazione)## Gestazione
## 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
## 1 1 2 4 3 5 8 9 18 16 33 62 192 437 581 741 329 56 2summary(Gestazione)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 38.00 39.00 38.98 40.00 43.00hist(Gestazione,
col = "green",
main = "Numero di settimane di gestazione",
xlab = "Numero settimane di gestazione")
Come era prevedibile attendersi, il numero di settimane di gestazione si
concentra tra le 38 e le 41 settimane, pur presentando casi sporadici di
parti prematuri con un minimo di 25 settimane
LUNGHEZZA E CRANIO Esamineremo insieme le due variabili in quanto si ipotizzano altamente correlate tra loro: a meno che nel neonato non sia insorta un qualche tipo di menomazione, la distribuzione dei due caratteri dovrebbe risutare presocchè simile.
summary(Lunghezza)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0summary(Cranio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390ggplot(dati,
aes(Lunghezza)) +
geom_histogram(fill = "orange",
col = "black",) +
labs(title = "Altezza del neonato",
x = "Altezza in cm") +
scale_x_continuous(breaks = seq(300, 600, 50))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(dati,
aes(Cranio)) +
geom_histogram(fill = "lightblue",
col = "blue",) +
labs(title = "Misura della circonferenza del cranio dei neonati",
x = "Circonferenza in cm") +
scale_x_continuous(breaks = seq(200, 400, 20))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
TIPO.PARTO Variabile qualitativa nominale trasformata in variabile dummy e indicante la tipologia di parto avvenuta, ovvero: naturale (indicato con valore = 1) o con intervento cesareo (indicato con valore = 2)
table(Tipo.parto)## Tipo.parto
## Ces Nat
## 728 1772table(Tipo.parto) / length(Tipo.parto) *100## Tipo.parto
## Ces Nat
## 29.12 70.88ggplot(dati,
aes(x = Tipo.parto)) +
geom_histogram( stat = "count",
fill = c("orange", "yellow"),
col = "black") +
labs(title = "Tipologia di parto",
x = "Tipo parto") ## Warning in geom_histogram(stat = "count", fill = c("orange", "yellow"), :
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
Il campione presenta un 70% di madri che hanno avuto un aprto naturale ed un 30% che hanno avuto bisogno di ricorrere ad un intervento cesareo.
SESSO Infine analizziamo la variabile sesso, ovvero variabile qualitativa nominale dicotomica indicante il sesso del neonato, anche questa trasformata in variabile dummy: maschio (M=1), femmina (F=2)
table(Sesso)## Sesso
## F M
## 1256 1244table(Sesso) / length(Sesso) * 100## Sesso
## F M
## 50.24 49.76ggplot(dati,
aes(x = Sesso))+
geom_histogram(stat = "count",
fill = c("blue", "pink"),
col = "black") +
ggtitle("Sesso dei neonati") ## Warning in geom_histogram(stat = "count", fill = c("blue", "pink"), col =
## "black"): Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
Il campione vede una sostanziale parità tra i nascituri a livello di sesso.
library(corrplot)## corrplot 0.95 loadedlibrary(dplyr)##
## Caricamento pacchetto: 'dplyr'## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionmatrix_corr = dati %>%
dplyr::select(Anni.madre, N.gravidanze, Fumatrici, Gestazione, Peso,
Lunghezza, Cranio) %>%
cor()
matrix_corr %>% round(2)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## Anni.madre 1.00 0.38 0.01 -0.14 -0.02 -0.06
## N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06
## Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02
## Gestazione -0.14 -0.10 0.03 1.00 0.59 0.62
## Peso -0.02 0.00 -0.02 0.59 1.00 0.80
## Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00
## Cranio 0.02 0.04 -0.01 0.46 0.70 0.60
## Cranio
## Anni.madre 0.02
## N.gravidanze 0.04
## Fumatrici -0.01
## Gestazione 0.46
## Peso 0.70
## Lunghezza 0.60
## Cranio 1.00matrix_corr %>% corrplot()
Considerazioni sulla correlation matrix per quanto concerne le sole
varaibili quantitative:
Sulla variabile dipendente Peso possiamo eseguire due importanti osservazioni: la prima è che essa presenta una forte e logica correlazione positiva con le variabili Gestazione, Lunghezza e Cranio; la seconda è inerente, invece, ad un’assenza di correlazione con variabili come Fumatrici, tipologi di parto e numero di gravidanze;
La variabile Gestazione presenta una correlazione negativa con le variabili anni madre e numero di gravidanze: seppur il valore assoluto della correlazione non sia elevata, sembra che l’età della madre e l’aver avuto già altre gravidanze influiscano negativamente sul numero di settimane di gestazione;
La variabile Fumatrici risulta decorrelata da tutte le variabili del dataset
1.1.a) Saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei
ggplot(dati)+
geom_bar(aes(x = Ospedale, fill = Tipo.parto), position = "dodge")+
scale_y_continuous(breaks = seq(0,600,50))+
labs(title = "Tipologie di parto eseguite nei tre ospedali")
Da una prima analisi grafica non sembrerebbe esserci una grande differenza nei numeri di interventi cesarei effettuati nei tre ospedali. Per confermare tali ipotesi eseguiamo un test del chi-quadro tra le variabili categoriche Ospedale e Tipo.parto.
Per prima cosa eseguiamo una tabella di contingenza con le modalità delle due variabili qualitative
tabella = table(Ospedale, Tipo.parto)
tabella## Tipo.parto
## Ospedale Ces Nat
## osp1 242 574
## osp2 254 595
## osp3 232 603eseguiamo il test del chi-quadrato e saggiamo l’ipotesi di indipendenza tra le due variabili
x2_test = chisq.test(tabella)
x2_test##
## Pearson's Chi-squared test
##
## data: tabella
## X-squared = 1.0972, df = 2, p-value = 0.5778x2_test$expected## Tipo.parto
## Ospedale Ces Nat
## osp1 237.6192 578.3808
## osp2 247.2288 601.7712
## osp3 243.1520 591.8480il p-value è abbondandemente al di sopra del 5% e pertanto siamo portati a non rifiutare l’ipotesi nulla di indipendenza tra le due variabili ovvero i cesarei non dipendono dal tipo di ospedale e anche dall’analisi grafica possiamo verificare quanto appena affermato. Inoltre una verifica delle frequesnze attese, mostrano tutti valori superiori a 5 e ciò ci consente di affermare la bontà della statistica test e che essa non induca a bias.
1.1.b) La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quella della popolazione
Considerando che non conosciamo la varianza della popolazione, eseguiamo un t-test per saggiare l’ipotesi che le medie campionarie di Peso e Altezza siano uguali a quella della popolazione
t.test(Peso)##
## One Sample t-test
##
## data: Peso
## t = 312.75, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(Lunghezza)##
## One Sample t-test
##
## data: Lunghezza
## t = 939.81, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692In entrambi i casi si rifiuta l’ipotesi nulla di ugualianza delle medie campionarie rispetto a quelle della popolazione per i caratteri di peso e lunghezza.
1.1.c) Le misure antropometriche sono significativamente diverse tra i due sessi
Selezioniamo quel che possiamo considerare le variabili antropometriche del campione: Peso, Lunghezza, Cranio. Cerchiamo di saggiare l’ipotesi che vi sia una relazione significativa tra questi caratteri antropomorfici e il sesso dei neonati.
ggplot(dati)+
geom_histogram(aes(x = Peso,
fill = Sesso))+
labs(title = "Confronto del peso tra maschi e femmine") +
scale_x_continuous(breaks = seq(0,5500,500)) ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(dati)+
geom_histogram(aes(x = Lunghezza,
fill = Sesso))+
labs(title = "Confronto dell'altezza tra maschi e femmine") +
scale_x_continuous(breaks = seq(0,5500,500)) ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(dati)+
geom_histogram(aes(x = Cranio,
fill = Sesso))+
labs(title = "Confronto della circonferenza cranio tra maschi e femmine") +
scale_x_continuous(breaks = seq(0,5500,500)) ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Eseguiamo un t-test per confronti multipli tra medie a due a due
pairwise.t.test(Peso, Sesso,
paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Peso and Sesso
##
## F
## M <2e-16
##
## P value adjustment method: bonferroniconsiderando che il p-value è nettamente inferiore al 5%, rifiutiamo l’ipotesi nulla e affermiamo che il peso sia significativamente diverso secondo il sesso del neonato, cosi come è possibile osservare dalla precedente analisi grafica. Inoltre dalla correlation matrix abbiamo osservato come i caratteri antropomorfici del campione sia strettamente correlati positivamente, pertanto dovremmo attendere un medesimo risultato nei test con la variabile cranio e lunghezza
pairwise.t.test(Lunghezza, Sesso,
paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Lunghezza and Sesso
##
## F
## M <2e-16
##
## P value adjustment method: bonferronipairwise.t.test(Cranio, Sesso,
paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Cranio and Sesso
##
## F
## M 1.7e-13
##
## P value adjustment method: bonferroniTesi confermata. Il p-value presenta valori ben al di sotto del livello di significatività, il che ci permette di rifiutare l’ipotesi nulla e di affermare che vi sia una significativa diversità dei caratteri lunghezza e cranio in relazione al sesso dei neonati.
1.2) CREAZIONE DEL MODELLO DI REGRESSIONE
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(dati, upper.panel = panel.smooth, lower.panel = panel.cor)
Da questa analisi grafica è possibile notare come la variabile peso non
abbia una relazione lineare con tutte le altre variabili, anzi con la
variabile lunghezza la relazione sembra essere per lo più esponenziale,
mentre con la variabile gestazione e cranio invece sembra vi sia una
relazione di tipo logaritmica.
Come primo passo, procediamo con l’ipotesi di dipendenza lineare del modello, poi proveremo a definire anche un modello con relazioni non lineari e li confronteremo tra loro.
mod1 = lm(Peso~., data = dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1124.40 -181.66 -14.42 160.91 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6738.4762 141.3087 -47.686 < 2e-16 ***
## Anni.madre 0.8921 1.1323 0.788 0.4308
## N.gravidanze 11.2665 4.6608 2.417 0.0157 *
## Fumatrici -30.1631 27.5386 -1.095 0.2735
## Gestazione 32.5696 3.8187 8.529 < 2e-16 ***
## Lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## Cranio 10.4707 0.4260 24.578 < 2e-16 ***
## Tipo.partoNat 29.5254 12.0844 2.443 0.0146 *
## Ospedaleosp2 -11.2095 13.4379 -0.834 0.4043
## Ospedaleosp3 28.0958 13.4957 2.082 0.0375 *
## SessoM 77.5409 11.1776 6.937 5.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.2 on 10 and 2489 DF, p-value: < 2.2e-16Dalla creazione di un primo modello è possibile affermare quanto detto in precedenza, ovvero che le variabili antropometriche sono quelle che incidono in modo maggiormente significativo sul peso del neonato; allo stesso modo varaibili come Sesso e Gestazione assumonon anche esse un valore molto significa- tivo; al contrario, invece, le variabili come Fumatrici e Anni.madre sembrano non influire sul peso del bambino. Infine possiamo escludere la variabile Ospedale in quanto ha solo funzione di variabile di controllo. Il valore dell’R^2 adjusted riporta un valore del 73% circa come bontà del modello.
Passiamo quindi alla creazione di un secondo modello, eliminando, oltre alla variabile Ospedale, le variabili che il modello non considera per nulla significative.
mod2 = update(mod1, ~. - Ospedale - Anni.madre - Fumatrici )
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.31 -181.70 -16.31 161.07 2638.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.2971 135.9911 -49.322 < 2e-16 ***
## N.gravidanze 12.7558 4.3366 2.941 0.0033 **
## Gestazione 32.2713 3.7941 8.506 < 2e-16 ***
## Lunghezza 10.2864 0.3007 34.207 < 2e-16 ***
## Cranio 10.5057 0.4260 24.659 < 2e-16 ***
## Tipo.partoNat 30.0342 12.0969 2.483 0.0131 *
## SessoM 77.9285 11.1905 6.964 4.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1110 on 6 and 2493 DF, p-value: < 2.2e-16in questo secondo modello, lo sfoltimento di alcune variabili considerati non significative ha permesso alla variabile N.gravidanze di aumentare la propria significatività; la bontà del modello rimane tuttavia la medesima seppur con un valore molto elevato. Siccome in questo secondo modello è il tipo di parto ad assumere la minor significatività, proviamo ad eliminarla e a vedere il risultato che ne consegue
mod3 = update(mod2, ~.- Tipo.parto)
summary(mod3)##
## 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-16sostanzialmente il risultato rimane invariato
1.3) SELEZIONE DEL MIGLIOR MODELLO
Dopo una prima disamina dei diversi modelli, eseguiamo sia un test ANOVA che un BIc test per selezionare il miglior modello
anova(mod3, mod2)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2494 188065546
## 2 2493 187601677 1 463870 6.1643 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod2, mod3)## df BIC
## mod2 8 35221.75
## mod3 7 35220.10il test anova afferma che l’aggiuntaa della variabile tipo.parto non apporta una maggiore varianza spiegata in quanto il p-value dell’analisi risulta essere pocco significativo; inoltre la tecnica della minimizzazione di informazione di Bayes indica il mod3 come quello con il valore più basso. Pertanto saremo portati a scegliere il mod3. Infine per evitare problemi di multicollinearità, eseguiamo anche un test VIF sil modello scelto
car::vif(mod3)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023475 1.669189 2.074689 1.624465 1.040054i risultati sono tutti inferiori a 5 e pertanto possiamo affermare che non vi sono problemi di multicollinearità.
Creazione di un modello con relazioni non lineari Proviamo ora a creare un modello con relazioni non lineari, basandoci sulla precedente matrice di correlazione
plot(Peso,Lunghezza)
plot(Peso,Cranio)
plot(Peso,Gestazione)
plot(Peso,N.gravidanze)
plot(Peso, Sesso)
mod4 = update(mod3, ~. + I(log(Lunghezza)) + I(log(Gestazione)))
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(log(Lunghezza)) + I(log(Gestazione)), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1190.78 -181.94 -11.91 164.71 1319.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.102e+04 6.626e+03 9.209 < 2e-16 ***
## N.gravidanze 1.454e+01 4.239e+00 3.429 0.000616 ***
## Gestazione -2.185e+02 5.490e+01 -3.980 7.08e-05 ***
## Lunghezza 4.795e+01 3.461e+00 13.857 < 2e-16 ***
## Cranio 1.043e+01 4.185e-01 24.930 < 2e-16 ***
## SessoM 7.199e+01 1.097e+01 6.559 6.55e-11 ***
## I(log(Lunghezza)) -1.815e+04 1.665e+03 -10.901 < 2e-16 ***
## I(log(Gestazione)) 9.846e+03 2.067e+03 4.763 2.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268 on 2492 degrees of freedom
## Multiple R-squared: 0.7402, Adjusted R-squared: 0.7395
## F-statistic: 1014 on 7 and 2492 DF, p-value: < 2.2e-16anova(mod4, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(log(Lunghezza)) + I(log(Gestazione))
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 178987036
## 2 2494 188065546 -2 -9078510 63.199 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod3, mod4)## df BIC
## mod3 7 35220.10
## mod4 9 35112.05car::vif(mod4)## N.gravidanze Gestazione Lunghezza Cranio
## 1.025391 366.196712 288.643492 1.643714
## Sesso I(log(Lunghezza)) I(log(Gestazione))
## 1.048093 303.741988 385.809972Il modello ha previsto l’inserimento di due relazioni logaritmiche sulle variabili lunghezza e gestazione. Secondo il summary eseguito, risulterebbe un buon miglioramento del modello rispetto a mod3, grazie all’acquisizione di una maggiore significatività di tutte le variaibli indipendenti e con un R^2 adjusted maggiorato di un punto. Inoltre il test anova definisce una maggior varianza spiegata del modello non lineare e il BIC test riporta un valore inferiore a quello del mod3. il problema s pone al momento della verifica della multicollinearità dove i valori rialsciati dai fattori di inflazione della varianza asumono tutti valori ampiamente superiori a 5 delineando pertanto la presenza di multicollinearità tra i regressori.
Pettanto possiamo selezionare il mod3 quale modello di riferimento.
1.3) ANALISI DELLA QUALITÀ DEL MODELLO
Considerando che l’indice R e R^2 sono stati già analizzati, passeremo direttamente all’analisi dei residui.
Verificheremo i seguenti 4 assunti: 1) Distribuzione normale dei residui 2) Linearità 3) Omoschedasticità 4) Assenza di outliers
Iniziamo con una visualizzazione grafica
par(mfrow = c(2,2))
plot(mod3)
#verifica di normalità della distribuzione dei residui
shapiro.test(residuals(mod3))##
## Shapiro-Wilk normality test
##
## data: residuals(mod3)
## W = 0.97408, p-value < 2.2e-16plot(density(residuals(mod3)))
#omoschedasticità mediante il test di BREUSCH-PAGAN
lmtest::bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 90.253, df = 5, p-value < 2.2e-16#verifica di autocorrelazione mediante test di DURBIN_WATSON
lmtest::dwtest(mod3)##
## Durbin-Watson test
##
## data: mod3
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0dai primi test effettuati è possibile definire che i residui non si distribuiscono normalmente in quanto lo shapiro test ha rifiutato l’ipotesi nulla e non presenta omoschedasticità in quanto anche qui l’ipotesi nulla dettata dal test di Breusch Pagan è stata rifutata; invece i residui non presentano autocorrelazione poichè il test di Durbin Watson non ha rifiutato l’ipotesi nulla.
Verifichiamo pertanto la presenza di leverages e outliers
#leverages
lev = hatvalues(mod3)
plot(lev)
p=sum(lev)
soglia = 2*p/n
abline(h = soglia, col = "darkblue")
punti.lev = lev[lev>soglia]
punti.lev## 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.005775173length(punti.lev)## [1] 152l’analisi rileva la presenza di 152 punti di leverages.
Passiamo ora all’analisi degli outliers e poi trarremo le nostre conclusioni a mezzo della distanza di cook
#outliers
plot(rstudent(mod3))
abline(h = c(-2,2), col = "green")
car::outlierTest(mod3)## 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-03anche se l’analisi grafica riporta la presenza di molti outliers, il relativo test riporta invece solo la presenza di 3 outliers
#distanza di cook
cook = cooks.distance(mod3)
plot(cook)
max(cook)## [1] 0.8300965Secondo l’analisi della distanza di cook esiste un unico outliers prossimo al valore di 1, ovvero 1551, anche se non supera tale soglia.
Pertanto, prendiamo atto che la parte erratica non è pulita e teniamo come buono il nostro modello.
2) PREVISIONI E RISULTATI
mod3$coefficients## (Intercept) N.gravidanze Gestazione Lunghezza Cranio SessoM
## -6681.14450 12.47497 32.33206 10.24857 10.54020 77.99271eseguiamo alcuni test predittivi con il nosro modello.
pred1 = predict(mod3, newdata = data.frame(N.gravidanze = 2,
Gestazione = 39,
Sesso = "F",
Lunghezza = median(Lunghezza),
Cranio = median(Cranio)))
pred1## 1
## 3312.706pred2_m = predict(mod3, newdata = data.frame(N.gravidanze = 0,
Gestazione = 30,
Sesso = "M",
Lunghezza = median(Lunghezza[Lunghezza<400]),
Cranio = median(Cranio[Cranio<300])))
pred2_m## 1
## 1110.035pred2_f = predict(mod3, newdata = data.frame(N.gravidanze = 0,
Gestazione = 30,
Sesso = "F",
Lunghezza = median(Lunghezza[Lunghezza<400]),
Cranio = median(Cranio[Cranio<300])))
pred2_f## 1
## 1032.042pred3_f = predict(mod3, newdata = data.frame(N.gravidanze = 1,
Gestazione = 40,
Sesso = "F",
Lunghezza = median(Lunghezza),
Cranio = median(Cranio)))
pred3_f## 1
## 3332.564pred3_m = predict(mod3, newdata = data.frame(N.gravidanze = 1,
Gestazione = 40,
Sesso = "M",
Lunghezza = median(Lunghezza),
Cranio = median(Cranio)))
pred3_m## 1
## 3410.5563) VISUALIZZAZIONI
ggplot(data=dati)+
geom_point(aes( x=Gestazione,
y = Peso, col = Sesso),
position = "jitter") +
geom_smooth(aes(x=Gestazione, y = Peso, col = Sesso), method = "lm", se = F)+
geom_smooth(aes(x=Gestazione, y = Peso), col = "black",se = F)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggplot(data = dati) +
geom_point(aes(x= Gestazione, y = Peso, col = Fumatrici),
position = "jitter") +
geom_smooth(aes( x= Gestazione, y = Peso, col = Fumatrici), method = "lm")## `geom_smooth()` using formula = 'y ~ x'## Warning: The following aesthetics were dropped during statistical transformation:
## colour.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
ggplot(data = dati) +
geom_point(aes( x = Lunghezza, y = Peso, col = Sesso),
position = "jitter") +
geom_smooth(aes( x= Lunghezza, y = Peso, col = Sesso), method = "lm", se = F)+
geom_smooth(aes(x=Lunghezza, y = Peso), col = "black",se = F)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggplot(data = dati) +
geom_point(aes( x = Gestazione, y = Peso, col = N.gravidanze),
position = "jitter") +
geom_smooth(aes( x = Gestazione, y = Peso, col = N.gravidanze), method = "lm",
se = F) +
geom_smooth(aes(x=Gestazione, y = Peso), col = "black",se = F)## `geom_smooth()` using formula = 'y ~ x'## Warning: The following aesthetics were dropped during statistical transformation:
## colour.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
CONCLUSIONI Secondo i modelli anova e BIC abbiamo potuto scegliere il mod3 come miglior modello di regressione multipla del dataset Neonatal Health Solution, nella quale abbiamo impostato la variabile Peso come variabile dipendente e scelto come regressori signiicativi le variabili Gestazione, N.gravidanze, Sesso, Lunghezze e Cranio. Il modello proposto presenta una diretta propozionalità positiva tra regressori e variabile risposta. Il dataset presenta un’alta variabilità di casistiche che ne determina una disribuzione non di tipo normale e pertanto il modello fatica ad interpolare i dati sopratutto nei casi dove le frequenze sono minori, mentre migliora quando le casistiche presentao una maggioe frequenza. Il nostro modello descrive che il peso del neonato tende ad aumentare a mano a mano che la gestazione arriva al suo naturale compimento ovvero intorno alla 39esima settimana; come diretta conseguenza anche la lunghezza e la circonferenza del cranio tendono ad aumentare con l’aumentare del peso del neonato. In questa proporzionalità, esiste una certa differenza tra il sesso dei neonati, ovvero i maschi tendono ad avere un maggiore peso rispetto ai neonati di genere femminile a parità di settimane di gestazione. Infine, anche il numero di gravidanze pregresse sembra avere un effetto di normalità sulla madre, ovvero tende a stabilizzare il numero di settimane di gestazione verso il periodo target previsto e di conseguenza anche il peso e l’altezza del neonato tenderanno ad aumentare.