dati<-read.csv(“C:/Users/angel/Desktop/Data Science/R statistica inferenziale/neonati.csv”)
summary(dati)
#analisi variaibili
library(psych) describe(dati) # Mostra statistiche estese per tutte le variabili
colSums(is.na(dati))
par(mfrow = c(3, 2)) #finestra 3x2 hist(dati\(Peso, main = "Distribuzione Peso Neonatale", xlab = "Peso (grammi)", col = "lightblue") boxplot(dati\)Peso, main = “Boxplot Peso Neonatale”, col = “lightblue”, horizontal = TRUE) hist(dati\(Gestazione, main = "Distribuzione Settimane di Gestazione", xlab = "Settimane", col = "lightgreen") boxplot(dati\)Gestazione, main = “Boxplot Settimane di Gestazione”, col = “lightgreen”, horizontal = TRUE) hist(dati\(Lunghezza, main = "Distribuzione Lunghezza Neonati", xlab = "Lunghezza (mm)", col = "lightpink") boxplot(dati\)Lunghezza, main = “Boxplot Lunghezza Neonati”, col = “lightpink”, horizontal = TRUE)
par(mfrow = c(1, 1))
peso_outlier <- boxplot.stats(dati\(Peso)\)out peso_outlier # valori anomali del peso
gestazione_outlier <- boxplot.stats(dati\(Gestazione)\)out gestazione_outlier
library(ggplot2)
ggplot(dati, aes(x = Sesso, fill = Sesso)) + geom_bar() + labs(title = “Distribuzione per Sesso”, x = “Sesso”, y = “Frequenza”) + scale_fill_manual(values = c(“lightblue”, “pink”)) + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, size = 16, face = “bold”), axis.title = element_text(size = 12), axis.text = element_text(size = 10)) + geom_text(stat = ‘count’, aes(label = ..count..), vjust = -0.5, size = 5)
ggplot(dati, aes(x = Tipo.parto, fill = Tipo.parto)) + geom_bar() + labs(title = “Distribuzione per Tipo di Parto”, x = “Tipo di Parto”, y = “Frequenza”) + scale_fill_brewer(palette = “Set2”) + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, size = 16, face = “bold”), axis.title = element_text(size = 12), axis.text = element_text(size = 10)) + geom_text(stat = ‘count’, aes(label = ..count..), vjust = -0.5, size = 5)
ggplot(dati, aes(x = Ospedale, fill = Ospedale)) + geom_bar() + labs(title = “Distribuzione per Ospedale”, x = “Ospedale”, y = “Frequenza”) + scale_fill_viridis_d() + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, size = 16, face = “bold”), axis.title = element_text(size = 12), axis.text = element_text(size = 10)) + geom_text(stat = ‘count’, aes(label = ..count..), vjust = -0.5, size = 5)
ggplot(dati, aes(x = Fumatrici, fill = Fumatrici)) + geom_bar() + labs(title = “Distribuzione per Fumo Materno”, x = “Fumatrici”, y = “Frequenza”) + scale_fill_manual(values = c(“lightgreen”, “lightcoral”)) + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, size = 16, face = “bold”), axis.title = element_text(size = 12), axis.text = element_text(size = 10)) + geom_text(stat = ‘count’, aes(label = ..count..), vjust = -0.5, size = 5)
ggplot(dati, aes(x = Tipo.parto, fill = Ospedale)) + geom_bar(position = “dodge”) + labs(title = “Frequenze per Tipo di Parto per Ospedale”, x = “Tipo di Parto”, y = “Frequenza”) + scale_fill_manual(name = “Ospedale”, values = c(“lightblue”, “lightgreen”, “lightcoral”)) + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, size = 16, face = “bold”), axis.title = element_text(size = 12), axis.text = element_text(size = 10)) + geom_text(stat = ‘count’, aes(label = ..count..), position = position_dodge(0.9), vjust = -0.5, size = 5)
plot(dati\(Gestazione, dati\)Peso, main = “Peso vs Settimane di Gestazione”, xlab = “Settimane di Gestazione”, ylab = “Peso (grammi)”, col = “darkblue”, pch = 19)
boxplot(dati\(Peso ~ dati\)Fumatrici, main = “Peso Neonatale e Fumo Materno”, xlab = “Fumo Materno (0 = No, 1 = Sì)”, ylab = “Peso (grammi)”, col = c(“lightblue”, “lightcoral”))
tabella_parti <- table(dati\(Ospedale, dati\)Tipo.parto)
barplot( tabella_parti, beside = TRUE, main = “Distribution of Delivery Types by Hospital”, # Titolo del grafico col = c(“lightgreen”, “lightblue”, “lightgrey”), # Colori per i tipi di parto legend.text = c(“Hospital A”, “Hospital B”, “Hospital C”), # Modifica leggenda args.legend = list(title = “Hospitals”, x = 3.5, y = 650, bty = “n”), xlab = “Type of Delivery”,
ylab = “Frequency”, ylim = c(0, 650)
)
bar_positions <- barplot(tabella_parti, beside = TRUE, plot = FALSE) # posizioni delle barre text( x = as.vector(bar_positions), y = as.vector(tabella_parti), labels = as.vector(tabella_parti), pos = 3, cex = 0.8, col = “black” # Etichette personalizzate )
#valutazioni sul campione in esame ——
chisq_test <- chisq.test(tabella_parti) chisq_test
#NON Si può rifiutare l’ipotesi nulla di ugualianza. il tipo di parto non è legato all’ospedale
t_test_peso <- t.test(dati$Peso, mu = 3300, alternative = “two.sided”) t_test_peso # Risultato del test #la media del campione non è significativamente diversa da quella della popolazione
t_test_lunghezza <- t.test(dati$Lunghezza, mu = 500, alternative = “two.sided”) t_test_lunghezza
#la media del campione è significativamente diversa da quella della popolazione per la lunghezza anche se quasi sempre entro circa 7mm poichè l’estremo inferiore è 493.6
#test condizionati dal Sesso t_test_peso_sesso <- t.test(Peso ~ Sesso, data = dati) t_test_peso_sesso
t_test_lunghezza_sesso <- t.test(Lunghezza ~ Sesso, data = dati) t_test_lunghezza_sesso
t_test_cranio_sesso <- t.test(Cranio ~ Sesso, data = dati) t_test_cranio_sesso # sono tutti significativi indicando differenze tra i due generi per peso lunghezza e cranio, come previsto dalla letteratura
shapiro.test(Peso) # Poiché il p-value è inferiore a 0.05, rifiutiamo l’ipotesi nulla di normalità #i dati relativi al Peso non seguono una distribuzione normale
dati\(Ospedale <- as.factor(dati\)Ospedale) dati\(Tipo.parto <- as.factor(dati\)Tipo.parto) dati\(Sesso <- as.factor(dati\)Sesso)
str(dati)
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) } #correlazioni pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)
#Tenendo presente lo scopo reale del modello di previsione farei delle considerazioni preliminari: #sembra che con il peso incidano le variabili GESTAZIONE LUNGHEZZA, CRANIO, SESSO #studi in letteratura indica che il fumo può incidere negativamente sul peso, lunghezza e parto prematura a #causa del minore apporto di ossigeno (fonte: salute.gov) #inoltre è noto in letteratura come una madre fumatrice possa innescare o favorire altre patologie al bambino. #Per queste considerazioni anche se in questo specifico campione non risulta incisivo l’apporto del fumo terrei #lo stesso tale variabile nel modello in quanto mantere la variabile può dare maggiore generalizzabilità al modello #Oltre il campione specifico, altra motivazione della non incidenza del fumo potrebbe essere che la variabile è stata caratterizzata solo #in due livelli. Valutare i mesi o anni di fumo precedenti oppure la quanitità di fumo potrebbe essere utile #per future ricerche per maggiore prevedibilità del peso # Visto che lo scopo è utilizzare al meglio le TIV (terapie intensive) e le risorse umane terrei gli anni della madre
#visivamente la relazione PESO - LUNGHEZZA E PESO-CRANIO sembrano lineari, la relazione PESO-GESTAZIONE appare curvilinea # per questo proverei a valutare tale relazione per una maggiore prevedibilità (in letteratura è noto che nascere con molte #settimane di anticipo compromette in modo maggiore rispetto a differenze più vicine al completamento della gestazione)
t.test(Peso~Sesso) t.test(Peso~Tipo.parto)
#il genere come prevedibile incide sul peso di circa 150g mentre il tipo di parto no è significativo
#anova per gli ospedali anova_result <- aov(Peso ~ Ospedale, data = dati) summary(anova_result)
boxplot(Peso ~ Ospedale, data = dati, pch=20, main = “Peso per Ospedale”, xlab = “Ospedale”, ylab = “Peso”, col = “lightblue”) #come si evince dal p value e dai box plot non ci sono differenze significative tra peso e Ospedale
modello_peso <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dati) summary(modello_peso)
library(MASS) modello_aic <- stepAIC(modello_peso, direction = “both”) #calcolo in automatico il modello migliore partendo da quello completo summary(modello_aic) #con il metodo AIC si ottiene il modello composto N.gravidanze Gestazione Lunghezza Cranio Tipo.partoNat Ospedale Sesso come variabili esplicative del peso #R quadro aggiustato risualta 0.7265 #il fumo è stato eliminato poichè pare non dare contributo a PESO
n <- nrow(dati) modello_bic <- stepAIC(modello_peso, direction = “both”, k = log(n)) summary(modello_bic) #con il metodo BIC si ottiene il modello composto N.gravidanze Gestazione Lunghezza Cranio Sesso come variabili esplicative del peso #R quadro aggiustato risualta 0.7265 identico al metodo AIC anche come p value #dalle considerazioni fatte precedentemente questo metodo sembra dare il risultato #più simile a ciò che ci si aspetta dalla letteratura , tranne per il fumo che è stato nuovamente eliminato. #Forse come spiegato sopra per il campione o per il metodo dicotomico di valutarlo
#Riprendo tale modello indagando possibili interazioni modello_interazioni <- lm(Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio)^2, data = dati) summary(modello_interazioni)
#in combinazione con la letteratura, delle interazioni significative potremmo # considerare Anni.madre:Gestazione Anni.madre:Lunghezza Fumatrici:Gestazione Fumatrici:Lunghezza # però è possibile notare che seppur significativi i contributi sono in ordine: +2.2 e 0.24 grammi per cui trascurabili #l’essere fumatrice modifica significativamente l’effetto della gestazione sul peso del bambino #e il coeff beta non è trascurabile per una possibile previsione essendo di 53g in più
modello_nonlineare <- lm(Peso ~ poly(Gestazione, 2) + Lunghezza + Anni.madre + Fumatrici, data = dati) summary(modello_nonlineare) #Il termine quadratico dela Gestazione non è significativo e potrebbe essere rimosso per semplificare il modello.
modello_peso_ridotto <- lm(Peso ~ N.gravidanze + Gestazione +Lunghezza + Cranio + Sesso, data = dati) summary(modello_peso_ridotto)
modello_peso_int <- lm(Peso ~ N.gravidanze + Gestazione +Lunghezza + Cranio + Sesso + Fumatrici:Gestazione, data = dati) summary(modello_peso_int) # inserito nel modello ridotto Gestazione:Fumatrici risulta non significativo per cui eliminabile
AIC(modello_peso, modello_peso_ridotto, modello_peso_int) BIC(modello_peso, modello_peso_ridotto, modello_peso_int) #il metodo BIC che predilige modelli più semplici conferma che che il modello migliore è modello_peso_ridotto
#c’è da decidere se includere nel modello di previsione il fumo (e/o sue interazioni) o meno per dare maggiore generalizzabilità anche se questi dati non #suggeriscono che vada inserita tale variabile. Questa considerazione dovrebbe essere fatta con un confronto #adeguato con esperti del settore per valutazioni cliniche. Per quanto esposto nei dati non dovrebbe essere inserito
#scelto il modello valutiamo le altre info sui residui #analisi residui par(mfrow=c(2,2)) plot(modello_peso_ridotto) # tutti vanno bene. viene visualizzata la distanza di cook e sembra indicare che il valore 1551 peggiori il modello # lo stesso valore come è evidente dai residui e dalla distribuzione dei quartili appare anomalo #verificiamolo con i test lmtest::bptest(modello_peso_ridotto) #Poiché il p-value è < 0.01, rifiutiamo l’ipotesi nulla di omoschedasticità lmtest::dwtest(modello_peso_ridotto) #I residui non presentano una correlazione significativa tra di loro e #soddisfano l’assunzione di indipendenza. Non possiamo rifiutare l’ipotesi nulla di autocorrelazione dei residui shapiro.test(modello_peso_ridotto$residuals) #rifiutiamo l’ipotesi nulla di normalità dei residui. #I residui non seguono una distribuzione normale, il che potrebbe indicare che il modello lineare non è #perfettamente appropriato (ed esempio effetti non lineari) o che nei dati vi siano outlier. plot(density(residuals(modello_peso_ridotto)))
#mi aspetto che l’osservazione 1551 sia oltre la soglia dello 0.5 dell’indice di Cool #leverage :quanto un’osservazione è distante dallo spazio medio delle variabili indipendenti. #Un’osservazione con alta leverage può influenzare notevolmente i coefficienti del modello lev<-hatvalues(modello_peso_ridotto) plot(lev) p<-sum(lev) n<-length(lev) soglia=2*p/n abline(h=soglia,col=2) #la soglia individua molti dati oltre però leverage alto non implica necessariamente che l’osservazione sia problematica lev[lev>soglia]
#outliers #Gli outlier nei residui standardizzati mostrano osservazioni che deviano significativamente dal modello plot(rstudent(modello_peso_ridotto)) abline(h=c(-2,2)) car::outlierTest(modello_peso_ridotto) # anche in questo caso emergono molte osservazioni sospette secondo il test, vediamo se è vero con la distanza di cook
#distanza di cook #combina gli effetti di leveare e outlier per valutare se il dato è anomalo rispetto al modello cook<-cooks.distance(modello_peso_ridotto) plot(cook,ylim = c(0,1)) abline(h = 0.5, col = “red”) # solo una osservazione è oltre 0.5 osservazioni_influenti <- which(cook > 0.5) print(osservazioni_influenti) # Mostra gli indici delle osservazioni influenti
#è lecito pensare di elimianre l’osservazione 1551 dati_modificati <- dati[-osservazioni_influenti,] nrow(dati) - nrow(dati_modificati) # verifico che il db abbia giusto numero di dati
modello_peso_ridotto_mod <- lm(Peso ~ N.gravidanze + Gestazione +Lunghezza + Cranio + Sesso, data = dati_modificati) summary(modello_peso_ridotto_mod)
AIC(modello_peso_ridotto, modello_peso_ridotto_mod) BIC(modello_peso_ridotto, modello_peso_ridotto_mod) #entrambi gli indici migliorano per il modello dopo la modifica
#ripeto i test par(mfrow=c(2,2)) plot(modello_peso_ridotto_mod) # tutti vanno bene. viene visualizzata la distanza di cook e sembra indicare che il valore 1551 peggiori il modello # lo stesso valore come è evidente dai residui e dalla distribuzione dei quartili appare anomalo #verificiamolo con i test lmtest::bptest(modello_peso_ridotto_mod) #Poiché il p-value è < 0.05, rifiutiamo l’ipotesi nulla di omoschedasticità ma non la possiamo rifiutare allo 0.01 lmtest::dwtest(modello_peso_ridotto_mod) #continuo a poter rifiutare l’autocorrelazione shapiro.test(modello_peso_ridotto_mod$residuals) #rifiutiamo ancora l’ipotesi nulla di normalità dei residui.
pred <- predict(modello_peso_ridotto_mod, newdata = dati_modificati) actual <- dati_modificati$Peso
R2 <- round(cor(pred, actual)^2,3) RMSE <- round(sqrt(mean((pred - actual)^2)),2) cat(“R-squared:”, R2, “:”, RMSE) #il modello spiega il 73.7% della varianza dei dati
nuovo_neonato <- data.frame( Anni.madre = 30, #informazione superflua al modello N.gravidanze = 3, Fumatrici = 0, #informazione superflua al modello Gestazione = 39, Lunghezza = 500, Cranio = 340, Tipo.parto = “Naturale”, #informazione superflua al modello Ospedale = “Ospedale 1”, #informazione superflua al modello Sesso = “F” )
peso_previsto <- predict(modello_peso_ridotto_mod, newdata = nuovo_neonato) cat(“Peso previsto del neonato con i dati inseriti:”, round(peso_previsto,1), “grammi”)
coeff <- summary(modello_peso_ridotto_mod)$coefficients; coeff round(coeff[“N.gravidanze”, ],2) # Impatto del numero di gravidanze round(coeff[“Gestazione”, ],2) # Impatto delle settimane di gestazione round(coeff[“Lunghezza”, ],2) round(coeff[“Cranio”, ],2) round(coeff[“SessoM”, ],2)
cat(“Ogni settimana aggiuntiva di gestazione aumenta in media il peso del neonato di:”, round(coeff[“Gestazione”, “Estimate”], 2), “grammi.”)
cat(“Una lunghezza unitaria aggiuntiva aumenta in media il peso del neonato di:”, round(coeff[“Lunghezza”, “Estimate”], 2), “grammi.”)
cat(“Un aumento di circonferenza del Cranio aggiuntiva aumenta in media il peso del neonato di:”, round(coeff[“Cranio”, “Estimate”], 2), “grammi.”)
cat(“Ogni Gravidanza aggiuntiva modifica in media il peso del neonato di:”, round(coeff[“N.gravidanze”, “Estimate”], 2), “grammi.”)
if (“SessoM” %in% rownames(coeff)) { cat(“Il genere maschile modifica in media il peso del neonato di:”, abs(round(coeff[“SessoM”, “Estimate”], 2)), “grammi.”) }
#— Visualizzazione
ggplot(dati, aes(x = Gestazione, y = Peso, color = factor(Sesso))) + geom_point(position = position_jitter(width = 0.2, height = 0)) + geom_smooth(method = “lm”, color = “blue”) + labs(title = “Relazione tra Settimane di Gestazione e Peso Neonatale”, x = “Settimane di Gestazione”, y = “Peso Neonatale (g)”, color = “Sesso”) + theme_minimal() #si nota dalla zona grigia intorno alla retta che l’accuratezza si riduce quando le settimane di gestazione si riducono
ggplot(dati, aes(x = factor(Fumatrici, labels = c(“Non Fumatrici”, “Fumatrici”)), y = Peso, fill = factor(Fumatrici))) + geom_boxplot() + labs(title = “Impatto del Fumo Materno sul Peso Neonatale”, x = “Fumo Materno”, y = “Peso Neonatale (g)”) + scale_fill_manual(values = c(“green”, “blue”)) + theme_minimal()
ggplot(dati, aes(x = Ospedale, y = Peso, fill = Ospedale)) + geom_violin(alpha = 0.7) + labs(title = “Distribuzione del Peso Neonatale nei Tre Ospedali”, x = “Ospedale”, y = “Peso Neonatale (g)”) + theme_minimal()
ggplot(dati, aes(x = N.gravidanze, y = Gestazione, fill = Peso)) + geom_tile() + labs(title = “Interazione tra Numero di Gravidanze e Settimane di Gestazione”, x = “Numero di Gravidanze”, y = ““, fill =”Peso (g)“) + scale_fill_gradient(low =”yellow”, high = “red”) + theme_minimal() #da un analisi verticale dei colori appare anche qui evidente come non c’è una correlazione evidente tra le variabili
ggplot(data=dati)+ geom_point(aes(x=Cranio, y=Peso, col=Sesso),position = “jitter”)+ geom_smooth(aes(x=Cranio, y=Peso, col=Sesso),se=F,method = “lm”)+ labs(title = “Relazione tra circonferenza cranio e peso (g)”, x = “Circonferenza cranio”, y = ““, fill =”Peso (g)“) #si nota come la retta del genere maschile è sempre maggiore rispetto a quella femminile nel 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=F,method = “lm”)+ labs(title = “Relazione tra Lunghezza (cm) e peso (g)”, x = “Lunghezza (cm)”, y = ““, fill =”Peso (g)“)