Previsione del Peso di Neonati
Avendo ricevuto la richiesta dalla Neonatal Health Solutions di un modello per poter prevedere il peso dei neonati e ricevuti i dati, ho deciso di iniziare con lo studio dei dati a mio dire più importanti. Questo, includendo la visualizzazione delle combinazioni dei dati per capire l’influenza degli stessi sul peso dei neonati, oltre che comprendere la situazione per le madri fumatrici e notare le differenze di gestazione.
Preparazione ed esplorazione dei dati
Variabili del dataset:
Età madre : quantitativa continua
N. gravidanze : quantitativa discreta
Fumatrice : qualitativa nominale
Gestazione : quantitativa continua
Peso : quantitativa continua
Diametro/lunghezza del cranio : quantitativa continua
Tipo di parto : qualitativa nominale
Ospedale di nascita : qualitativa nominale
Sesso del neonato : qualitativa nominale
Inzio con il preparare delle funzioni utili per il calcolo degli indici più importanti per le varie variabili.
library(ggplot2)
library(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, unionlibrary(moments)
library(car)## Caricamento del pacchetto richiesto: carData##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recodelibrary(lmtest)## Caricamento del pacchetto richiesto: zoo##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericlibrary(rgl)
library(knitr)
library(MASS)##
## Caricamento pacchetto: 'MASS'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## selectlibrary(Metrics)
df = read.csv("neonati.csv", sep = ",")
n = nrow(df)
#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.normalizzato = gini/((j-1)/j)
return(gini.normalizzato)
}
#curtosi
kurtosi.index = function(x){
mu = mean(x)
sigma = sd(x)
n = length(x)
m4 = sum((x-mu)**4)/n
kurtosi.index = m4/sigma**4 - 3
return(kurtosi.index)
}
quant = function(x) {
sd_val = sd(x)
var_val = var(x)
kurt_val = kurtosi.index(x)
return(c(sd = sd_val, var = var_val, kurt = kurt_val))
}In seguito sono andato a calcolare il summary dell’intero dataset, per avere un idea generale dei dati presenti, poi, sono andato a calcolare indici extra che possono essere util, pulisco i dati dalle variabili Na, e le età più estreme delle madri.
df = df[!df$Anni.madre <= 10, ]
df = na.omit(df)
summary(df)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. :13.00 Min. : 0.0000 Min. :0.00000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.00000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.00000 Median :39.00
## Mean :28.19 Mean : 0.9816 Mean :0.04163 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.00000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.00000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2498
## 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:2498 Length:2498
## Class :character Class :character
## Mode :character Mode :character
##
##
## peso = quant(df$Peso)
lungh = quant(df$Lunghezza)
cranio = quant(df$Cranio)
kable(peso)| x | |
|---|---|
| sd | 5.252294e+02 |
| var | 2.758659e+05 |
| kurt | 2.024728e+00 |
kable(lungh)| x | |
|---|---|
| sd | 26.328847 |
| var | 693.208160 |
| kurt | 6.473341 |
gini_gest = gini.index(df$Gestazione)
gest = c(quant(df$Gestazione), gini = gini_gest)
kable(gest)| x | |
|---|---|
| sd | 1.8689503 |
| var | 3.4929751 |
| kurt | 8.2465059 |
| gini | 0.8475314 |
kable(cranio)| x | |
|---|---|
| sd | 16.429469 |
| var | 269.927460 |
| kurt | 2.940112 |
Vedendo questi risultati, ho consultato un grafico per avere una visione chiara delle distribuzioni tra peso, sesso, numero di gravidanze passate della madre ed il tipo di parto necessario per il neonato. In più, per la gestazione possiamo notare un gini index quasi uniforme.
ggplot(data.frame(x = rnorm(2498, mean(df$Peso), peso[1]), y = df$Sesso), aes(x = x, fill = y))+
geom_density(alpha = .5, position = "identity")+
geom_vline(aes(xintercept = mean(df$Peso)), color = "red")+
theme_minimal()+
labs(
title = "Distribuzione normale del peso",
x = "Peso",
y = "Densità",
fill = "Sesso"
)
ggplot(df)+
geom_bar(aes(x = N.gravidanze, fill = Tipo.parto))+
theme_minimal()+
labs(
title = "Tipo di parto per numero di gravidanze già fatte",
x = "Numero Gravidanze",
y = "Conto",
fil = "Tipo Parto"
)
Possiamo notare che tra maschi e femmine c’è una buona differenza nel peso, graficamente, la distribuzione maschile si può notare essere, infatti, più alta e sottile ai latim, invece quella femminile è più piccola e stretta, entrambe presentano delle lievi distorsioni lunga la distribuzione.
Capiamo quindi che i maschi, hanno una maggioranza di peso di media 3300 grammi, con il minimo che si aggira sulla media dei 1000 grammi, ed un massimo attorno alla media dei 5500 grammi.
Le femmine, invece, notiamo che la loro media di picco si aggira alla media di 3600 grammi, con il minimo al di sotto della media dei 1000 grammi ed il massimo, come per i maschi, attorno alla media dei 5500 grammi.
Il secondo grafico rappresenta il rapporto tra il tipo di parto ed il numero di gravidanze che la madre ha già avuto. Notiamo che il parto cesareo all’inizio con le prime gravidanze arriva ad essere molto più alto rispetto al parto cesareo, ma con un numero di gravidanze più alto arriva a scomparire del tutto.
In alcuni ospedali si fanno più parti cesarei
Pensando che in alcuni ospedali ci sia una maggioranza di parti cesarei possiamo calcolare con il test di Chisq il tipo di parto in relazione all’ospedale, ed in seguito, associare questo ad un grafico a barre che mostra la differenza in altezza delle due colonne che rappresentano il tipo di parto per ogni ospedale.
Il test ci risulta un P-Value di ben 0.5778, dimostra che effettivamente la nostra ipotesi è corretta. Osservando il grafico otteniamo anche la conferma visiva, notando come nel secondo ospedali ci sia effettivamente un numero più elevato di parti cesarei rispetto agli altri due.
csq = chisq.test(df$Tipo.parto, df$Ospedale)
csq##
## Pearson's Chi-squared test
##
## data: df$Tipo.parto and df$Ospedale
## X-squared = 1.083, df = 2, p-value = 0.5819ggplot(df)+
geom_bar(aes(x = Ospedale, fill = Tipo.parto))+
theme_minimal()+
labs(
title = "Quantità dei tipi di parto negli ospedali",
x = "Ospedali",
y = "Conto",
fill = "Tipo Parto"
)
La media del peso e lunghezza del nostro campione è significativamente uguale a quella della popolazione
Per capire se il nostro campione può essere un esempio corretto rispetto alla popolazione, possiamo calcolare con un t-test se la media del campione è significativamente uguale a quella della popolazione, e poi creiamo ancora un grafico a barre per mostrare la differenza in altezza delle due classi peso e lunghezza tra campione e popolazione. Ho recuperato le medie di peso e lunghezza da questo sito web che ci garantisce un peso medio di 3300 grammi ed una lunghezza media di 500 millimetri.
media_pop_peso = 3300
media_pop_lung = 500
media_peso = mean(df$Peso)
media_lung = mean(df$Lunghezza)
kable(cbind(media_peso, media_lung))| media_peso | media_lung |
|---|---|
| 3284.184 | 494.6958 |
t.test(df$Peso, mu = media_pop_peso)##
## One Sample t-test
##
## data: df$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184t.test(df$Lunghezza, mu = media_pop_lung)##
## One Sample t-test
##
## data: df$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958ggplot(df)+
geom_boxplot(aes(y = Peso))+
geom_hline(aes(yintercept = media_pop_peso), color = "red")+
theme_minimal()+
labs(
title = "T.test tra il peso medio della popolazione e del campione",
y = "Peso"
)
ggplot(df)+
geom_boxplot(aes(y = Lunghezza))+
geom_hline(aes(yintercept = media_pop_lung), color = "red")+
theme_minimal()+
labs(
title = "T.test tra la lunghezza media della popolazione e del campione",
y = "Lunghezza"
)
Possiamo vedere con questi t-test, che ci ritornano un P-value di 0.1296 e per la lunghezza <0.00000000000000022, che la differenza tra la media campione e media popolazione non è significante per poter rendere inaffidabile il campione. Nei grafici dedicati, infatti, possiamo avere un’ulteriore sicurezza su questi risultati, vedendo come la linea rossa, che rappresenta la media in grammi della popolazione, sia quasi perfettamente nella stessa posizione della media del nostro campione.
Le misure antropometriche sono significativamente diverse tra i due sessi
Ora dobbiamo controllare se le misure delle varie parti del corpo tra maschi e femmine siano significativamente diverse. Usando la misura del cranio possiamo calcolare con il T-test se tra maschi e femmine ci sono differenze.
t_test_cranio = t.test(df$Cranio ~ df$Sesso)
t_test_cranio##
## Welch Two Sample t-test
##
## data: df$Cranio by df$Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M
## 337.6231 342.4586ggplot(df)+
geom_bar(aes(Cranio, fill = Sesso))+
theme_minimal()+
scale_y_continuous(breaks = seq(10, 250, 20))+
labs(
title = "Confronto tra le misure del cranio maschili e femminili",
x = "Misurazioni Cranio",
y = "Conto",
fill = "Sesso"
)
Anche con questo test, che ci riporta un P-value veramente insignificante, del valore di 0,001, ci fa capire, che non ci sono differenze significanti tra i due sessi. Come con le precedenti due domande, possiamo avere la sicurezza visiva utilizzando un grafico a barre, che ci mostra, che nonostante siano registrate delle misure del cranio molto più piccole per le femmine rispetto che per i maschi, queste non sono di significatività elevata per far pensare che le misure antropometriche possano essere differenti tra i due sessi.
Distribuzione di Frequenze per la classe del peso
Considerando che il nostro modello dovrà prevedere la classe del peso, ho deciso di ristrutturare il dataset permettendo alle varie variabili di rientrare in delle categorie specifiche di peso.
Ho diviso il dataset in 5 classi di peso, queste varianti da sotto i 1000 grammi fino ad arrivare sotto ai 5000 grammi. Con questo, ho deciso di andare a controllare la variabile della gestazione, per avere la conferma sulla sua significatività di interazione sulla nostra classe peso, e viceversa, con un grafico diviso per sesso.
df$peso_cut = cut(df$Peso, breaks = 5)
ggplot(df)+
geom_boxplot(aes(x = peso_cut, y = Gestazione, group = Sesso, fill = Sesso))+
theme_minimal()+
scale_x_discrete(labels = c("< 1000", "< 2000", "< 3000", "< 4000", "< 5000"))+
labs(
title = "Gestazione per peso in base al sesso",
x = "Classi di peso (in grammi)",
y = "Gestazione"
)
ggplot(df)+
geom_point(aes(x = format(peso_cut, scientific = FALSE), y = Gestazione, fill = Sesso, color = Sesso), position = "jitter")+
scale_x_discrete(labels = c("< 1000", "< 2000", "< 3000", "< 4000", "< 5000"))+
labs(
title = "Gestazione per peso in base al sesso",
x = "Classi di peso (in grammi)",
y = "Gestazione"
)+
theme_minimal()
Notiamo che nonostante si possa pensare il contrario, un peso più elevato non è strettamente legato ad una gestazione più lunga, al contrario, si registrano gestazioni più grandi per le classi di peso tra i <2000 grammi e i <3000 grammi, con una grande differenza fra la gestazione di quelle fasce di peso tra maschi e femmine, con le femmine che richiedono una gestazione molto più alta e sono molto più dense nel grafico.
Distribuzione di Frequenze per la classe della lunghezza
Ora, faccio lo stesso processo anche per la lunghezza, anche qui in 5 classi, a partire dai 310mm - 361mm ad arrivare alla classe dei 514mm - 565mm. Per parità rispetto alla classe del peso, ho deciso di registrare la gestazione divisa per la lunghezza.
df$lung_cut = cut(df$Lunghezza, breaks = 5)
ggplot(df)+
geom_boxplot(aes(x = lung_cut, y = Gestazione, group = Sesso, fill = Sesso))+
theme_minimal()+
labs(
title = "Gestazione per lunghezza in base al sesso",
x = "Classi di Lunghezza",
y = "Gestazione"
)
Possiamo vedere che nonostante per la classe peso non serviva essere di una classe alta per avere una gestazione elevata, per la lunghezza invece la gestazione più alta è presente nella lunghezza compresa tra i 463mm - 514mm, e la classe più alta a seguire. Anche qui, notiamo che le femmine hanno una gestazione molto più alta rispetto al genere maschile.
Come piccolo extra, ho deciso di andare a controllare se le madri fumatrici avessero una gestazione più o meno elevata rispetto alle madri non fumatrici.
ggplot(df, aes(x = Fumatrici, y = Gestazione, fill = Sesso))+
stat_summary(fun = "mean", geom = "bar", position = "dodge")+
theme_minimal()+
scale_x_continuous(breaks = seq(0, 1, 1))+
labs(
title = "Gestazione per donne fumatrici in base al sesso",
x = "Fumatrici",
y = "Gestazione",
fill = "Sesso"
) 
Notiamo che questo grafico le madri che non fumano tendono ad aver una gestazione molto più alta rispetto alle madri che fumano. Con una superiorità per il sesso maschile nelle settimane totali di gestazioni in entrambi i casi.
Modello di Regressione
Prima di creare il modello voglio andare a calcolare la correlazione tra le varie variabili del dataset con la classe peso.
cor.test(df$Peso, df$Gestazione)##
## Pearson's product-moment correlation
##
## data: df$Peso and df$Gestazione
## t = 36.694, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5658780 0.6168568
## sample estimates:
## cor
## 0.5919592cor.test(df$Peso, df$N.gravidanze)##
## Pearson's product-moment correlation
##
## data: df$Peso and df$N.gravidanze
## t = 0.11377, df = 2496, p-value = 0.9094
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.03694459 0.04149182
## sample estimates:
## cor
## 0.002277118cor.test(df$Peso, df$Lunghezza)##
## Pearson's product-moment correlation
##
## data: df$Peso and df$Lunghezza
## t = 65.71, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7812121 0.8099730
## sample estimates:
## cor
## 0.7960415cor.test(df$Peso, df$Cranio)##
## Pearson's product-moment correlation
##
## data: df$Peso and df$Cranio
## t = 49.642, df = 2496, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6845483 0.7240475
## sample estimates:
## cor
## 0.7048438cor.test(df$Peso, df$Anni.madre)##
## Pearson's product-moment correlation
##
## data: df$Peso and df$Anni.madre
## t = -1.1885, df = 2496, p-value = 0.2348
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06294109 0.01545144
## sample estimates:
## cor
## -0.02378138Possiamo notare che:
Per la gestazione abbiamo una correlazione media
Per la lunghezza e la misura del cranio abbiamo una correlazione grande
Per il numero di gravidanze, la variabile delle fumatrici e gli anni della madre abbiamo una correlazione nulla
In seguito vado a rappresentare queste correlazioni usando la funzione pairs().
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
num_data = cbind(Anni=df$Anni.madre, N.Grav = df$N.gravidanze, Fumatrici = df$Fumatrici, Gestazione = df$Gestazione, Peso=df$Peso,Lungh = df$Lunghezza, Cranio = df$Cranio)
pairs(num_data,upper.panel = panel.smooth, lower.panel = panel.cor)
Una volta confermate, andiamo finalmente a creare il nostro modello.
Creazione del modello
mod1 = lm(Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + Cranio + N.gravidanze + Sesso, data = df)
best = stepAIC(mod1, direction = "both", trace= TRUE)## Start: AIC=28064.05
## Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + Cranio +
## N.gravidanze + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 44292 187949505 28063
## - Fumatrici 1 91474 187996688 28063
## <none> 187905214 28064
## - N.gravidanze 1 446756 188351970 28068
## - Sesso 1 3658879 191564093 28110
## - Gestazione 1 5587942 193493156 28135
## - Cranio 1 45789523 233694736 28607
## - Lunghezza 1 87128339 275033553 29014
##
## Step: AIC=28062.64
## Peso ~ Gestazione + Lunghezza + Fumatrici + Cranio + N.gravidanze +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 92548 188042054 28062
## <none> 187949505 28063
## + Anni.madre 1 44292 187905214 28064
## - N.gravidanze 1 643981 188593487 28069
## - Sesso 1 3666800 191616305 28109
## - Gestazione 1 5544825 193494331 28133
## - Cranio 1 46056754 234006260 28608
## - Lunghezza 1 87116561 275066067 29012
##
## Step: AIC=28061.87
## Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188042054 28062
## + Fumatrici 1 92548 187949505 28063
## + Anni.madre 1 45366 187996688 28063
## - N.gravidanze 1 621053 188663107 28068
## - Sesso 1 3650790 191692844 28108
## - Gestazione 1 5477493 193519547 28132
## - Cranio 1 46098547 234140601 28608
## - Lunghezza 1 87532691 275574744 29015summary(best)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze +
## Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16anova(best, mod1)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + Cranio +
## N.gravidanze + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2490 187905214 2 136840 0.9067 0.404BIC(best, mod1)## df BIC
## best 7 35193.65
## mod1 9 35207.48Utilizzando la funzione stepAIC con entrambe le direzioni impostate, abbiamo automaticamente testato diversi tipi di modello andando a rimuovere e riaggiungere le varie variabili finchè non abbiamo trovato il modello più ottimizzato. Per andare a verificarlo manualmente, possiamo controllare usando il BIC e l’ANOVA, vedendo quale delle due tra il modello che ci ha dato la funzione, ed il modello originario ha lo score più basso. In più, con la funzione stepAIC possiamo già vedere nelle varie prove che ha fatto quale score ha ottenuto ogni modello.
Guardando invece i coefficienti del modello, notiamo che per ogni settimana di gestazione in più il peso aumenta in media di 32.38 grammi , per i maschi, 77.98 grammi in più delle femmine, la lunghezza di 10.24mm, ed il diametro del cranio di 10.54, in più, a seconda del n. di gravidanze il peso medio sarà di 12.45 grammi in più.
Valutazione del modello
Ora che abbiamo il modello possiamo fare le valutazioni, partiamo con il fare un plot per i 4 grafici principali:
La linearità tra residui e valori predetti
Controllare se i dati residui sono normali
Controllare l’omoschedasticità dei dati
Individuare dati estremi influenti
par(mfrow = c(2, 2))
plot(best)
Per il primo grafico, vediamo che la distribuzione dei residui rispetto alle predizioni del modello non è una dispersione casuale, ma è più una macchia di dati con una “coda”, cosa che ci fa capire esserci non-linearità tra i dati.
Nel secondo grafico, per osservare la normalità dei dati residui, ci mostra una linea retta di dati curva nelle code. Nel caso di questo dataset penso sia legato alla distribuzione dei dati non perfettamente gaussiana.
Nel terzo grafico, vediamo che i punti del grafico seguono più o meno, la stessa distribuzione del primo grafico, in questo caso però sono segno di omoschedasticità tra i dati.
Nell’ultimo notiamo che maggior parte dei dati rimane nel range della distanza di cook’s con un paio di outliers che vanno oltre la distanza minima, mentre un singolo dato rimane nella soglia di accettazione, in alto a destra, nonostante questo però, non mi sento di bocciare il modello.
lev = hatvalues(best)
par(mfrow = c(1, 1))
plot(lev)
p = sum(lev)
soglia = 2* p/n
abline(h = soglia, col = 2)
In seguito abbiamo il grafico dedicato all’influenza delle osservazioni sulla stima del modello, possiamo notare come non ci sono troppi dati estremi che ci potrebbero far pensare ad una brutta uscita del modello.
lev[lev>soglia]## 13 15 34 67 89 96
## 0.005632888 0.007056483 0.006748974 0.005896189 0.012817563 0.005351869
## 101 106 131 134 151 155
## 0.007528094 0.014487904 0.007237755 0.007553514 0.010889886 0.007209736
## 161 189 190 204 205 206
## 0.020341133 0.004894598 0.005366831 0.014494112 0.005351828 0.009482503
## 220 294 305 310 312 315
## 0.007403130 0.005914307 0.005445189 0.028815300 0.013173723 0.005386836
## 378 440 442 445 486 492
## 0.015942080 0.005405733 0.007725708 0.007511227 0.005165714 0.008274392
## 497 516 582 587 592 614
## 0.005167522 0.013080707 0.011667393 0.008415872 0.006385013 0.005300091
## 638 656 657 684 697 702
## 0.006693153 0.005934494 0.005323517 0.008825948 0.005864920 0.005203163
## 729 748 750 757 765 805
## 0.005024448 0.008567254 0.006943755 0.008146487 0.006076412 0.014358658
## 828 893 895 913 928 946
## 0.007180179 0.005076266 0.005297643 0.005574016 0.022745493 0.006910965
## 947 949 956 985 1008 1014
## 0.008409518 0.004801229 0.007791530 0.007040035 0.005343994 0.008474086
## 1049 1067 1091 1106 1130 1166
## 0.004956561 0.008467035 0.008940030 0.005967511 0.031739177 0.005513581
## 1181 1188 1200 1219 1238 1248
## 0.005678043 0.006481533 0.005493009 0.030697778 0.005912496 0.014631359
## 1273 1291 1293 1311 1321 1325
## 0.007089553 0.006118589 0.006074423 0.009626391 0.009295311 0.004857340
## 1356 1357 1385 1395 1400 1402
## 0.005306940 0.006967713 0.012641828 0.005129491 0.005932524 0.004816704
## 1411 1420 1428 1429 1450 1505
## 0.008049539 0.005156670 0.008195421 0.021758961 0.015106684 0.013334256
## 1551 1553 1556 1573 1593 1606
## 0.048802841 0.008507417 0.005923458 0.005049368 0.005624961 0.005009312
## 1610 1617 1619 1628 1686 1693
## 0.008725821 0.004869997 0.015069038 0.005070717 0.009356578 0.005079185
## 1701 1712 1718 1727 1735 1780
## 0.010846383 0.006993461 0.006961226 0.013303625 0.004886226 0.025544997
## 1781 1809 1827 1868 1892 1962
## 0.016833696 0.008711082 0.006068655 0.005206137 0.005333985 0.005541287
## 1967 1977 2037 2040 2046 2086
## 0.005339716 0.006928794 0.004890003 0.011504028 0.005471894 0.013194769
## 2089 2098 2114 2115 2120 2140
## 0.006293791 0.005097113 0.013318873 0.011779730 0.018667270 0.006244737
## 2146 2148 2149 2157 2175 2200
## 0.005804705 0.007930323 0.013589469 0.005910248 0.032531736 0.011679031
## 2215 2216 2220 2221 2224 2225
## 0.004894044 0.008120291 0.005415586 0.021633907 0.005841119 0.005593576
## 2244 2257 2307 2317 2318 2337
## 0.006929530 0.006171101 0.013972901 0.007677230 0.004834368 0.005230839
## 2359 2408 2422 2436 2437 2452
## 0.010068259 0.009702475 0.021536269 0.004986609 0.023951342 0.023848222
## 2458 2471 2478
## 0.008509142 0.020905930 0.005777111plot(rstudent(best))
abline(h = c(-2, 2), col = 2)
outlierTest(best)## rstudent unadjusted p-value Bonferroni p
## 1551 10.046230 2.6345e-23 6.5810e-20
## 155 5.025345 5.3818e-07 1.3444e-03
## 1306 4.824963 1.4848e-06 3.7092e-03Con questo invece calcoliamo se i valori di leverage rimangono all’interno del nostro range dei dati residui dalla funzione rstudent(), questo ci aiuta ad identificare osservazioni “strane” andando a visualizzare la varianza del modello, nel nostro caso possiamo accettare il modello, non avendo troppe osservazioni estreme all’infuori della nostra soglia.
Per avere una sicurezza finale, andiamo a fare il bptest, che ci presenta un p-value basso, quindi notiamo che il modello ha un pò di problemi di omoschedasticità, il dwtest, che ci presenta dei valori che ci fanno notare come il modello non ha problemi di autocorrelazione tra i valori residui, il test di shapiro, che ci presenta un p-value basso, che spiega che i nostri residui non seguono una distribuzione normale, ed infine, calcoliamo l’rmse, che ci presenta un valore basso rispetto alla scala del peso dei nostri dati.
bptest(best)##
## studentized Breusch-Pagan test
##
## data: best
## BP = 90.297, df = 5, p-value < 2.2e-16dwtest(best)##
## Durbin-Watson test
##
## data: best
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(best))##
## Shapiro-Wilk normality test
##
## data: residuals(best)
## W = 0.97414, p-value < 2.2e-16summary(best)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze +
## Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16kable(rmse(df$Peso, predict(best, df)))| x |
|---|
| 274.3666 |
Nonostante il piccolo problema di omoschedasticità e i residui che non seguono perfettamente una scala normale, posso sentirmi comunque sicuro sul fatto che il modello non ha problemi nella sua previsione del peso. Infatti, se ora proviamo a fare una previsione considerando una madre con gestazione di 39 settimane 2 gravidanze passate, con il neonato di lunghezza media ed il cranio di misura media e di sesso femminile, riceviamo come previsione dal modello un peso di media 3258.615 grammi.
media_cranio = mean(df$Cranio)
new_data = data.frame(Gestazione = 39, Lunghezza = media_lung, Cranio = media_cranio, N.gravidanze = 2, Sesso = "F")
results = predict.lm(best, newdata = new_data)
kable(results)| x |
|---|
| 3258.727 |
Grafici extra per l’interazione tra variabili
ggplot(data = df,
aes(x = factor("N. gravidanze"), y = Peso, fill = factor("N. gravidanze")))+
geom_violin()+
geom_boxplot(width = 0.1, outlier.shape = NA, fill = "white", alpha = 0.6)+
labs(
title = "Distribuzione del peso in base al numero di gravidanze",
x = "N. Gravidanze",
y = "Peso in grammi",
fill = "N. Gravidanze"
)+
theme_minimal()
ggplot(data = df)+
geom_boxplot(aes(x = Fumatrici, y = Peso, fill = Sesso))+
labs(
title = "Peso in relazione alla madre fumatrice",
x = "Madre Fumatrice",
y = "Peso in grammi",
fill = "Sesso"
)+
theme_minimal()
ggplot(data = df)+
geom_point(aes(x = Lunghezza,
y = Peso,
col = Sesso), position = "jitter")+
geom_smooth(aes(x = Lunghezza,
y = Peso,
col = Sesso), se = F, method = "lm")+
geom_smooth(aes(x = Lunghezza,
y = Peso), col = "black", se = F, method = "lm")+
theme_minimal()+
labs(
title = "Influenza della Lunghezza sul Peso",
x = "Lunghezza",
y = "Peso",
fill = "Sesso"
)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
ggplot(data = df)+
geom_point(aes(x = Gestazione,
y = Peso,
col = Sesso), position = "jitter")+
geom_smooth(aes(x = Gestazione,
y = Peso,
col = Sesso), se = F, method = "lm")+
geom_smooth(aes(x = Gestazione,
y = Peso), col = "black", se = F, method = "lm")+
theme_light()## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Con questi ultimi grafici abbiamo infine un ultima visuale sulle variabili che secondo me possono avere un influenza maggiore sul peso del neonato. Notiamo come la gestazione attorno alle 40 settimane ha un’influenza elevata sul peso ma che comunque nelle “prime” settimane di gestazione questa non influenza troppo il peso, ma nelle settimane più avanti invece il peso riceve molta più influenza, questo perchè il bambino è ora quasi o totalmente formato. Vediamo che le donne fumatrici non hanno una grande influenza sul peso mentre le donne che non fumano si. ed infine vediamo che (come dovremmo aspettarci normalmente), la lunghezza ed il peso sono molto influenti l’una sull’altra.