Previsione Peso 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 discreta
Peso : quantitativa continua
Diametro/lunghezza del cranio : quantitativa continua
Tipo di parto : qualitativa nominale
Ospedale di nascita : qualitativa nominale
Sesso del neonato : qualitativa nominale
Come primo step ho preso ogni singolo dato e cercato di capire quali variabili ed indici sarebbe opportuno calcolare per ognuno per capire meglio la loro interferenza nel dataset. Quindi, ho deciso di preparare inizialmente le funzioni per il calcolo veloce dell’indice di eterogeneità di Gini e la Curtosi.
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)
}In seguito, ho calcolato gli indici principali per le varie variabili utili.
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
##
##
## sd_peso = sd(df$Peso)
var_peso = var(df$Peso)
kurt_peso = kurtosi.index(df$Peso)
kable(cbind(kurt_peso, sd_peso, var_peso))| kurt_peso | sd_peso | var_peso |
|---|---|---|
| 2.024728 | 525.2294 | 275865.9 |
sd_lung = sd(df$Lunghezza)
var_lung = var(df$Lunghezza)
kurt_lung = kurtosi.index(df$Lunghezza)
kable(cbind(kurt_lung, sd_lung, var_lung))| kurt_lung | sd_lung | var_lung |
|---|---|---|
| 6.473341 | 26.32885 | 693.2082 |
sd_gest = sd(df$Gestazione)
var_gest = var(df$Gestazione)
gini_gest = gini.index(df$Gestazione)
kurt_gest = kurtosi.index(df$Gestazione)
kable(cbind(gini_gest, kurt_gest, sd_gest, var_gest))| gini_gest | kurt_gest | sd_gest | var_gest |
|---|---|---|---|
| 0.8475314 | 8.246506 | 1.86895 | 3.492975 |
sd_cranio = sd(df$Cranio)
var_cranio = var(df$Cranio)
kurt_cranio = kurtosi.index(df$Cranio)
kable(cbind(kurt_cranio, sd_cranio, var_cranio))| kurt_cranio | sd_cranio | var_cranio |
|---|---|---|
| 2.940112 | 16.42947 | 269.9275 |
Notiamo che tutte queste variabili hanno una una distribuzione Leptocurtica, designata dalla kurtosi superiore a 0, la variabile della gestazione ha un gini index quasi massimo, essendo molto vicino a 1, mentre la deviazione standard e la variazione portano ad indicare nel caso del peso, una grande dispersione e variabilità dei dati, per la lunghezza si nota una bella differenza tra le osservazioni, ma comunque nella norma, la gestazione suggerisce invece che i valori sono abbastanza concentrati nella media e senza molta varianza, infine per la variabile del cranio, notiamo che i dati non si dispergono più di tanto, ma che hanno comunque una grande varianza.
In seguito ho deciso di rappresentare i dati in dei grafici.
ggplot(data.frame(x = rnorm(2498, mean(df$Peso), sd_peso), 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"
)
Con questi due grafici, notiamo nel primo che il peso fra maschi e femmine ha una differenza abbastanza visibile, con la distribuzione del sesso maschile che si può vedere essere più alta e sottile ai lati e con le code leggermente più ristrette rispetto al sesso femminile. Questo ci fa capire che la buona maggioranza dei neonati maschi ha un peso che si aggira attorno ai 3300 grammi, un minimo di attorno ai 1000 grammi ed un massimo circa 5500 grammi, mentre per le femmine notiamo la grande maggioranza attorno ai 3600 grammi, un minimo al di sotto dei 1000 grammi ed un massimo attorno ai 5500 come per i maschietti.
Con il secondo grafico abbiamo invece una visuale più chiara sul rapporto di tipo di parto necessario rispetto al numero di gravidanze passate della madre, con la grande maggioranza, sia per parto naturale che cesareo sulla prima gravidanza, però notiamo che all’aumentare del numero di gravidanze diminuisce quasi fino allo zero la presenza dei parti cesarei, molto meno necessari.
In alcuni ospedali si fanno più parti cesarei
Abbiamo modo di pensare che a seconda degli ospedali si presenta la maggioranza di parti cesarei rispetto ai parti naturali, per capire questo, ho deciso di avere un grafico associato ad un test di Chisq, che ci mostra un P-Value di ben 0.5778, che ci dimostra che c’è effettivamente che la nostra ipotesi sia 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 garantire l’efficacia del nostro modello abbiamo anche bisogno di garantire che le nostre medie dei dati campione siano rispecchiate anche al di fuori dei 2500 campioni registrati in questo dataset, possiamo garantircene usando un test di T. Student e un grafico. Per paragonare le nostre medie alle medie della popolazione, sono andato a recuperare 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 <0.0….22, che significa che per il peso molto probabilmente, le differenze presenti non sono abbastanza grandi da poter dire che siano significativamente lontane tra loro, mentre per la lunghezza non ci sono differenze tra loro. 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
In quest’ultimo quesito che ci viene richiesto, dobbiamo controllare se le misure delle varie parti del corpo tra maschi e femmine siano significativamente diverse. Noi utilizzeremo la misura del cranio, quindi possiamo calcolare con il T-test se tra maschi e femmine ci sono differenze significanti di quest’ultima.
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))+
scale_x_continuous(breaks = seq(100, 400, 20))
labs(
title = "Confronto tra le misure del cranio maschili e femminili",
x = "Misurazioni Cranio",
y = "Conto",
fill = "Sesso"
)## $x
## [1] "Misurazioni Cranio"
##
## $y
## [1] "Conto"
##
## $fill
## [1] "Sesso"
##
## $title
## [1] "Confronto tra le misure del cranio maschili e femminili"
##
## attr(,"class")
## [1] "labels"Notiamo che come risultato ci da un valore P veramente molto basso, il che ci fa pensare che non ci sono delle differenze notevoli tra i due sessi.
Il grafico ci mostra che a seconda del sesso però c’è una quantità maggiore di osservazioni sulle varie misure del cranio.
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_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
Oltre che al peso, secondo me una variabile a cui si dovrebbe prestare attenzione in questo genere di modelli è la lunghezza del neonato, per cui, ho deciso di dividere anche questa metrica per classi, 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_point(aes(x = lung_cut, y = Gestazione, fill = Sesso, color = Sesso),
position = "jitter")+
labs(
title = "Gestazione per lunghezza in base al sesso",
x = "Classi di Lunghezza",
y = "Gestazione"
)+
theme_minimal()
Possiamo vedere che nonostante per la classe peso non serviva essere di una classe alta per avere una gestazione elevata, per la lunghezza invece, notiamo che 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)+
geom_bar(aes(x = Fumatrici, y = Gestazione, fill = Sesso), stat = "identity")+
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.
Modello di Regressione
Prima di effettivamente creare un modello di regressione lineare multipla, ho bisogno di andare a vedere che correlazione ed influenza hanno le varie variabili sul peso.
cor.test(df$Peso, df$Gestazione) # Correlazione media##
## 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) # Corr quasi nulla##
## 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) # Correlazione alta##
## 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) # Correlazione alta##
## 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) # Correlazione quasi nulla##
## 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.02378138Come commentato, possiamo notare che per la gestazione abbiamo una correlazione media, quindi non influisce ne troppo ma neanche poco. Al contrario la lunghezza e la misura del cranio che hanno invece una correlazione più grande, quindi, influiscono molto sul peso (come ci si dovrebbe aspettare). Infine abbiamo variabili come il numero di gravidanze, se la madre è fumatrice e gli anni della madre, che non influiscono praticamente di nulla. Con queste abbiamo più o meno idea di che variabili integrare nel nostro modello.
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, Gestazione = df$Gestazione, Peso=df$Peso,Lungh = df$Lunghezza, Cranio = df$Cranio)
pairs(num_data,upper.panel = panel.smooth, lower.panel = panel.cor)
mod1 = lm(Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + Cranio + N.gravidanze + Sesso, data = df)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici +
## Cranio + N.gravidanze + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.6 -181.3 -15.7 163.6 2630.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6712.2405 141.3339 -47.492 < 2e-16 ***
## Gestazione 32.9472 3.8288 8.605 < 2e-16 ***
## Lunghezza 10.2316 0.3011 33.979 < 2e-16 ***
## Anni.madre 0.8803 1.1491 0.766 0.444
## Fumatrici -30.3958 27.6080 -1.101 0.271
## Cranio 10.5198 0.4271 24.633 < 2e-16 ***
## N.gravidanze 11.3789 4.6767 2.433 0.015 *
## SessoM 78.0787 11.2132 6.963 4.24e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared: 0.7272, Adjusted R-squared: 0.7264
## F-statistic: 948.3 on 7 and 2490 DF, p-value: < 2.2e-16anova(mod1)## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Gestazione 1 241379330 241379330 3198.6049 < 2.2e-16 ***
## Lunghezza 1 206092456 206092456 2731.0057 < 2.2e-16 ***
## Anni.madre 1 1307180 1307180 17.3219 3.262e-05 ***
## Fumatrici 1 71365 71365 0.9457 0.330915
## Cranio 1 47900452 47900452 634.7462 < 2.2e-16 ***
## N.gravidanze 1 522265 522265 6.9207 0.008573 **
## Sesso 1 3658879 3658879 48.4851 4.242e-12 ***
## Residuals 2490 187905214 75464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod1)## [1] 35207.48cat("Dal mod1 togliamo la var Cranio\n")## Dal mod1 togliamo la var Craniomod2 = update(mod1, ~. - Cranio)
summary(mod2)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici +
## N.gravidanze + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1188.1 -192.3 -16.5 184.3 3597.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5400.8362 145.9774 -36.998 < 2e-16 ***
## Gestazione 47.2644 4.2196 11.201 < 2e-16 ***
## Lunghezza 13.5730 0.2997 45.282 < 2e-16 ***
## Anni.madre 2.3317 1.2796 1.822 0.0685 .
## Fumatrici -36.3380 30.7813 -1.181 0.2379
## N.gravidanze 20.5723 5.1978 3.958 7.77e-05 ***
## SessoM 87.9930 12.4944 7.043 2.43e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 306.3 on 2491 degrees of freedom
## Multiple R-squared: 0.6607, Adjusted R-squared: 0.6599
## F-statistic: 808.6 on 6 and 2491 DF, p-value: < 2.2e-16anova(mod2, mod1)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + 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 2491 233694736
## 2 2490 187905214 1 45789523 606.77 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod2, mod1) # Aumenta tutto quindi non ne vale la pena rimuoverlo## df BIC
## mod2 8 35744.41
## mod1 9 35207.48cat("Dal mod1 togliamo la var Anni.Madre\n")## Dal mod1 togliamo la var Anni.Madremod3 = update(mod1, ~. - Anni.madre)
summary(mod3)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Fumatrici + Cranio +
## N.gravidanze + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.24 -181.32 -15.73 162.95 2635.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6682.2637 135.7983 -49.207 < 2e-16 ***
## Gestazione 32.6437 3.8079 8.573 < 2e-16 ***
## Lunghezza 10.2309 0.3011 33.979 < 2e-16 ***
## Fumatrici -30.5728 27.6048 -1.108 0.26818
## Cranio 10.5366 0.4265 24.707 < 2e-16 ***
## N.gravidanze 12.6996 4.3470 2.921 0.00352 **
## SessoM 78.1596 11.2117 6.971 4.01e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265
## F-statistic: 1106 on 6 and 2491 DF, p-value: < 2.2e-16anova(mod3, mod2)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Fumatrici + Cranio + N.gravidanze +
## Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Anni.madre + Fumatrici + N.gravidanze +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 187949505
## 2 2491 233694736 0 -45745231BIC(mod3, mod1) # Diminuisce quindi va bene rimuoverlo## df BIC
## mod3 8 35200.24
## mod1 9 35207.48cat("Dal mod3 togliamo la var Fumatrici\n")## Dal mod3 togliamo la var Fumatricimod4 = update(mod3, ~. - Fumatrici)
summary(mod4)##
## 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(mod4, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Fumatrici + Cranio + N.gravidanze +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2491 187949505 1 92548 1.2266 0.2682BIC(mod4, mod3)## df BIC
## mod4 7 35193.65
## mod3 8 35200.24cat("Dal mod4 togliamo la var Lunghezza\n")## Dal mod4 togliamo la var Lunghezzamod5 = update(mod4, ~. - Lunghezza)
summary(mod5)##
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + N.gravidanze + Sesso,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1318.65 -218.40 -12.46 212.64 1537.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6229.0089 163.5785 -38.080 <2e-16 ***
## Gestazione 93.2296 4.0604 22.961 <2e-16 ***
## Cranio 17.1029 0.4605 37.140 <2e-16 ***
## N.gravidanze 5.0793 5.2482 0.968 0.333
## SessoM 117.9239 13.4947 8.739 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 332.5 on 2493 degrees of freedom
## Multiple R-squared: 0.5999, Adjusted R-squared: 0.5993
## F-statistic: 934.6 on 4 and 2493 DF, p-value: < 2.2e-16anova(mod5, mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Cranio + N.gravidanze + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 275574744
## 2 2492 188042054 1 87532691 1160 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod5, mod4) # Aumenta quindi non lo facciamo## df BIC
## mod5 6 36140.54
## mod4 7 35193.65cat("Dal mod4 togliamo la var Sesso\n")## Dal mod4 togliamo la var Sessomod6 = update(mod4, ~. - Sesso)
anova(mod6, mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 191692844
## 2 2492 188042054 1 3650790 48.382 4.466e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod6, mod4) # Anche qua aumenta rispetto al mod4## df BIC
## mod6 6 35233.86
## mod4 7 35193.65cat("Dal mod4 togliamo la var Fumatrici\n")## Dal mod4 togliamo la var Fumatricimod7 = update(mod4, ~. - N.gravidanze)
anova(mod7, mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 188663107
## 2 2492 188042054 1 621053 8.2304 0.004154 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod7, mod4) # Diminuisce di così poco che non ne vale la pena## df BIC
## mod7 6 35194.06
## mod4 7 35193.65AIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7) # Il 4° è il migliore## df AIC
## mod1 9 35155.07
## mod2 8 35697.83
## mod3 8 35153.66
## mod4 7 35152.89
## mod5 6 36105.60
## mod6 6 35198.92
## mod7 6 35159.12BIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7) # Idem qua## df BIC
## mod1 9 35207.48
## mod2 8 35744.41
## mod3 8 35200.24
## mod4 7 35193.65
## mod5 6 36140.54
## mod6 6 35233.86
## mod7 6 35194.06cat("Il modello con minor p-value, AIC e BIC è il modello 4\n")## Il modello con minor p-value, AIC e BIC è il modello 4In sostanza, in queste righe sono andato manualmente a testare diverse variazioni del modello base, in cui ho inserito tutte le variabili, piano piano sottraendo o ri-aggiungendo (come nel caso del modello 3 in cui ho ri-aggiunto la variabile cranio rimossa nel modello 2) variabili rimosse in precedenza.
Ad ogni modello che creo vado a calcolare il BIC e l’ANOVA e andiamo a rimuovere o aggiungere variabili ad un nuovo modello in base all’aumento o discesa di questi due valori.
Infine, calcoliamo queste ultime due con tutti i modelli creati per ricevere in output quale modello secondo le due formule è il migliore. Notiamo come il modello che viene ritenuto migliore (con i risultati di entrambe più basso) è il modello 4 con, ad esempio, il BIC risultante di 35193,65.
Andando a studiare meglio il modello 4 con summary notiamo un paio di cose.
summary(mod4)##
## 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-16Notiamo che per ogni settimana di gestazione in più, il peso aumenta di circa 32.38g, la lunghezza aumenta di 10.24mm ed il diametro del cranio di 10.24mm.
A seconda del numero di gravidanze da parte della madre, il peso medio aumenta di circa 12.45g, infine per i neonati di sesso maschile, il peso medio è di 77.98g in più rispetto alle femmine.
Valutazione del modello
Siamo riusciti a trovare il nostro modello, garantendoci sia il migliore possibile nella nostra situazione, ora, possiamo procedere con la sua valutazione, per farlo dobbiamo pensare a visualizzare e controllare se i dati estremi riportati nel modello non siano troppo elevati o troppo piccoli rispetto al resto dei dati. Per fare questo possiamo partire con il visualizzare quattro grafici che ci mostrano, in ordine per riga:
- 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(mod4)
Il primo possiamo notare come i dati rimangono tutti abbastanza compattati, seguendo una linea curva ed in fondo al grafico si annicchiano in una “macchia”, eccetto qualche dato estremo.
Questo ci fa capire che abbiamo un paio di dati estremi che si discostano dal resto dei dati, ma nulla di significanza importante.
Nel secondo grafico, possiamo notare che i dati residui sono pressochè normali poichè seguono una linea retta che si curva solo agli estremi.
Nel terzo possiamo notare che i dati sono abbastanza sparsi per tutto il grafico, suggerendo che i dati avrebbero bisogno di una trasformazione o modifica dei dati per migliorare ancora di più, secondo me però in questo caso non è necessario applicarne.
Infine, nell’ultimo grafico notiamo che la 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(mod4)
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(mod4))
abline(h = c(-2, 2), col = 2)
outlierTest(mod4)## 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.
cook = cooks.distance(mod4)
plot(cook)
Per la distanza di cook possiamo invece notare che quasi tutti i dati non hanno una grande influenza sul peso, tranne un paio di outliers che invece possono essere più o meno influenti tra loro sul peso.
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(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 90.297, df = 5, p-value < 2.2e-16dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(mod4))##
## Shapiro-Wilk normality test
##
## data: residuals(mod4)
## W = 0.97414, p-value < 2.2e-16summary(mod4)##
## 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(mod4, 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.
Per avere qualche sostegno extra per il 100% di sicurezza che sia il modello migliore che possiamo avere, proverei a testare su una variabile che mi sembra più adatta, se un modello con l’effetto quadratico oppure di interazione risultano migliori.
media_cranio = mean(df$Cranio)
interaz_gest_sesso = update(mod4, ~. + Gestazione:Sesso)
summary(interaz_gest_sesso)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze +
## Sesso + Gestazione:Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.70 -181.22 -15.19 163.39 2637.64
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6573.1588 167.0900 -39.339 < 2e-16 ***
## Gestazione 29.5202 4.5863 6.437 1.46e-10 ***
## Lunghezza 10.2494 0.3008 34.071 < 2e-16 ***
## Cranio 10.5422 0.4265 24.721 < 2e-16 ***
## N.gravidanze 12.4098 4.3416 2.858 0.00429 **
## SessoM -183.6605 234.8934 -0.782 0.43435
## Gestazione:SessoM 6.7047 6.0124 1.115 0.26490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7272, Adjusted R-squared: 0.7265
## F-statistic: 1106 on 6 and 2491 DF, p-value: < 2.2e-16anova(interaz_gest_sesso)## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Gestazione 1 241379330 241379330 3199.1571 < 2.2e-16 ***
## Lunghezza 1 206092456 206092456 2731.4773 < 2.2e-16 ***
## Cranio 1 48941070 48941070 648.6478 < 2.2e-16 ***
## N.gravidanze 1 731441 731441 9.6943 0.001869 **
## Sesso 1 3650790 3650790 48.3863 4.456e-12 ***
## Gestazione:Sesso 1 93827 93827 1.2435 0.264897
## Residuals 2491 187948227 75451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(interaz_gest_sesso, mod4)## df AIC
## interaz_gest_sesso 8 35153.64
## mod4 7 35152.89BIC(interaz_gest_sesso, mod4)## df BIC
## interaz_gest_sesso 8 35200.22
## mod4 7 35193.65quad_gestazione = update(mod4, ~. + I(Gestazione ^ 2))
summary(quad_gestazione)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze +
## Sesso + I(Gestazione^2), data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.0 -181.5 -12.9 165.8 2661.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4646.7158 898.6322 -5.171 2.52e-07 ***
## Gestazione -81.2309 49.7402 -1.633 0.10257
## Lunghezza 10.3502 0.3040 34.045 < 2e-16 ***
## Cranio 10.6376 0.4282 24.843 < 2e-16 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## SessoM 75.7563 11.2435 6.738 1.99e-11 ***
## I(Gestazione^2) 1.5168 0.6621 2.291 0.02206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16anova(quad_gestazione)## Analysis of Variance Table
##
## Response: Peso
## Df Sum Sq Mean Sq F value Pr(>F)
## Gestazione 1 241379330 241379330 3204.2973 < 2.2e-16 ***
## Lunghezza 1 206092456 206092456 2735.8660 < 2.2e-16 ***
## Cranio 1 48941070 48941070 649.6900 < 2.2e-16 ***
## N.gravidanze 1 731441 731441 9.7098 0.001854 **
## Sesso 1 3650790 3650790 48.4640 4.286e-12 ***
## I(Gestazione^2) 1 395323 395323 5.2479 0.022057 *
## Residuals 2491 187646731 75330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(quad_gestazione, mod4)## df AIC
## quad_gestazione 8 35149.63
## mod4 7 35152.89BIC(quad_gestazione, mod4)## df BIC
## quad_gestazione 8 35196.21
## mod4 7 35193.65Notiamo che per l’effetto d’interazione tra gestazione e sesso, l’aic e bic specificatamente aumentano in maniera notevole, facendoci capire che per lo scopo di prevedere il peso dei neonati, avere un modello più complicato con variabili simili a queste due non è necessario, ed anzi, peggioreremmo la previsione.
Per il modello con in aggiunta l’effetto quadratico sulla gestazione notiamo che l’AIC diminuisce di poco, mentre il BIC aumenta in maniera notevole, per cui io mi sento comunque di accettare il modello 4 originale così come è grazie alla sua semplicità.
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 circa 3258.615 grammi.
new_data = data.frame(Gestazione = 39, Lunghezza = media_lung, Cranio = media_cranio, N.gravidanze = 2, Sesso = "F")
results = predict.lm(mod4, newdata = new_data)
kable(results)| x |
|---|
| 3258.727 |
Grafici extra per l’interazione tra variabili
ggplot(data = df)+
geom_bar(aes(x = Gestazione, y = Peso), stat = "identity")+
theme_minimal()+
labs(
title = "Influenza della gestazione sul peso",
x = "Gestazione",
y = "Influenza sul Peso"
)
ggplot(data = df)+
geom_bar(aes(x = Fumatrici, y = Peso), stat = "identity")+
theme_minimal()+
scale_x_continuous(breaks = seq(0, 1))+
labs(
title = "Influenza dell'essere fumatrici sul peso",
x = "Fumatrici",
y = "Influenza sul Peso"
)
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 nelle “prime” settimane di gestazione questa non influenza troppo il peso, nelle settimane successive invece il peso riceve molta più influenza.
Vediamo che le donne fumatrici non hanno una grande influenza sul peso mentre le donne che non fumano si. Infine vediamo che (come dovremmo aspettarci normalmente), la lunghezza ed il peso sono molto influenti l’una sull’altra.