Modello Statistico per la Previsione del Peso Neonatale
Daniela Galliani
2024-11-19
Analisi e Modelizzazione
Analisi Preliminare
Il dataset dello studio è così composto:
- Età della madre: variabile quantitativa continua (tempo)
- Numero di gravidanze: variabile quantitativa discreta
- Fumo materno: variabile dicotomica (0= non fumatrice, 1=fumatrice)
- Durata della gravidanza: variabile quantitativa continua (tempo)
- Peso del neonato: variabile quantitativa continua
- Tipo di parto: variabile categorica nominale
- Sesso del neonato: variabile categorica nominale
Per eseguire un’analisi descrittiva preliminare del dataset, vengono
di seguito riportate le tabelle che sintetizzano gli indici di
posizione, variabilità e forma per le variabili quantitative.
library(kableExtra)
library(moments)
library(ggplot2)
dati <- read.csv("neonati.csv",
stringsAsFactors = T)
attach(dati)
calcola_moda <- function(x) {
freq <- table(x)
max_freq <- max(freq)
moda <- as.numeric(names(freq[freq == max_freq]))
return(moda)
}
metrics <- list(
"Anni Madre" = Anni.madre,
"N gravidanze" = N.gravidanze,
"Gestazione [settimane]" = Gestazione,
"Peso [gr]" = Peso,
"Lunghezza [cm]" = Lunghezza,
"Cranio [cm]" = Cranio
)
indici.posizione <- data.frame(
Misure = c("Min", "Moda", "Mediana", "Media","Max")
)
for (name in names(metrics)) {
data <- metrics[[name]]
indici.posizione[[name]] <- c(
formatC(min(data), format = "f", digits = 2),
formatC(calcola_moda(data), format = "f", digits = 2),
formatC(median(data), format = "f", digits = 2),
formatC(mean(data), format = "f", digits = 2),
formatC(max(data), format = "f", digits =2)
)
}
kable(indici.posizione, caption = "Indici di Posizione")%>%
kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Min | 0.00 | 0.00 | 25.00 | 830.00 | 310.00 | 235.00 |
| Moda | 30.00 | 0.00 | 40.00 | 3300.00 | 500.00 | 340.00 |
| Mediana | 28.00 | 1.00 | 39.00 | 3300.00 | 500.00 | 340.00 |
| Media | 28.16 | 0.98 | 38.98 | 3284.08 | 494.69 | 340.03 |
| Max | 46.00 | 12.00 | 43.00 | 4930.00 | 565.00 | 390.00 |
Si può notare come il minimo per la variabile Anni
Madre sia 0, questo è evidentemente un errore di battitura.
Vengono quindi identificate quelle osservazioni per cui l’età della
madre è inferiore a 13 anni e vengono eliminate dal dataset.
Viene di seguito riportato il summary aggiornato degli indici di
posizione.
dati<- dati[Anni.madre>=13,]
attach(dati)## I seguenti oggetti sono mascherati da dati (pos = 3):
##
## Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
## Ospedale, Peso, Sesso, Tipo.partometrics <- list(
"Anni Madre" = Anni.madre,
"N gravidanze" = N.gravidanze,
"Gestazione [settimane]" = Gestazione,
"Peso [gr]" = Peso,
"Lunghezza [cm]" = Lunghezza,
"Cranio [cm]" = Cranio
)
indici.posizione <- data.frame(
Misure = c("Min", "Moda", "Mediana", "Media","Max")
)
for (name in names(metrics)) {
data <- metrics[[name]]
indici.posizione[[name]] <- c(
formatC(min(data), format = "f", digits = 2),
formatC(calcola_moda(data), format = "f", digits = 2),
formatC(median(data), format = "f", digits = 2),
formatC(mean(data), format = "f", digits = 2),
formatC(max(data), format = "f", digits =2)
)
}
kable(indici.posizione, caption = "Indici di Posizione")%>%
kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Min | 13.00 | 0.00 | 25.00 | 830.00 | 310.00 | 235.00 |
| Moda | 30.00 | 0.00 | 40.00 | 3300.00 | 500.00 | 340.00 |
| Mediana | 28.00 | 1.00 | 39.00 | 3300.00 | 500.00 | 340.00 |
| Media | 28.19 | 0.98 | 38.98 | 3284.18 | 494.70 | 340.03 |
| Max | 46.00 | 12.00 | 43.00 | 4930.00 | 565.00 | 390.00 |
indici.variabilita <- data.frame(
Misure = c("Range", "IQR", "Varianza", "Deviazione Standard", "Coefficiente di variazione [%]")
)
CV <- function(x){
return(sd(x)/mean(x)*100)
}
for (name in names(metrics)) {
data <- metrics[[name]]
indici.variabilita[[name]] <- c(
formatC(diff(range(data)), format = "f", digits = 2),
formatC(IQR(data), format = "f", digits = 2),
formatC(var(data), format = "f", digits = 2),
formatC(sd(data), format = "f", digits = 2),
formatC(CV(data), format = "f", digits = 2)
)
}
# Calcola indici di forma
indici.forma <- data.frame(
Misure = c("Skewness", "Kurtosis")
)
for (name in names(metrics)) {
data <- metrics[[name]]
indici.forma[[name]] <- c(
formatC(skewness(data),format = "f",digits=3),
formatC(kurtosis(data)-3,format = "f",digits=3)
)
}
kable(indici.variabilita, caption = "Indici di Variabilità")%>%
kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Range | 33.00 | 12.00 | 18.00 | 4100.00 | 255.00 | 155.00 |
| IQR | 7.00 | 1.00 | 2.00 | 630.00 | 30.00 | 20.00 |
| Varianza | 27.22 | 1.64 | 3.49 | 275865.90 | 693.21 | 269.93 |
| Deviazione Standard | 5.22 | 1.28 | 1.87 | 525.23 | 26.33 | 16.43 |
| Coefficiente di variazione [%] | 18.51 | 130.50 | 4.79 | 15.99 | 5.32 | 4.83 |
kable(indici.forma, caption = "Indici di Forma")%>%
kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Skewness | 0.151 | 2.513 | -2.065 | -0.647 | -1.515 | -0.785 |
| Kurtosis | -0.106 | 10.982 | 8.256 | 2.029 | 6.481 | 2.945 |
frequenza_percentuale <- table(N.gravidanze) / length(N.gravidanze) * 100
frequenza_percentuale_arrotondata <- sprintf("%.2f", frequenza_percentuale)
frequenza_df <- data.frame(
Gravidanze = names(frequenza_percentuale),
Percentuale = frequenza_percentuale_arrotondata
)
kable(frequenza_df, col.names = c("Numero di Gravidanze", "Percentuale (%)"), caption = "Distribuzione di frequenza del numero di gravidanze") %>%
kable_styling(position = "center")| Numero di Gravidanze | Percentuale (%) |
|---|---|
| 0 | 43.84 |
| 1 | 32.71 |
| 2 | 13.61 |
| 3 | 6.00 |
| 4 | 1.92 |
| 5 | 0.84 |
| 6 | 0.44 |
| 7 | 0.04 |
| 8 | 0.32 |
| 9 | 0.08 |
| 10 | 0.12 |
| 11 | 0.04 |
| 12 | 0.04 |
Dalle tabelle sopra riportate si possono fare le seguenti considerazioni
per variabile:
- Anni madre: ha una variabilità non trascurabile
(circa il 18.5%) con un range di 33 anni. La moda indica che l’età più
comune tra le madri sia 30 anni, mentre la mediana e la media sono
leggermente inferiori (28 anni). Sia la Skewness che la Curtosi si
aggirano intorno allo zero indicando una distribuzione della variabile
pressocchè normale e simmetrica;
- Numero gravidanze: il numero di gravidanze mostra
una variazione del 130.5% con un range di 12 che va da 0 gravidanze
precedenti a 12. Dal valore della moda, della mediana e della media si
deduce che il numero di partecipanti allo studio per cui questa è la
prima gravidanza è preponderante. Come si vede dalla distribuzione di
frequenza sopra riportata di questa variabile, quasi il 44% dei soggetti
è alla prima gravidanza e circa il 90% dei soggetti ha avuto massimo due
gravidanze. Numeri molto elevati di gravidanze (11 e 12) compaiono nel
dataset con una sola osservazione ciascuno, nonostante sia raro, questi
dati sono credibili quindi queste osservazioni vengono mantenute nel
dataset. Questa variabile ha un’asimmetria positiva e una distribuzione
leptocurtica;
- Gestazione: le settimane di gestazione presentano
una variabilità relativamente bassa (intorno al 5%), con range che va
dalle 25 settimane (circa 6 mesi, nascita prematura) a 43 settimane
(quasi 11 mesi, parto tardivo) di gestazione, con una moda di 40
settimane e una mediana e una media comparabili tra loro e leggermente
inferiori (circa 39 settimane). La variabile ha un’asimmetria negativa e
una distribuzione leptocurtica;
- Peso: il peso del neonato ha una variabilità non
trascurabile (circa il 16%) e un range di circa 4,1 kg, da 830 g
(neonato sottopeso, sotto il 3° percentile delle curve di crescita per
entrambi i sessi) a 4930 gr( neonato sovrappeso, sopra il 97° percentile
delle curve di crescita per entrambi i sessi). Ha una distribuzione
lievemente asimmetrica negativa e leptocurtica;
- Lunghezza: ha una bassa variabilità (coefficiente
di variazioni di circa il 5%). Presenta un valore minimo di 310 cm
(inferiore al 3° percentile delle curve di crescita per entrambi i
sessi) e un valore massimo di 565 cm (sopra il 97° percentile delle
curve di crescita per entrambi i sessi). Ha una distribuzione
asimmetrica negativa e leptocurtica;
- Cranio: ha una bassa variabilità (coefficiente di
variazioni di circa il 5%). Presenta un valore minimo di 235 cm e un
valore massimo di 390 cm. Ha una distribuzione che presenta una leggera
asimmetria negativa e leptocurtica.
freq_distr <- function(x){
ni <- table(x)
fi <- formatC(ni/length(x), format="f", digits=2)
distr_freq <- as.data.frame(cbind(ni,fi))
return(distr_freq)
}
distr_freq_fumo <- freq_distr(Fumatrici)
distr_freq_parto <- freq_distr(Tipo.parto)
distr_freq_osp <- freq_distr(Ospedale)
distr_freq_sesso <- freq_distr(Sesso)
kable(distr_freq_fumo, caption = "Distribuzione di frequenza per la variabile Fumo materno")%>%
kable_styling(position = "center")| ni | fi | |
|---|---|---|
| 0 | 2394 | 0.96 |
| 1 | 104 | 0.04 |
kable(distr_freq_parto, caption = "Distribuzione di frequenza per la variabile Tipo di parto")%>%
kable_styling(position = "center")| ni | fi | |
|---|---|---|
| Ces | 728 | 0.29 |
| Nat | 1770 | 0.71 |
kable(distr_freq_osp, caption = "Distribuzione di frequenza per la variabile Ospedale di nascita")%>%
kable_styling(position = "center")| ni | fi | |
|---|---|---|
| osp1 | 816 | 0.33 |
| osp2 | 848 | 0.34 |
| osp3 | 834 | 0.33 |
kable(distr_freq_sesso, caption = "Distribuzione di frequenza per la variabile Sesso del neonato")%>%
kable_styling(position = "center")| ni | fi | |
|---|---|---|
| F | 1255 | 0.50 |
| M | 1243 | 0.50 |
Per le variabili qualitative sono state invece riportate le
distribuzioni di frequenza:
- Fumo materno : la variabile assume valore 0 quando
la madre è non fumatrice come nel 96% dei casi, 1 altrimenti;
- Tipo di parto: si ha una prevalenza di parto
naturale con il 71% delle osservazioni;
- Ospedale di nascita: le partecipanti allo studio
sono state selezionate in maniera eguale tra i 3 ospedali;
- Sesso del neonato: si ha un numero uguale di
neonati per ambo i sessi.
Per valutare come il Peso del neonato varia in
relazione alle variabili materne, vengono di seguito riportati dei
grafici.
Viene indagato come il peso del neonato varia all’aumentare delle
settimane di gestazione, prendendo in considerazione
anche la variabile Fumatrici.
Si nota come la relazione tra le settimane di gestazione e il peso del
neonato sia prevedibilmente abbastanza lineare, all’aumentare del tempo
di gestazione aumenta anche il peso del neonato. Si vede anche che a
parità di settimane di gestazione se la madre è una fumatrice il peso
del neonato è inferiore.
Dal boxplot non si nota invece una dipendenza tra il numero di settimane
di gestazione e la variabile fumo.
library(dplyr)##
## Caricamento pacchetto: 'dplyr'## Il seguente oggetto è mascherato da 'package:kableExtra':
##
## group_rows## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(scales)
dati_aggregati_gest <- dati %>%
group_by(Gestazione,Fumatrici) %>%
summarise(PesoMedio = mean(Peso, na.rm = TRUE))## `summarise()` has grouped output by 'Gestazione'. You can override using the
## `.groups` argument.ggplot(dati_aggregati_gest, aes(x = Gestazione, y = PesoMedio, color = as.factor(Fumatrici))) +
geom_line() +
labs(title = "Peso Medio del Neonato rispetto alle settimane di gestazione",
x = "Settimane di gestazione",
y = "Peso Medio del Neonato (g)") +
scale_color_manual(values = c("0" = "darkblue", "1" = "red"),
labels = c("Non Fumatrici", "Fumatrici")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
legend.title = element_blank(),
legend.text = element_text(size = 7)) +
theme(panel.grid.major = element_line(color = "lightgrey", size = 0.5),
panel.grid.minor = element_blank())## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
cat("<div style='text-align: center;'>")## <div style='text-align: center;'>cat("</div>")## </div>ggplot(data = dati, aes(x = factor(Fumatrici), y = Gestazione)) +
geom_boxplot() +
labs(
title = "Settimane di Gestazione in base al Fumo Materno",
x = "Fumatrice (0 = No, 1 = Sì)",
y = "Settimane di Gestazione",
fill = "Fumatrice"
) +
theme_minimal() + # Cambia il tema per un aspetto più moderno
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), # Centra e formatta il titolo
axis.title.x = element_text(size = 12, face = "bold"), # Formatta il titolo dell'asse X
axis.title.y = element_text(size = 12, face = "bold"), # Formatta il titolo dell'asse Y
legend.title = element_text(size = 12, face = "bold") # Formatta il titolo della legenda
)
Di seguito lo scatter plot riporta la distribuzione del peso del
neonato in funzione dell’età della madre. Non
si evidenzia un particolare andamento, i pesi maggiori e minori si
concentrano nelle età tra 25 anni e 35 anni dove si hanno però anche il
maggior numero di osservazioni. Anche la variabile
fumatrici viene presa in considerazione, ma i pesi dei
neonati da madri fumatrici o non non sembrano differenziarsi
significativamente.
ggplot(dati, aes(x = Anni.madre, y = Peso, color = as.factor(Fumatrici))) +
geom_point() +
labs(title = "Peso del Neonato vs Età della Madre",
x = "Età della Madre",
y = "Peso del Neonato (g)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.title = element_blank(),
legend.text = element_text(size = 7))+
theme(panel.grid.major = element_line(color = "grey", size = 0.5),
panel.grid.minor = element_blank())+
scale_color_manual(values = c("0" = "darkblue", "1" = "red"),
labels = c("Non Fumatrici", "Fumatrici"))
cat("<div style='text-align: center;'>")## <div style='text-align: center;'>cat("</div>")## </div>
Di seguito viene invece mostrato l’andamento del peso del
neonato in relazione al numero di gravidanze.
Il peso del neonato non sembra variare significativamente in relazione
al numero di gravidanze.
ggplot(dati, aes(x = as.factor(N.gravidanze), y = Peso)) +
geom_boxplot() +
labs(title = "Distribuzione del Peso del Neonato per Numero di Gravidanze",
x = "Numero di Gravidanze",
y = "Peso del Neonato (g)") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
cat("<div style='text-align: center;'>")## <div style='text-align: center;'>cat("</div>")## </div>
In ultimo viene mostrato il boxplot che mette in relazione il
peso del neonato con la presenza o l’assenza del
fumo materno. Anche in questo caso non si vede una
chiara relazione tra le due variabili.
ggplot(dati, aes(x = as.factor(Fumatrici), y = Peso)) +
geom_boxplot() +
labs(title = "Distribuzione del Peso del Neonato rispetto allo Stato di Fumatrice",
x = "Fumatrici",
y = "Peso del Neonato (g)") +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
cat("<div style='text-align: center;'>")## <div style='text-align: center;'>cat("</div>")## </div>
Prima di proseguire valutiamo se le medie del peso e della lunghezza del
campione oggetto di studio sono significativamente diverse da quelle
della popolazione attraverso un test t di student. Per i valori di
riferimento per la popolazione si prendono le medie per neonati
italiani:
- Peso: M = 3200 g e F = 3100 g;
- Lunghezza: M = 49.9 cm e F = 49.1 cm;
peso_medio_M <- 3200
maschi <- dati %>% filter(Sesso== "M")
t.test(x= maschi$Peso,
mu = peso_medio_M)##
## One Sample t-test
##
## data: maschi$Peso
## t = 14.883, df = 1242, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 3200
## 95 percent confidence interval:
## 3381.012 3435.979
## sample estimates:
## mean of x
## 3408.496peso_medio_F <- 3100
femmine <- dati %>% filter(Sesso== "F")
t.test(x= femmine$Peso,
mu = peso_medio_F)##
## One Sample t-test
##
## data: femmine$Peso
## t = 4.1085, df = 1254, p-value = 4.241e-05
## alternative hypothesis: true mean is not equal to 3100
## 95 percent confidence interval:
## 3131.904 3190.219
## sample estimates:
## mean of x
## 3161.061lunghezza_media_M <- 499
t.test(x= maschi$Lunghezza,
mu = lunghezza_media_M)##
## One Sample t-test
##
## data: maschi$Lunghezza
## t = 0.98965, df = 1242, p-value = 0.3225
## alternative hypothesis: true mean is not equal to 499
## 95 percent confidence interval:
## 498.3369 501.0131
## sample estimates:
## mean of x
## 499.675lunghezza_media_F <- 491
t.test(x= femmine$Lunghezza,
mu = lunghezza_media_F)##
## One Sample t-test
##
## data: femmine$Lunghezza
## t = -1.5894, df = 1254, p-value = 0.1122
## alternative hypothesis: true mean is not equal to 491
## 95 percent confidence interval:
## 488.2387 491.2896
## sample estimates:
## mean of x
## 489.7641per_sovrappeso <- nrow(dati[dati$Peso>4000,])/nrow(dati)*100
per_sottopeso <- nrow(dati[dati$Peso<2500,])/nrow(dati)*100I test t di student per la variabile lunghezza non rifiutano
l’ipotesi nulla quindi la media del campione in esame non è
significativamente diversa da quella della popolazione.
Per i t test eseguiti sul peso viene rifiutata l’ipotesi nulla, quindi
le medie sono significativamente diverse da quelle prese di riferimento
per la popolazione. Considerando i seguenti valori di riferimento:
- sottopeso : 2,5 kg
- sovrappeso : 4 kg
e calcolando le percentuali per nascite premature e tardive nel
nostro dataset queste risultano essere rispettivamente intorno al 5.5% e
6.6%. Mentre la percentuale di nascite premature è in linea rispetto a
quella italiana (7%-9%) quella per nascite tardive è nettamente
superiore (media italiana ~ 1-2%). Questa potrebbe essere la ragione per
cui il test t di student evidenzia una differenza significativa dalla
media della popolazione.
Potrebbe essere utile anche valutare se ci sono differenze
significative per peso, lunghezza e diametro del cranio tra maschi e
femmine. Anche quest’analisi viene fatta con test t di student.
t.test(x= femmine$Peso,y=maschi$Peso)##
## Welch Two Sample t-test
##
## data: femmine$Peso and maschi$Peso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -287.4841 -207.3844
## sample estimates:
## mean of x mean of y
## 3161.061 3408.496t.test(x= femmine$Lunghezza,y=maschi$Lunghezza)##
## Welch Two Sample t-test
##
## data: femmine$Lunghezza and maschi$Lunghezza
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.939001 -7.882672
## sample estimates:
## mean of x mean of y
## 489.7641 499.6750t.test(x= femmine$Cranio,y=maschi$Cranio)##
## Welch Two Sample t-test
##
## data: femmine$Cranio and maschi$Cranio
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean of x mean of y
## 337.6231 342.4586In questo caso tutti i t test danno lo stesso risultato, ovvero che
per le tre variabili prese in considerazione c’è una differenza
statisticamente significativa tra le medie dei valori per i due
sessi.
L’ultima analisi prima di proseguire con la creazione del modello
riguarda l’associazione tra tipo di parto e ospedale, quest’analisi
viene eseguita con un test chi-quadro.
tabella_contingenza <- table(dati$Ospedale, dati$Tipo.parto)
print(tabella_contingenza)##
## Ces Nat
## osp1 242 574
## osp2 254 594
## osp3 232 602test_chi_quadro <- chisq.test(tabella_contingenza)
print(test_chi_quadro)##
## Pearson's Chi-squared test
##
## data: tabella_contingenza
## X-squared = 1.083, df = 2, p-value = 0.5819Il test chi quadro non evidenzia un’associazione statisticametne significativa tra tipo di parto e ospedale.
Creazione del Modello di Regressione
Visualizziamo la matrice di correlazione tra le variabili.
dati <- dati %>% select(Peso, everything())
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- (cor(x, y))
txt <- format(c(r, 1), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 1.5)
}
pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
Tralasciando per ora le variabili qualitative
(Fumatrici, Tipo parto,
Ospedale e Sesso) analizziamo i
coefficienti di correlazione lineare e i relativi scatterplot per le
variabili restanti.
- Anni madre : ha un coefficiente di correlazione lineare leggermente negativo, ma molto basso, infatti anche dallo scatterplot si vede che la nuvola di punti non ha una chiara tendenza;
- N.gravidanze : anche in questo caso lo scatterplot mostra una nuvola di punti molto sparsa e il coefficiente di correlazione è prossimo allo 0;
- Gestazione : questa variabile ha una correlazione positiva abbastanza importante con il peso, si vede infatti sia dal coefficiente di correlazione positivo che dall’andamento crescente dello scatterplot. In questo caso si intuisce anche un possibile andamento quadratico;
- Lunghezza : questa variabile ha una forte correlazione con il peso, infatti mostra un coefficiente di correlazione positivo alto (0.8) e una nuvola di punti concentrata intorno alla bisettrice con andamento lineare crescente;
- Cranio : anche in questo caso si ha una forte correlazione positiva, ma lo scatterplot mostra un andamento leggermente quadratico.
Per le variabili qualitative vengono invece riportati dei boxplot.
Non viene fatta alcuna analisi per le variabili Tipo di
parto ed Ospedali in quanto la loro influenza
sulla variabile output peso può essere esclusa logicamente a
priori.
par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Fumatrici)
t.test(Peso~Fumatrici)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.0362, df = 114.12, p-value = 0.3023
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.5076 145.3399
## sample estimates:
## mean in group 0 mean in group 1
## 3286.262 3236.346Dal boxplot non sembra che la variabile variabile Fumatrici abbia influenza sulla variabile peso e questo è confermato dal test t di student che con un p-value di 0.3 non rifiuta l’ipotesi nulla che le due medie non siano significativamente diverse.
par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Sesso)
t.test(Peso~Sesso)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496
Relativamente al peso in funzione del sesso, si vede invece dai boxplot
che i due pesi sono diversi tra maschi e femmine, confermato dal test t
di student che rifiuta l’ipotesi nulla. C’è quindi una differenza
statisticamente significativa tra la media del peso dei neonati maschi e
quella del peso delle neonate femmine.
Creiamo ora il primo modello che considera tutte le variabili, a parte
Tipo di parto e Ospedali escluse
logicamente a priori.
mod1 <- lm(Peso~.-Tipo.parto-Ospedale, data= dati )
summary(mod1)##
## Call:
## lm(formula = Peso ~ . - Tipo.parto - Ospedale, data = dati)
##
## 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 ***
## Anni.madre 0.8803 1.1491 0.766 0.444
## N.gravidanze 11.3789 4.6767 2.433 0.015 *
## Fumatrici -30.3958 27.6080 -1.101 0.271
## Gestazione 32.9472 3.8288 8.605 < 2e-16 ***
## Lunghezza 10.2316 0.3011 33.979 < 2e-16 ***
## Cranio 10.5198 0.4271 24.633 < 2e-16 ***
## 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-16Ottimizzazione del modello
Dal primo modello calcolato sembra che le variabili Anni
madre e Fumatrici non abbiano impatto sulla
variabile output. La variabile del fumo in gravidanze è una delle
variabili di interesse nello studio, quindi per il momento rimuovo solo
la variabile anni madre e calcolo un nuovo modello.
mod2 <- update(mod1, ~.-Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso, data = dati)
##
## 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 ***
## N.gravidanze 12.6996 4.3470 2.921 0.00352 **
## Fumatrici -30.5728 27.6048 -1.108 0.26818
## Gestazione 32.6437 3.8079 8.573 < 2e-16 ***
## Lunghezza 10.2309 0.3011 33.979 < 2e-16 ***
## Cranio 10.5366 0.4265 24.707 < 2e-16 ***
## 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-16La rimozione della variabile non influisce sul R2. La variabile
Fumatrici continua a non avere un effetto influente,
per fare un ulteriore prova prima della rimozione della variabile dal
modello provo a considerare l’effetto simultaneo di
Fumatrici e Gestazione.
mod3 <- update(mod2, ~.+Fumatrici:Gestazione)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + Fumatrici:Gestazione, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.83 -181.58 -16.95 163.76 2635.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6699.6947 136.7290 -49.000 < 2e-16 ***
## N.gravidanze 12.7452 4.3470 2.932 0.0034 **
## Fumatrici 795.7016 757.5315 1.050 0.2936
## Gestazione 33.2007 3.8418 8.642 < 2e-16 ***
## Lunghezza 10.2252 0.3011 33.957 < 2e-16 ***
## Cranio 10.5313 0.4265 24.693 < 2e-16 ***
## SessoM 78.7436 11.2241 7.016 2.94e-12 ***
## Fumatrici:Gestazione -21.0469 19.2830 -1.091 0.2752
## ---
## 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.7273, Adjusted R-squared: 0.7265
## F-statistic: 948.6 on 7 and 2490 DF, p-value: < 2.2e-16L’R2 non subisce variazioni, ma entrambe le variabili che considerano
l’effetto del fumo non risultano significative, questo potrebbe
dipendere dal numero di osservazioni relative a madri fumatrici molto
inferiore rispetto a quello delle madri non fumatrici. Si rimuovono
quindi entrambe le variabili.
mod4 <- update(mod3, ~.-Fumatrici-Fumatrici:Gestazione)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## 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 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## 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 ***
## 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-16L’R2 non subisce variazioni e in questo modello tutti regressori sono
significativi.
Per la variabile gestazione potrebbe aver senso considerare il suo
effetto quadratico, come evidenziato dalla matrice di
correlazione.
mod5 <- update(mod4, ~.+I(Gestazione^2))
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2), data = dati)
##
## 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 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## 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 ***
## 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-16Si vede un aumento di R2 trascurabile e la variabile
Gestazione perde la sua influenza sul modello,
l’effetto quadratico di Gestazione non viene quindi considerato.
Anche la variabile Cranio potrebbe avere un effetto
quadratico, provo quindi a considerarla nel modello.
mod6 <- update(mod5, ~.-I(Gestazione^2)+I(Cranio^2))
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Cranio^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.6 -179.4 -14.8 163.4 2622.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.10118 1151.77280 0.073 0.94180
## N.gravidanze 12.76356 4.31259 2.960 0.00311 **
## Gestazione 38.90540 3.93291 9.892 < 2e-16 ***
## Lunghezza 10.48745 0.30157 34.776 < 2e-16 ***
## Cranio -31.79371 7.16973 -4.434 9.63e-06 ***
## SessoM 73.10236 11.16590 6.547 7.11e-11 ***
## I(Cranio^2) 0.06262 0.01059 5.915 3.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared: 0.7308, Adjusted R-squared: 0.7301
## F-statistic: 1127 on 6 and 2491 DF, p-value: < 2.2e-16Anche in questo modello tutti i regressori considerati sono
significativi e l’R2 aumenta leggermente.
Selezione del modello ottimale
Procedo ora alla selezione del modello migliore tra i sei modelli
creati.
Utilizzo i criteri BIC e AIC per scegliere il modello migliore.
AIC(mod1,mod2,mod3,mod4,mod5,mod6)## df AIC
## mod1 9 35155.07
## mod2 8 35153.66
## mod3 9 35154.46
## mod4 7 35152.89
## mod5 8 35149.63
## mod6 8 35120.04BIC(mod1,mod2,mod3,mod4,mod5,mod6)## df BIC
## mod1 9 35207.48
## mod2 8 35200.24
## mod3 9 35206.87
## mod4 7 35193.65
## mod5 8 35196.21
## mod6 8 35166.63Per entrambi i criteri il modello 6 è quello migliore (AIC e BIC
minore), seguito dal modello 4. Gli R2 dei due modelli differivano di un
meno di un punto percentuale, utilizzo il test anova per verificare se
c’è stato un aumento significativo della varianza spiegata con
l’aggiunta di una variabile nel modello 6.
anova(mod6,mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Cranio^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 185437522
## 2 2492 188042054 -1 -2604532 34.987 3.775e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C’è un aumento significativo di varianza spiegata.
Verifico ora la multicollinearità tra le variabili.
library(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.4.2## Caricamento del pacchetto richiesto: carData## Warning: il pacchetto 'carData' è stato creato con R versione 4.4.2##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recodevif(mod6)## N.gravidanze Gestazione Lunghezza Cranio Sesso I(Cranio^2)
## 1.023611 1.812253 2.114666 465.421862 1.045890 452.902551vif(mod4)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023462 1.669779 2.075747 1.624568 1.040184In presenza del termine quadratico i VIF aumentano leggermente per le
variabili comuni a entrambi i modelli, ma hanno un valore >> 5 per
la variabile cranio e il suo termine quadratico. Il modello 6 ha un
aumento di meno di un punto percentuale su R2, quindi viene preferito il
più semplice (modello 4).
Analisi qualità del modello
Viene eseguita ora un’analisi dei residui.
par(mfrow=c(2,2))
plot(mod4)
I punti dovrebbero presentarsi sparsi casualmente attorno alla media di
0 per il primo grafico, questo è abbastanza vero anche se si nota una
leggera curvatura e c’è l’osservazione 1551 che si distanzia
notevolmente dalle altre.
Nel secondo grafico è mostrata la relazione tra i residui e i quantili
di una normale. I residui di dispongono sulla bisettrice a parte per le
code che si distanziano leggermente. Anche questo caso l’osservazione
1551 si distanzia dalle altre.
Nel terzo grafico la nuvola di punti si dispone abbastanza casualmente
intorno ad un certo valore, a parte per l’osservazione 1551, indicando
una varianza abbastanza costante. Anche qui si nota una leggera
curvatura.
Nel quarto grafico si vede come l’osservazione 1551 soprassi la soglia
di avvertimento a 0.5, rimanendo comunque sotto la soglia di allarme a
1.
Si eseguono ora i test di Breusch-Pagan, Durbin-Watson e
Shapiro-Wilk.
lmtest::bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 90.297, df = 5, p-value < 2.2e-16lmtest::dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(mod4$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4$residuals
## W = 0.97414, p-value < 2.2e-16plot(density(residuals(mod4)))
Il test di Breusch-Pagan rifiuta l’ipotesi nulla quindi si rifiuta
l’ipotesi nulla di omoschedasticità, questo potrebbe influenzare
l’efficienza delle stime dei coefficienti.
Non viene invece rifiutata l’ipotesi nulla di Durbin-Watson, quindi i
residui non sono autocorrelati.
Anche il test di Shapiro-Wilk rifiuta l’ipotesi nulla di normalità, ma
dalla visualizzazione grafica la distribuzione dei residui si avvicina a
quella normale.
Consideriamo ora i valori di leva e gli outliers.
#leverage
lev<-hatvalues(mod4)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
plot(lev)
abline(h=soglia,col=2)
n_lev <- length(lev[lev>soglia])
#outliers
plot(rstudent(mod4))
abline(h=c(-2,2))
car::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-03#distanza di cook
cook<-cooks.distance(mod4)
plot(cook,ylim = c(0,1)) 
Sono stati identificati sia valori di leva (152) che valori outlier (3),
ma l’unica osservazione che supera il primo valore di avvertimento della
distanza di Cook è l’osservazione 1551. Facendo una valutazione
sull’osservazione questa non può essere esclusa dal modello perchè non
presenta dei valori logicamente impossibili.
Si eseguono i test statistici sui residui del modello 6 per un’analisi
comparativa.
lmtest::bptest(mod6)##
## studentized Breusch-Pagan test
##
## data: mod6
## BP = 90.755, df = 6, p-value < 2.2e-16lmtest::dwtest(mod6)##
## Durbin-Watson test
##
## data: mod6
## DW = 1.9523, p-value = 0.1165
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(mod6$residuals)##
## Shapiro-Wilk normality test
##
## data: mod6$residuals
## W = 0.97436, p-value < 2.2e-16lev6<-hatvalues(mod6)
p6<-sum(lev6)
n6<-length(lev6)
soglia6=2*p6/n6
plot(lev6)
abline(h=soglia6,col=2)
n_lev6 <- length(lev6[lev6>soglia])
#outliers
plot(rstudent(mod6))
abline(h=c(-2,2))
car::outlierTest(mod6)## rstudent unadjusted p-value Bonferroni p
## 1551 10.052081 2.4894e-23 6.2185e-20
## 155 5.271647 1.4682e-07 3.6676e-04
## 1306 4.805789 1.6331e-06 4.0795e-03
## 1694 4.277127 1.9646e-05 4.9076e-02#distanza di cook
cook6<-cooks.distance(mod6)
plot(cook6,ylim = c(0,1)) 
Anche in questo caso non viene rifiutata l’ipotesi nulla per il test di
Durbin-Watson, ma viene rifiutata per i test di Breusch-Pagan e
Shapiro-Wilk. Con il modello 6 ottengo 4 outliers e 212 valori di leva,
l’osservazione 1551 supera ancora una volta il valore di 0.5 per la
distanza di cook ma è inferiore al valore di allarme.
L’aumento della complessità del modello non migliora sensibilmente l’R2
e non ha effetto sull’analisi dei residui, viene quindi confermata la
scelta del modello 4.
Previsioni e risultati
Facciamo una previsione considerando questi dati:
- Sesso = F
- Numero gravidanze = 2 (terza gravidanza)
- Settimane Gestazione = 39
- Fumatrice = NO
Non essendoci forniti dati per la lunghezza o il diametro del cranio
si usano come dati le medie di lunghezza e diametro per le
neonate.
dati_lunghezza <- dati %>%
group_by(Sesso) %>%
summarise(LunghezzaMedia = mean(Lunghezza, na.rm = TRUE))
dati_diametro <- dati %>%
group_by(Sesso) %>%
summarise(DiametroMedio = mean(Cranio, na.rm = TRUE))La lunghezza media per le neonate è 489,7641 mm e il diametro del
cranio 337,6231 mm.
new_data <- data.frame(N.gravidanze=2,Gestazione=39,Lunghezza=489.7641, Cranio=339.6231,Sesso="F")
peso_new <- predict(mod4,new_data)Il modello prevede un peso della neonata di 3203.92 g.
Visualizzazioni
Visualizziamo graficamente i risultati del modello per mostrare le
relazioni tra le variabili.
Viene di seguito mostrata la relazione tra settimane di gestazione e
peso del neonato, tenendo in considerazione anche il sesso.
ggplot(data = dati) +
geom_point(aes(x = Gestazione, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Gestazione, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Settimane di Gestazione e Peso del Neonato",
x = "Settimane di Gestazione",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Si può notare come la retta di regressione per maschi e femmine abbia la
stessa tendenza crescente positiva, ma che i pesi delle neonate siano
inferiori a quelli dei neonati. Le osservazioni per una durata minore
della gestazione si trovano al di sotto delle retta di regressione
quindi il modello potrebbe sovrastimare il peso del neonato quando
questo nasce prematuramente. La nuvola di punti sembra invece più
concentrata sulla retta di regressioni per gravidanze più lunghe.
Indaghiamo ora la relazione tra il numero di gravidanze e il peso del
neonato.
ggplot(data = dati) +
geom_point(aes(x = N.gravidanze, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = N.gravidanze, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Numero di Gravidanze e Peso del Neonato",
x = "Numero di gravidanze",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Il numero di gravidanze non sembra avere influenza sul peso del neonato.
Ancora una volta il peso delle neonate risulta essere inferiore a quello
dei neonati.
Visualizziamo la relazione tra lunghezza e peso.
ggplot(data = dati) +
geom_point(aes(x = Lunghezza, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Lunghezza, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Lunghezza e Peso del Neonato",
x = "Lunghezza (cm)",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Anche in questo caso, come per le settimane di gestazione, si ha una
dipendenza lineare positiva della variabile output per ambo i sessi. Si
nota come sia meno visibile la differenza di peso tra maschi e femmine
al variare della lunghezza.
Analizziamo anche la relazione tra il diametro del cranio e il peso del
neonato.
ggplot(data = dati) +
geom_point(aes(x = Cranio, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Cranio, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Diametro del Cranio e Peso del Neonato",
x = "Diametro del cranio",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Si può osservare in questo caso la stessa tendenza evidenziata per le
settimane di gestazione. Si nota una crescita lineare, ma il modello
sembra sovrastimare i pesi per diametri del cranio minori. Nuovamente si
nota la differenza in peso tra maschi e femmine.