progetto_neonati
filo
2024-11-17
Carico il dataset dei neonati dal percorso specificato
neonati <- read.csv("/Users/filomena/Desktop/1 esercizio con r/neonati.csv")Installazione e caricamento delle librerie necessarie per l’analisi
necessary_packages <- c("ggplot2", "broom", "caret", "dplyr", "corrplot", "tidyr", "MASS",
"mgcv", "lattice", "rpart", "randomForest", "xgboost", "car")
installed_packages <- rownames(installed.packages())
for (pkg in necessary_packages) {
if (!pkg %in% installed_packages) {
install.packages(pkg, dependencies = TRUE)
}
}
lapply(necessary_packages, library, character.only = TRUE)## Loading required package: lattice##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## corrplot 0.95 loaded##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## margin##
## Attaching package: 'xgboost'## The following object is masked from 'package:dplyr':
##
## slice## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## [[1]]
## [1] "ggplot2" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "broom" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "caret" "lattice" "broom" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "dplyr" "caret" "lattice" "broom" "ggplot2" "stats"
## [7] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "corrplot" "dplyr" "caret" "lattice" "broom" "ggplot2"
## [7] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[6]]
## [1] "tidyr" "corrplot" "dplyr" "caret" "lattice" "broom"
## [7] "ggplot2" "stats" "graphics" "grDevices" "utils" "datasets"
## [13] "methods" "base"
##
## [[7]]
## [1] "MASS" "tidyr" "corrplot" "dplyr" "caret" "lattice"
## [7] "broom" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [13] "datasets" "methods" "base"
##
## [[8]]
## [1] "mgcv" "nlme" "MASS" "tidyr" "corrplot" "dplyr"
## [7] "caret" "lattice" "broom" "ggplot2" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "mgcv" "nlme" "MASS" "tidyr" "corrplot" "dplyr"
## [7] "caret" "lattice" "broom" "ggplot2" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "rpart" "mgcv" "nlme" "MASS" "tidyr" "corrplot"
## [7] "dplyr" "caret" "lattice" "broom" "ggplot2" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "randomForest" "rpart" "mgcv" "nlme" "MASS"
## [6] "tidyr" "corrplot" "dplyr" "caret" "lattice"
## [11] "broom" "ggplot2" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "xgboost" "randomForest" "rpart" "mgcv" "nlme"
## [6] "MASS" "tidyr" "corrplot" "dplyr" "caret"
## [11] "lattice" "broom" "ggplot2" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "car" "carData" "xgboost" "randomForest" "rpart"
## [6] "mgcv" "nlme" "MASS" "tidyr" "corrplot"
## [11] "dplyr" "caret" "lattice" "broom" "ggplot2"
## [16] "stats" "graphics" "grDevices" "utils" "datasets"
## [21] "methods" "base"Imposto un seme per la riproducibilitÃ
set.seed(123)Questa riga di codice serve per garantire la riproducibilità dei risultati. Impostando un valore di seme (in questo caso 123), si garantisce che i risultati siano gli stessi ogni volta che il codice viene eseguito. È una pratica fondamentale, soprattutto in ambiti di analisi statistica o machine learning, per ottenere risultati consistenti e confrontabili. 1. Analisi Preliminare dei Dati
str(neonati) # Struttura del dataset## 'data.frame': 2500 obs. of 10 variables:
## $ Annimadre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ Ngravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipoparto : chr "Nat" "Nat" "Nat" "Nat" ...
## $ Ospedale : chr "osp3" "osp1" "osp2" "osp2" ...
## $ Sesso : chr "M" "F" "M" "M" ...names(neonati) # Nomi delle variabili## [1] "Annimadre" "Ngravidanze" "Fumatrici" "Gestazione" "Peso"
## [6] "Lunghezza" "Cranio" "Tipoparto" "Ospedale" "Sesso"head(neonati) # Prime righe del dataset## Annimadre Ngravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipoparto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 FcolSums(is.na(neonati)) # Controllo valori mancanti## Annimadre Ngravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipoparto Ospedale Sesso
## 0 0 0 0Converto le variabili categoriche in fattori
variabili_categoriche <- c("Fumatrici", "Tipoparto", "Ospedale", "Sesso")
neonati <- neonati %>%
mutate(across(all_of(variabili_categoriche), as.factor))variabili_categoriche è un vettore che contiene i nomi delle colonne del dataset neonati che devono essere trasformate. La a funzione mutate() di dplyr permette di modificare o aggiungere colonne al dataset I fattori sono essenziali per la corretta gestione delle variabili categoriali in R 2. Identificazione e Rimozione Outlier Visualizzazione outlier
neonati %>%
select_if(is.numeric) %>%
gather(key = "Variabile", value = "Valore") %>%
ggplot(aes(x = Variabile, y = Valore)) +
geom_boxplot(fill = "lightblue") +
labs(title = "Identificazione di outlier nelle variabili numeriche") +
theme_minimal()
Una volta identificati gli outlier, si possono decidere azioni
successive, come la rimozione o la trasformazione dei dati. Rimozione
degli outlier per la variabile Peso
neonati <- neonati %>%
filter(Peso > quantile(Peso, 0.25) - 1.5 * IQR(Peso) &
Peso < quantile(Peso, 0.75) + 1.5 * IQR(Peso))Eliminare i valori anomali (outlier) dalla variabile Peso basandosi sulla regola dell’IQR (Interquartile Range). 3. Controllo di Normalità e Trasformazioni
shapiro.test(neonati$Peso) # Test Shapiro-Wilk##
## Shapiro-Wilk normality test
##
## data: neonati$Peso
## W = 0.99826, p-value = 0.01038qqnorm(neonati$Peso) # Q-Q plot
qqline(neonati$Peso, col = "red")
Questo test verifica se la distribuzione della variabile Peso segue una
distribuzione normale. Trasformazione logaritmica di Peso
neonati <- neonati %>%
mutate(Peso_log = log(Peso + 1))Applicare una trasformazione logaritmica alla variabile Peso per affrontare problemi di non normalità , ridurre l’influenza degli outlier e stabilizzare la varianza 4. Feature Engineering Creazione di variabili derivate
neonati <- neonati %>%
mutate(Annimadre_sq = Annimadre^2,
Gestazione_sq = Gestazione^2,
Annimadre_fasce = cut(Annimadre,
breaks = c(0, 20, 30, 40, Inf),
labels = c("Teen", "Young", "Middle-aged", "Older")))Standardizzazione delle variabili numeriche
numerical_vars <- c("Annimadre", "Gestazione", "Peso", "Ngravidanze")
neonati <- neonati %>%
mutate(across(all_of(numerical_vars), ~ scale(.) %>% as.vector()))Creazione di nuove interazioni
neonati <- neonati %>%
mutate(Interaction1 = Annimadre * Gestazione,
Interaction2 = Ngravidanze * Gestazione_sq)- Analisi delle Correlazioni e Distribuzioni Matrice di correlazione
cor_matrix <- neonati %>%
select_if(is.numeric) %>%
cor(use = "complete.obs")
corrplot(cor_matrix, method = "color", type = "upper", tl.cex = 0.8, tl.col = "black")
Istogrammi delle variabili
ggplot(neonati, aes(x = Annimadre)) +
geom_histogram(binwidth = 1, fill = "skyblue", color = "black") +
labs(title = "Distribuzione dell'età delle madri")
ggplot(neonati, aes(x = Peso)) +
geom_histogram(binwidth = 1, fill = "lightgreen", color = "black") +
labs(title = "Distribuzione del peso dei neonati", x = "Peso (grammi)", y = "Frequenza") +
theme_minimal()
ggplot(neonati, aes(x = Ngravidanze)) +
geom_histogram(binwidth = 1, fill = "lightcoral", color = "black") +
labs(title = "Distribuzione del numero di gravidanze")
Boxplot del fumo materno rispetto al peso del neonato
ggplot(neonati, aes(x = Fumatrici, y = Peso)) +
geom_boxplot(fill = "lightblue", color = "black") +
labs(title = "Relazione tra fumo materno e peso del neonato", x = "Fumatrici", y = "Peso del neonato (grammi)")
Scatter plot dell’età della madre rispetto al peso del neonato
ggplot(neonati, aes(x = Annimadre, y = Peso)) +
geom_point(color = "darkblue") +
labs(title = "Relazione tra età della madre e peso del neonato", x = "Età della madre", y = "Peso del neonato (grammi)")
6. Modellazione e Validazione Impostazioni per la validazione
incrociata
train_control <- trainControl(method = "repeatedcv", number = 10, repeats = 3)Modello di regressione lineare con validazione incrociata
modello_cv <- train(Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione,
data = neonati, method = "lm", trControl = train_control)
print(modello_cv)## Linear Regression
##
## 2431 samples
## 4 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 2187, 2188, 2188, 2188, 2188, 2189, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.8944999 0.2023882 0.7194992
##
## Tuning parameter 'intercept' was held constant at a value of TRUECalcolo delle metriche di valutazione
predicted_values <- predict(modello_cv, newdata = neonati)
actual_values <- neonati$Peso
MAE <- mean(abs(predicted_values - actual_values))
MSE <- mean((predicted_values - actual_values) ^ 2)
RMSE <- sqrt(MSE)
RSS <- sum((predicted_values - actual_values) ^ 2)
TSS <- sum((actual_values - mean(actual_values)) ^ 2)
R_squared <- 1 - (RSS / TSS)
cat("MAE:", MAE, "\n")## MAE: 0.7182185cat("MSE:", MSE, "\n")## MSE: 0.7977624cat("R-squared:", R_squared, "\n")## R-squared: 0.2019093- Selezione Stepwise del Modello
modello_base <- lm(Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione, data = neonati)
modello_step <- stepAIC(modello_base, direction = "both", trace = FALSE)
summary(modello_step)##
## Call:
## lm(formula = Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9762 -0.6165 -0.0391 0.6217 3.7935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.01118 0.01853 0.604 0.5461
## Annimadre 0.04967 0.01969 2.522 0.0117 *
## Ngravidanze 0.04763 0.01968 2.421 0.0156 *
## Fumatrici1 -0.26656 0.09062 -2.941 0.0033 **
## Gestazione 0.45143 0.01831 24.656 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8941 on 2426 degrees of freedom
## Multiple R-squared: 0.2019, Adjusted R-squared: 0.2006
## F-statistic: 153.4 on 4 and 2426 DF, p-value: < 2.2e-16- Modelli Avanzati Modello ad Alberi di Decisione
modello_tree <- rpart(Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione, data = neonati, method = "anova")
pred_tree <- predict(modello_tree, newdata = neonati)
MAE_tree <- mean(abs(pred_tree - neonati$Peso))Modello Random Forest
modello_rf <- randomForest(Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione, data = neonati, ntree = 500, mtry = 2)
pred_rf <- predict(modello_rf, newdata = neonati)
MAE_rf <- mean(abs(pred_rf - neonati$Peso))Modello XGBoost
dati_xgb <- model.matrix(Peso ~ Annimadre + Ngravidanze + Fumatrici + Gestazione - 1, data = neonati)
peso_xgb <- neonati$Peso
dtrain <- xgb.DMatrix(data = dati_xgb, label = peso_xgb)
modello_xgb <- xgboost(data = dtrain, params = list(objective = "reg:squarederror", eta = 0.1, max_depth = 3), nrounds = 100, verbose = FALSE)
pred_xgb <- predict(modello_xgb, newdata = dati_xgb)
MAE_xgb <- mean(abs(pred_xgb - peso_xgb))Confronto tra modelli
comparison <- data.frame(Modello = c("Albero di Decisione", "Random Forest", "XGBoost"),
MAE = c(MAE_tree, MAE_rf, MAE_xgb))
ggplot(comparison, aes(x = Modello, y = MAE, fill = Modello)) +
geom_bar(stat = "identity") +
geom_text(aes(label = round(MAE, 2)), vjust = -0.5) +
labs(title = "Confronto delle Performance dei Modelli", x = "Modello", y = "MAE") +
theme_minimal()
9. Modello con Interazioni ed Effetti Non Lineari Modello di regressione
lineare con interazioni
modello_interaction <- lm(Peso ~ Annimadre * Fumatrici + Ngravidanze * Gestazione, data = neonati)
summary(modello_interaction)##
## Call:
## lm(formula = Peso ~ Annimadre * Fumatrici + Ngravidanze * Gestazione,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9841 -0.6160 -0.0341 0.6219 3.7353
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.010227 0.018601 0.550 0.58249
## Annimadre 0.051933 0.020096 2.584 0.00982 **
## Fumatrici1 -0.263613 0.090709 -2.906 0.00369 **
## Ngravidanze 0.047360 0.019687 2.406 0.01622 *
## Gestazione 0.452222 0.018361 24.630 < 2e-16 ***
## Annimadre:Fumatrici1 -0.056853 0.087327 -0.651 0.51509
## Ngravidanze:Gestazione -0.009321 0.015416 -0.605 0.54550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8943 on 2424 degrees of freedom
## Multiple R-squared: 0.2022, Adjusted R-squared: 0.2002
## F-statistic: 102.4 on 6 and 2424 DF, p-value: < 2.2e-16Calcolo del VIF per verificare la collinearitÃ
vif_interaction <- vif(modello_interaction)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vifprint(vif_interaction)## Annimadre Fumatrici Ngravidanze
## 1.226933 1.005322 1.177521
## Gestazione Annimadre:Fumatrici Ngravidanze:Gestazione
## 1.024220 1.048758 1.008372Modello aggiornato con interazioni e termini quadratici
modello_interaction_non_lineare <- lm(Peso ~ Annimadre * Fumatrici + Ngravidanze * Gestazione +
Annimadre_sq + Gestazione_sq, data = neonati)
summary(modello_interaction_non_lineare)##
## Call:
## lm(formula = Peso ~ Annimadre * Fumatrici + Ngravidanze * Gestazione +
## Annimadre_sq + Gestazione_sq, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9301 -0.6061 -0.0338 0.6206 3.8715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.2318384 7.4500828 3.521 0.000438 ***
## Annimadre 0.1906301 0.1268635 1.503 0.133063
## Fumatrici1 -0.2680175 0.0905112 -2.961 0.003095 **
## Ngravidanze 0.0504973 0.0197180 2.561 0.010498 *
## Gestazione 2.3933788 0.5601097 4.273 2e-05 ***
## Annimadre_sq -0.0004788 0.0004266 -1.122 0.261829
## Gestazione_sq -0.0168447 0.0048551 -3.469 0.000531 ***
## Annimadre:Fumatrici1 -0.0593220 0.0872208 -0.680 0.496483
## Ngravidanze:Gestazione -0.0099435 0.0153875 -0.646 0.518207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8922 on 2422 degrees of freedom
## Multiple R-squared: 0.2066, Adjusted R-squared: 0.2039
## F-statistic: 78.81 on 8 and 2422 DF, p-value: < 2.2e-16Calcolo del VIF per il modello aggiornato
vif_non_lineare <- vif(modello_interaction_non_lineare)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vifprint(vif_non_lineare)## Annimadre Fumatrici Ngravidanze
## 49.127892 1.005629 1.186809
## Gestazione Annimadre_sq Gestazione_sq
## 957.636818 49.644788 957.681918
## Annimadre:Fumatrici Ngravidanze:Gestazione
## 1.051118 1.009318Confronto tra modelli
AIC(modello_interaction, modello_interaction_non_lineare)## df AIC
## modello_interaction 8 6364.811
## modello_interaction_non_lineare 10 6355.431BIC(modello_interaction, modello_interaction_non_lineare)## df BIC
## modello_interaction 8 6411.180
## modello_interaction_non_lineare 10 6413.392Modello GAM
gam_model <- gam(Peso ~ s(Annimadre) + s(Gestazione) + Fumatrici + Ngravidanze, data = neonati)
summary(gam_model)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Peso ~ s(Annimadre) + s(Gestazione) + Fumatrici + Ngravidanze
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.01137 0.01848 0.615 0.53857
## Fumatrici1 -0.27094 0.09041 -2.997 0.00276 **
## Ngravidanze 0.04911 0.01965 2.500 0.01250 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Annimadre) 1.000 1.00 6.073 0.0138 *
## s(Gestazione) 2.655 3.37 184.528 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.205 Deviance explained = 20.6%
## GCV = 0.79764 Scale est. = 0.79546 n = 2431plot(gam_model, pages = 1, se = TRUE, shade = TRUE)
Confronto AIC/BIC con GAM
AIC(modello_interaction_non_lineare, gam_model)## df AIC
## modello_interaction_non_lineare 10.000000 6355.431
## gam_model 7.654729 6351.230BIC(modello_interaction_non_lineare, gam_model)## df BIC
## modello_interaction_non_lineare 10.000000 6413.392
## gam_model 7.654729 6395.597- Visualizzazione degli Effetti Non Lineari
ggplot(neonati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ x + I(x ^ 2), color = "blue") +
labs(title = "Effetto di Gestazione e Gestazione^2 sul Peso del Neonato",
x = "Settimane di Gestazione", y = "Peso del Neonato (g)") +
theme_minimal()
ggplot(neonati, aes(x = Annimadre, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", formula = y ~ x + I(x ^ 2), color = "green") +
labs(title = "Effetto di Annimadre e Annimadre^2 sul Peso del Neonato",
x = "Età della Madre", y = "Peso del Neonato (g)") +
theme_minimal()
Analisi dei residui
ggplot(data.frame(residuals = residuals(modello_step), fitted = fitted(modello_step)), aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Grafico dei Residui", x = "Valori Predetti", y = "Residui")
qqnorm(residuals(modello_step))
qqline(residuals(modello_step), col = "red")
# Visualizzazione dei coefficienti del modello
coefficients_df <- tidy(modello_step)
ggplot(coefficients_df, aes(x = term, y = estimate)) +
geom_bar(stat = "identity", fill = "lightblue", color = "black") +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error), width = 0.2) +
labs(title = "Effetto delle Variabili sul Peso del Neonato", x = "Variabili", y = "Stima dei Coefficienti") +
theme_minimal()
# Analisi di SensibilitÃ
gestazione_vals <- seq(min(neonati$Gestazione, na.rm = TRUE), max(neonati$Gestazione, na.rm = TRUE), by = 1)
fumatrici_vals <- levels(neonati$Fumatrici)
sens_data <- expand.grid(Gestazione = gestazione_vals, Fumatrici = fumatrici_vals,
Annimadre = mean(neonati$Annimadre, na.rm = TRUE),
Ngravidanze = mean(neonati$Ngravidanze, na.rm = TRUE))
sens_data$Predicted_Peso <- predict(modello_step, newdata = sens_data)
ggplot(sens_data, aes(x = Gestazione, y = Predicted_Peso, color = Fumatrici)) +
geom_line() +
labs(title = "Impatto del Fumo e della Gestazione sul Peso Previsto del Neonato",
x = "Settimane di Gestazione", y = "Peso del Neonato (g)") +
theme_minimal()