Report Statistico
Sara Greco
library(ggplot2)
library(corrplot)## Warning: il pacchetto 'corrplot' è stato creato con R versione 4.5.2## corrplot 0.95 loadedlibrary(dplyr)## Warning: il pacchetto 'dplyr' è stato creato con R versione 4.5.1##
## 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(cluster)
library(factoextra)## Warning: il pacchetto 'factoextra' è stato creato con R versione 4.5.2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.5.2## Caricamento del pacchetto richiesto: carData##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recodelibrary(lmtest)## Warning: il pacchetto 'lmtest' è stato creato con R versione 4.5.2## Caricamento del pacchetto richiesto: zoo## Warning: il pacchetto 'zoo' è stato creato con R versione 4.5.2##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericsetwd("C:/Users/Sara/Desktop/Report")
df = read.csv("Boston.csv")
dim(df)## [1] 506 14sapply(df, class)## CRIM ZN INDUS CHAS NOX RM AGE DIS
## "numeric" "numeric" "numeric" "integer" "numeric" "numeric" "numeric" "numeric"
## RAD TAX PTRATIO B LSTAT Price
## "integer" "integer" "numeric" "numeric" "numeric" "numeric"EDA
colSums((is.na(df))) #verifico se ci sono valori nulli## CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX
## 0 0 0 0 0 0 0 0 0 0
## PTRATIO B LSTAT Price
## 0 0 0 0sum(duplicated((df))) #verifico se ci sono valori duplicati## [1] 0Non sono presenti ne valori duplicati, ne varoli nulli
summary(df)## CRIM ZN INDUS CHAS
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## NOX RM AGE DIS
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## RAD TAX PTRATIO B
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## LSTAT Price
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00ANALISI UNIVARIATA
df$CHAS <- factor(df$CHAS, levels = c(0,1), labels = c("No", "Si"))
ggplot(df, aes(x=CHAS)) +
geom_bar(fill = "skyblue", color = "black") +
labs(title = "Distribuzione della variabile CHAS",
x = "Confina con il fiume Charles",
y = "Numero di abitazioni") +
theme_minimal()
La maggior parte delle abitazioni non confina con il fiume Charles (categoria = “No”) con oltre 400 casi. Solo un piccolo numero di abitazioni confina col fiume (categoria = “Si”) circa una quarantina di casi
par(mfrow = c(1, 4))
boxplot(df$CRIM, main='CRIM',col='Sky Blue')
boxplot(df$ZN, main='ZN',col='Sky Blue')
boxplot(df$RM, main='RM',col='Sky Blue')
boxplot(df$B, main='B',col='Sky Blue')
Dai boxplot emerge che la variabile CRIM presenta numerosi outlier verso l’alto, con valori che superano ampiamente 40 e si avvicinano a 90, mentre la maggior parte dei dati è concentrata su valori bassi, indicando una distribuzione fortemente asimmetrica.
La variabile ZN mostra anch’essa diversi outlier elevati, con valori che superano 60 fino quasi a 100, mentre la maggioranza dei valori è molto bassa o nulla, con numerosi zeri.
Per RM, la distribuzione appare più concentrata attorno alla media, circa 6–7 stanze, con outlier presenti ma meno estremi rispetto a CRIM e ZN, suggerendo una variabilità moderata.
Infine, B mostra valori molto alti per la maggior parte dei dati, con outlier verso il basso, e una distribuzione sbilanciata verso il valore massimo, che supera 300 e arriva fino a quasi 400.
par(mfrow = c(3, 3))
hist(df$INDUS, main = "INDUS", col = "lightblue", breaks = 20)
hist(df$TAX, main = "TAX", col = "lightgreen", breaks = 20)
hist(df$PTRATIO, main = "PTRATIO", col = "pink", breaks = 20)
hist(df$NOX, main = "NOX", col = "lightyellow", breaks = 20)
hist(df$AGE, main = "AGE", col = "lightgray", breaks = 20)
hist(df$DIS, main = "DIS", col = "orange", breaks = 20)
hist(df$RAD, main = "RAD", col = "lightcoral", breaks = 20)
hist(df$LSTAT, main = "LSTAT", col = "lightcyan", breaks = 20)
hist(df$Price, main = "Price", col = "lightpink", breaks = 20)
Dall’analisi preliminare della matrice di istogrammi emergono alcune evidenze rilevanti sulla distribuzione delle variabili considerate.
AGE presenta una distribuzione fortemente asimmetrica verso destra, con un picco molto pronunciato in prossimità del valore massimo. DIS mostra anch’essa un’asimmetria positiva, sebbene più graduale e distribuita.
RAD e TAX mostrano picchi netti su pochi valori.
LSTAT e Price hanno distribuzioni più regolari e continue.
ANALISI BIVARIATA
ggplot(df, aes(x = CHAS, y = Price)) +
geom_boxplot(fill = "lightgreen") +
labs(title = "Distribuzione del prezzo in base a CHAS",
x = "Confina con fiume Charles",
y = "Prezzo") +
theme_minimal()
Le abitazioni che confinano con il fiume Charles tendono ad avere prezzi mediamente più alti rispetto a quelle che non confinano. La mediana del gruppo (“Si”), infatti, è più alta di quella del gruppo “No”. La distribuzione del prezzo per il gruppo “No” sembra più stretta ma con diversi outlier molto alti.
num_vars <- df %>% select_if(is.numeric)
cor_matrix <- cor(num_vars)
corrplot(cor_matrix, method = "color", type = "upper",
addCoef.col = "black",
number.cex = 0.7,
tl.cex = 0.8,
diag = FALSE)
Il prezzo (Price) è fortemente correlato con alcune variabili numeriche
# Seleziona le variabili di interesse
vars <- df[, c("NOX", "RM", "LSTAT", "PTRATIO", "Price", "INDUS", "TAX", "AGE", "DIS")]
# Crea scatterplot
pairs(vars,
pch = 19,
col = rgb(0, 0, 0.6, 0.9),
cex = 0.6,
oma = c(4, 4, 4, 4),
labels = c("NOX", "RM", "LSTAT", "PTRATIO", "TAX", "Price", "INDUS", "AGE", "DIS"))
L’analisi mostra che RM è positivamente correlata con Price, mentre LSTAT, PTRATIO e TAX presentano relazioni negative, più o meno marcate, indicando che il prezzo delle abitazioni aumenta con più stanze e diminuisce all’aumentare dello status basso, del rapporto alunni/insegnanti o della tassazione. Tra le variabili indipendenti, si osserva una correlazione negativa tra RM e LSTAT e una positiva tra NOX e TAX, mentre altre relazioni risultano meno evidenti.
CLUSTERING
# Lista delle variabili da escludere
dropList <- c('CHAS', 'RAD', 'CRIM', 'ZN', 'B')
# Seleziona le variabili numeriche tranne quelle in dropList
vars_num <- df[, !(names(df) %in% dropList)]
# Standardizza
vars_scaled <- scale(vars_num)wss <- numeric(10)
set.seed(123)
for(k in 1:10){
wss[k] <- kmeans(vars_scaled, centers = k, nstart = 25)$tot.withinss
}
plot(1:10, wss, type="b",
xlab = "Numero di cluster",
ylab = "WSS",
main = "Metodo del gomito")
Il metodo del gomito mostra che passare da 1 a 2 cluster riduce molto l’errore, mentre aggiungerne altri non porta grandi miglioramenti.
avg_sil <- numeric(10) # vettore per salvare la silhouette media
for(k in 2:10){ # parte da 2 perché per k=1 non ha senso
kmeans_k <- kmeans(vars_scaled, centers = k, nstart = 25)
sil <- silhouette(kmeans_k$cluster, dist(vars_scaled))
avg_sil[k] <- mean(sil[, 3]) # media valori silhouette
}
plot(2:10, avg_sil[2:10], type = "b",
xlab = "Numero di cluster", ylab = "silhouette",
main = "Metodo silhouette")
L’indice silhouette conferma che 2 cluster sono i più appropriati per descrivere i dati.
set.seed(123)
kmeans_result <- kmeans(vars_scaled, centers = 2, nstart = 25)
kmeans_result$size## [1] 228 278variabilita <- kmeans_result$betweenss / kmeans_result$totss * 100
variabilita## [1] 42.65044fviz_cluster(kmeans_result, data = vars_scaled,
ellipse.type = "convex",
palette = "jco",
ggtheme = theme_minimal(),
main = "Cluster",
geom = "point",
labelsize = 0)
La figura mostra la proiezione bidimensionale dei dati nelle prime due componenti principali (Dim1 e Dim2), che insieme spiegano circa il 73% della varianza totale.
Il cluster 1 (in blu) appare compatto e localizzato maggiormente nella parte destra del grafico lungo la prima componente (Dim1).
Il cluster 2 (in giallo) è più ampio e concentrato verso la parte sinistra del grafico su Dim1.
La distinzione fra i due gruppi risulta ben definita e coerente con la precedente analisi quantitativa (metodo del gomito e indice silhouette), che aveva indicato k=2 come soluzione ottimale.
PCA
pca_result <- prcomp(vars_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.2507 1.2220 0.89267 0.69986 0.60997 0.47420 0.46144
## Proportion of Variance 0.5628 0.1659 0.08854 0.05442 0.04134 0.02499 0.02366
## Cumulative Proportion 0.5628 0.7288 0.81730 0.87172 0.91306 0.93804 0.96170
## PC8 PC9
## Standard deviation 0.44810 0.37931
## Proportion of Variance 0.02231 0.01599
## Cumulative Proportion 0.98401 1.00000fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 60),
main = "Varianza spiegata")
La prima componente spiega il 56,3% della varianza, la seconda aggiunge il 16,6%, coprendo insieme oltre il 70%. Le componenti successive contribuiscono ciascuna per meno del 10%.
fviz_pca_var(pca_result,
col.var = "contrib",
gradient.cols = c("blue", "yellow", "red"),
repel = TRUE,
title = "Mappa delle variabili - PCA")
La mappa delle variabili mostra chiaramente il contributo e la direzione delle diverse variabili originali nelle prime due componenti della PCA, che complessivamente spiegano circa il 73% della varianza totale.
La prima componente (Dim1), che spiega il 56,3% della varianza, è dominata da variabili come NOX, AGE, INDUS e Price, suggerendo che questo asse rappresenta fattori legati principalmente a caratteristiche ambientali e socioeconomiche.
La seconda componente (Dim2), che spiega il 16,6%, è influenzata maggiormente da variabili come PTRATIO, TAX e LSTAT, che potrebbero rappresentare dimensioni alternative del fenomeno in studio, come qualità della vita o pressione fiscale.
REGRESSIONE LINEARE MULTIPLA
df$RM_c <- scale(df$RM, center=T, scale=F)
df$LSTAT_c <- scale(df$LSTAT, center=T, scale=F)
model1 <- lm(Price ~ RM_c + LSTAT_c + RM_c:LSTAT_c, data=df)
summary(model1)##
## Call:
## lm(formula = Price ~ RM_c + LSTAT_c + RM_c:LSTAT_c, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.2349 -2.6897 -0.6158 1.9663 31.6141
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.04226 0.23449 89.73 <2e-16 ***
## RM_c 3.56524 0.39263 9.08 <2e-16 ***
## LSTAT_c -0.85371 0.04006 -21.31 <2e-16 ***
## RM_c:LSTAT_c -0.48494 0.03459 -14.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.701 on 502 degrees of freedom
## Multiple R-squared: 0.7402, Adjusted R-squared: 0.7387
## F-statistic: 476.9 on 3 and 502 DF, p-value: < 2.2e-16vif(model1)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## RM_c LSTAT_c RM_c:LSTAT_c
## 1.738798 1.869362 1.167616Il modello lineare analizza l’effetto di RM_c, LSTAT_c e della loro interazione sul prezzo delle abitazioni. Tutti i coefficienti risultano altamente significativi (p < 2e-16), confermando la rilevanza sia degli effetti principali sia dell’interazione.
Il coefficiente positivo di RM_c indica che, a parità di LSTAT_c, un aumento del numero medio di stanze è associato a un incremento del prezzo, mentre il coefficiente negativo di LSTAT_c mostra che una maggiore percentuale di popolazione a basso status socioeconomico è legata a una riduzione del prezzo. L’interazione negativa suggerisce che l’effetto positivo di RM_c sul prezzo si attenua all’aumentare di LSTAT_c, evidenziando una relazione non puramente additiva tra le due variabili.
Il modello spiega circa il 74% della variabilità del prezzo (R^2 = 0.74), con un errore standard residuo di circa 4.7, indicando un buon adattamento ai dati. I valori di VIF < 2 confermano inoltre l’assenza di problemi di multicollinearità, anche in presenza del termine di interazione.
residui <- residuals(model1)
residuals_df <- data.frame(residui = residui)
# Calcolo media e deviazione standard dei residui
mu <- mean(residui)
sigma <- sd(residui)
ggplot(residuals_df, aes(x = residui)) +
geom_histogram(aes(y = ..density..),
color = 'black', fill = 'skyblue', binwidth = 1) +
stat_function(fun = dnorm, args = list(mean = mu, sd = sigma),
color = "red", size = 1) +
ggtitle("Istogramma dei residui") +
ylab("Densità")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
L’istogramma dei residui mostra una distribuzione circa simmetrica e prossima alla normalità, suggerendo che l’assunzione di normalità è ragionevolmente soddisfatta. Sono tuttavia presenti alcuni residui estremi, che potrebbero indicare outlier o osservazioni influenti e che meritano ulteriori verifiche. Per completare la diagnostica del modello, è opportuno analizzare il Q-Q plot e i residui rispetto ai valori predetti, così da valutare linearità, omoschedasticità e presenza di anomalie.
par(mfrow = c(1, 2))
# Residui vs valori predetti
plot(model1$fitted.values, residuals(model1),
xlab = "Valori Predetti",
ylab = "Residui",
main = "Residui vs Valori Predetti")
abline(h = 0, col = "red")
# QQ-plot dei residui
qqnorm(residuals(model1))
qqline(residuals(model1), col = "red")
par(mfrow = c(1, 1))Grafico residui vs valori predetti: i residui sono distribuiti casualmente attorno allo zero, senza pattern evidenti rispetto ai valori predetti, suggerendo che le assunzioni di linearità e omoschedasticità sono sostanzialmente rispettate.
Q-Q plot: mostra una buona aderenza alla normalità nella parte centrale della distribuzione, con lievi deviazioni nelle code dovute a residui estremi.
bptest(model1)##
## studentized Breusch-Pagan test
##
## data: model1
## BP = 6.4671, df = 3, p-value = 0.09097Il test di Breusch-Pagan ha dato BP = 6.4671 con p-value 0.091, superiore al livello di significatività 0.05. Pertanto, non ci sono evidenze sufficienti di eteroschedasticità e la varianza degli errori del modello model1 può essere considerata costante, rendendo affidabili le stime degli errori standard.
predicted_test <- predict(model1)
ggplot(data = df, aes(x = Price, y = predicted_test)) +
geom_point(color = "mediumpurple", alpha = 0.6, size = 2) +
geom_abline(slope = 1, intercept = 0, color = "red", size = 1) +
labs(
title = "Confronto tra valori osservati e predetti",
x = "Prezzi osservati",
y = "Prezzi predetti"
) +
theme_minimal()
Lo scatter plot mostra che i punti sono generalmente vicini alla linea y = x, indicando che il modello predice bene i prezzi nella maggior parte dei casi. Tuttavia come osservato anche durante l’EDA, per i prezzi più alti (40-50), poco rappresentati nel dataset, si osserva maggiore dispersione e predizioni meno precise. Nel complesso, la relazione è lineare e positiva, senza evidenti bias sistematici.
residuals <- df$Price - predicted_test # residui
RMSE <- sqrt(mean(residuals^2))
MAE <- mean(abs(residuals))
RMSE## [1] 4.682768MAE## [1] 3.217945I valori ottenuti per le metriche di errore RMSE (circa 4.68) e MAE (circa 3.22) indicano che il modello di regressione multipla con interazione ha una buona capacità predittiva rispetto alla variabile Price.n media, le previsioni si discostano di circa 3.22 unità (MAE), mentre l’RMSE mostra che gli scostamenti più grandi sono contenuti nell’ordine di 4.68. Considerando che Price varia tra 5 e 50, il modello fornisce stime abbastanza precise, anche se i prezzi più alti risultano più difficili da predire.