0. Caricamento dataset e libreria

library(lubridate)
library(moments)
library(ggplot2)
library(patchwork)
library(dplyr)
library(knitr)
real_estate_df <- read.csv("realestate_texas.csv", stringsAsFactors = TRUE)
#print(head(real_estate_df,5))
kable(head(real_estate_df,5))
city year month sales volume median_price listings months_inventory
Beaumont 2010 1 83 14.162 163800 1533 9.5
Beaumont 2010 2 108 17.690 138200 1586 10.0
Beaumont 2010 3 182 28.701 122400 1689 10.6
Beaumont 2010 4 200 26.819 123200 1708 10.6
Beaumont 2010 5 202 28.833 123100 1771 10.9 
print(dim(real_estate_df))
## [1] 240   8

Il dataset ha 8 grandezze:

  1. city (qualitativa nominale)
  2. year (qualitativa su scala ordinale)
  3. month (qualitativa su scala ordinale)
  4. sales (quantitativa discreta)- numero di case vendute in quel mese
  5. volume (quantitativa continua)- ricavi totali in milioni di dollari in quel mese
  6. median_price (quantitativa continua)- prezzo mediano in dollari delle case
  7. listings (quantitativa discreta)- numero totale di annunci attivi
  8. months_inventory (quantitativa continua)- mesi necessari per vendere tutte le inserzioni attuali

1. Analisi delle variabili

levels(real_estate_df$city)
## [1] "Beaumont"              "Bryan-College Station" "Tyler"                
## [4] "Wichita Falls"
unique(real_estate_df$year)
## [1] 2010 2011 2012 2013 2014
unique(real_estate_df$month)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12

Dovendo analizzare il risultato di quattro città per quattro anni, è meglio avere una colonna che possa combinare anno e mese in una formato tipo data.

real_estate_df$date <- paste(real_estate_df$year, real_estate_df$month, "01", sep = "-") %>% ymd() 

real_estate_df$date <- format(real_estate_df$date, "%b-%Y")

Molteplici sono le analisi che posso essere fatte con queste variabili, in particolare per quelle quantitative (sales, volume, median_price, listings e months_inventory) studieremo prima le loro distribuzioni attraverso istogrammi e boxplot, determineremo per ciascuna di queste variabili anche gli indici di posizione, variabilità e forma, successivamente studieremo come cambiano in funzione della città e se presentano dei trend nel tempo costruendo alcuni scatterplots.

Per le variabili city, month e year costruiremo delle distribuzioni di frequenze per vedere se tutte le modalità delle variabili sono equidistribuite nel dataset.

2. Indici di posizione, variabilità e forma

kable(table(real_estate_df$city))
Var1 Freq
Beaumont 60
Bryan-College Station 60
Tyler 60
Wichita Falls 60 
kable(table(real_estate_df$year))
Var1 Freq
2010 48
2011 48
2012 48
2013 48
2014 48 
kable(table(real_estate_df$year, real_estate_df$city))
Beaumont Bryan-College Station Tyler Wichita Falls
2010 12 12 12 12
2011 12 12 12 12
2012 12 12 12 12
2013 12 12 12 12
2014 12 12 12 12 
kable(table(real_estate_df$month))
Var1 Freq
1 20
2 20
3 20
4 20
5 20
6 20
7 20
8 20
9 20
10 20
11 20
12 20

Le distribuzioni di frequenze assolute ci mostrano che i dati sono equidistribuite per anno e per città, quindi avremo un totale 60 entries per ogni città, una per ogni mese dei 5 anni di analisi.

Calcoliamo gli indici di posizione, di variabilità e di forma per le variabili: sales, volume, median_price, listing e months_inventory

# Definizione funzione per il calcolo del coefficiente di variazione
CV <- function(x){
  return (sd(x)/mean(x))*100
}

# Lista delle colonne da analizzare
columns <- c("sales", "volume", "median_price", "listings", "months_inventory")

# Creazione di una lista per salvare i risultati
results <- list()

# Ciclo sulle colonne
for (col in columns) {
  
  data <- real_estate_df[[col]] 
    
  # Calcola le statistiche base
  stats <- round(summary(data),2)  
    
  # Calcola ulteriori statistiche
  iqr <- round(IQR(data),2)  
  varianza <- round(var(data),2) 
  sd <- round(sd(data),2)  
  cv <- round(CV(data),2)
  skew <- round(skewness(data),2)  
  kurt <- round(kurtosis(data),2)-3  
    
    # Combina tutte le statistiche in una riga
  stats_row <- c(stats, IQR = iqr, Varianza = varianza, SD = sd, CV = cv, Skewness = skew, Kurtosis = kurt)
    
   
  results[[col]] <- as.data.frame(t(stats_row))

}

# Combina i risultati in un unico data.frame
results_table <- do.call(rbind, results)
results_table <- cbind(Variable = rownames(results_table), results_table)  
rownames(results_table) <- NULL  

# Aggiorna i nomi delle colonne per includere le nuove statistiche
colnames(results_table) <- c("Variable", "Min", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max", 
                             "IQR", "Variance", "SD", "CV", "Skewness", "Kurtosis")

kable(results_table)
Variable Min 1st Qu. Median Mean 3rd Qu. Max IQR Variance SD CV Skewness Kurtosis
sales 79.00 127.00 175.50 192.29 247.00 423.00 120.00 6.34430e+03 79.65 0.41 0.72 -0.31
volume 8.17 17.66 27.06 31.01 40.89 83.55 23.23 2.77270e+02 16.65 0.54 0.88 0.18
median_price 73800.00 117300.00 134500.00 132665.42 150050.00 180000.00 32750.00 5.13573e+08 22662.15 0.17 -0.36 -0.62
listings 743.00 1026.50 1618.50 1738.02 2056.00 3296.00 1029.50 5.66569e+05 752.71 0.43 0.65 -0.79
months_inventory 3.40 7.80 8.95 9.19 10.95 14.90 3.15 5.31000e+00 2.30 0.25 0.04 -0.17

Le 5 grandezze che abbiamo considerato hanno domini molto diversi questo a causa delle diverse unità di misura ma anche per l’informazione che rappresentano. Ci sono delle tendenze, infatti possiamo notare che tutte, tranne median_price, presentano un’assimetria positiva.

3. Identificazione delle variabili con maggiore variabilità e asimmetria

Le grandezze sales, volume, listings e months_inventory mostrano tutte unìasimmetria positiva in quanto la media è maggiore della mediana, mentre median_price mostra una asimmetria negativa.

Leggendo la colonna del coefficiente di variabilità possiamo notare che volume è la variabile con la più alta variabilità ed è anche quella che presente un coefficiente di asimmetria più alto, è anche l’unica con una Curtosi positiva (0.18) valore che ci dice come la variabile presenti la tendenza ad avere una distribuzione più stretta rispetto a quella della distribuzione normale. Validiamo queste osservazioni anche attraverso la rappresentazione grafica con istogrammi e boxplot.

Costruiamo gli istogrammi di queste sei grandezze:

# Definizione per il calcolo del numero di bin secondo il criterio di Rice
rice_bins <- function(data) {
  n <- length(data) # Numero di osservazioni
  return(2 * n^(1/3))
}
#Usiamo attach per accedere più facilmente alla variabili, dovendo costruire molteplici grafici

attach(real_estate_df)

# Lista per salvare i grafici
plots <- list()

# Ciclo per generare gli istogrammi
for (col in columns) {
  
  data <- real_estate_df[[col]]
  
  # Calcola il numero di bin e il binwidth
  bins <- rice_bins(data)
  binwidth <- diff(range(data, na.rm = TRUE)) / bins
  
  # Crea l'istogramma
  plot <- ggplot(real_estate_df) +
    geom_histogram(aes_string(x = col),
                   binwidth = binwidth,
                   col = "black",
                   fill = "lightblue") +
    labs(title = paste("HS per", col),
         x = col,
         y = "Count") +
    theme_minimal()
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

# Disposizione dei grafici in una griglia 2x3
final_plot <- (plots[[1]] | plots[[2]] | plots[[3]]) /
              (plots[[4]] | plots[[5]])

# Mostra i grafici
print(final_plot)

Costruiamo i boxplot di queste cinque grandezze:

# Lista per salvare i grafici
plots <- list()

# Ciclo per generare gli istogrammi
for (col in columns) {
  
  data <- real_estate_df[[col]]
  
  # Crea l'istogramma
  plot <- ggplot(real_estate_df) +
    geom_boxplot(aes_string(y = col),
                   col = "black",
                   fill = "lightblue") +
    labs(title = paste("BP per", col),
         y = col) +
    theme_minimal()+
    theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
          axis.ticks.x = element_blank())
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

# Disposizione dei grafici in una griglia 2x3
final_plot <- (plots[[1]] | plots[[2]] | plots[[3]]) /
              (plots[[4]] | plots[[5]])

# Mostra i grafici
print(final_plot)

Osservando i boxplot vediamo che volume è l’unica grandezza che presenta degli outliers, la presenza dei quali potrebbe giustificare l’asimetria registrata dal relativo coefficiente e la sua maggiore variabilità.

#4. Creazione di classi per una variabile quantitativa

Selezioniamo la variabile sales per costruire una distribuzione delle frequenze e determinare l’indice di eterogeneità di Gini.

#Definiamo la funzione per calcolare l'indice di eterogeneità di Gini
gini.index <- function(x){
  ni = table(x)
  fi = ni/length(x)
  fi2 = fi^2
  J = length(table(x))
  gini = 1 - sum(fi2)
  gini.normalizzato = gini/((J-1)/J)
  return(gini.normalizzato)
}

step_sales=(max(sales)-min(sales))/rice_bins(sales)
sales_cl <- cut(sales, seq(min(sales)-1, max(sales)+step_sales, step_sales))
n = length(sales_cl)

ni = table(sales_cl)
fi = table(sales_cl)/n
Ni = cumsum(ni)
Fi = Ni/n

dist_freq_sales <- as.data.frame(
  cbind(ni, fi, Ni, Fi)
)

dist_freq_sales$Intervalli <- rownames(dist_freq_sales)
dist_freq_sales$Intervalli <- factor(
  dist_freq_sales$Intervalli,  # Intervalli già presenti come stringhe
  levels = dist_freq_sales$Intervalli  # Ordine delle righe di dist_freq_sales
)

ggplot(dist_freq_sales, aes(x = Intervalli, y = fi)) +
  geom_col(fill = "steelblue") +
  xlab("Intervalli di Vendite") +
  ylab("Frequenza Relativa") +
  ggtitle("Distribuzione delle Vendite") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

round(gini.index(dist_freq_sales$fi),2)
## [1] 0.99

Studiando il Colplot della distribuzione di frequenze relative della variabile sales, osserviamo che 5 delle 13 classi in cui abbiamo suddiviso la variabile, usando il criterio di Rice, presentano una frequenza tra il 10% e il 20%. Il campione quindi si presenta molto etereogeneo, come confermato anche dall’indice di eterogeneità di Gini che presenta un valore pari allo 0.9941.

5. Calcolo della probabilità

Vogliamo determinare alcune probabilità: la probabilità di avere un dato dalla città di “Beaumont”, la probabilità di avere un dato estratto dal mese di Luglio e la probabilità di estrarre un dato di Dicembre 2012.

# Le tre probabilità che dobbiamo determinare possiamo tutte farle come casi favorevoli diviso casi possibile, cioè applicando la definisione di probabilità.

# Probabilità Beaumont
city_counts <- table(city)
p_Beaumont <- city_counts["Beaumont"]/length(city)
print(p_Beaumont)
## Beaumont 
##     0.25
# Probabilità Luglio
month_counts <- table(month)
p_July <- round(month_counts["7"]/length(month),2)
print(p_July)
##    7 
## 0.08
# Probabilità Dicembre 2012
month_year_counts <-table(date)
p_dic_2012 <- round(month_year_counts["Dec-2012"]/length(date),2)
print(p_dic_2012) 
## Dec-2012 
##     0.02

6. Creazione di nuove variabili

Vogliamo aggiungere una colonna al nostro dataset con il prezzo medio degli immobili. Per questo scopo useremo le features volumes e sales il loro rapporto moltiplicato per 10^6 ci permetterà di determinare il valore medio della vendita di una casa in dollari.

Dopo aver definito la nuova variabile, mostreremo il suo istogramma.

real_estate_df$mean_price <- volume*10^6/sales

bins <- rice_bins(real_estate_df$mean_price)
binwidth <- diff(range(real_estate_df$mean_price, na.rm = TRUE)) / bins

ggplot() +
  geom_histogram(aes(x = real_estate_df$mean_price),
                  binwidth = binwidth,
                  col = "black",
                  fill = "lightblue") +
  labs(title = paste("HS per Mean Price"),
        x = "Mean Price",
        y = "Count") +
  theme_minimal()

summary(real_estate_df$mean_price)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   97010  132939  156588  154320  173915  213234

Dall’istogramma possiamo osservare che c’è una zona centrale (intorno a 160000 $) dove si concentrano la maggior parte delle case ma che esiste anche un picco secondario attorno a 120000 $, nonostante questo la distribuzione sembra avere una discreta simmetria, il valore mediano e la media sono molto vicini, il coefficiente di asimmetria in valore assoluto è minore di 0.1.

round(skewness(real_estate_df$mean_price),2)
## [1] -0.07

Vogliamo creare una variabile che descriva l’efficacia degli annunci di vendita. Potremmo dire che un annuncio è efficace se la casa è venduta in un minore quantitativo di tempo. sostanzialmente una campagna publicitaria è efficacia quando a parità di numero di annunci c’è una frequenza più alta di case vendute al mese. Per i dati che abbiamo possiamo definire l’indice che misura l’efficacia degli annunci di vendita come il rapporto tra il numero di annunci attivati quel mese - listing - e quanti mesi sono stati necessari per venderli tutti - months_inventory, più alto sarà questo valore più la campagna publicitaria è stata efficace. Prima di definire la nuova variabile osserviamo se effettivamente listings e month_inventory mostrano una qualche dipendenza:

ggplot(data = real_estate_df, aes(x = listings, y = months_inventory)) +
  geom_point(color = "steelblue", size = 3, alpha = 0.7) +
  xlab("Listings") +
  ylab("Months Inventory") +
  ggtitle("Scatterplot di Listings vs Months Inventory") +
  theme_minimal()

Sono presenti tre serie di dati, segno di una dipendenza da qualche altro fattore probabilmente la città. Quindi riproduciamo il grafico precedente esplicitando la città di appartenza di ciascun dato (questo ci porta in modo naturale ad iniziare la nostra analisi condizionata).

7. Analisi condizionata

ggplot(data = real_estate_df, aes(x = listings, y = months_inventory, color = city)) +
  geom_point(size = 3, alpha = 0.7) +
  xlab("Listings") +
  ylab("Months Inventory") +
  ggtitle("Scatterplot di Listings vs Months Inventory (per città)") +
  theme_minimal()+
  theme(legend.position = "right")

Come sospettavano, esplicitando la città di appartenenza, le due grandezze mostrano una molto probabile dipendenza lineare, quindi procediamo a definire il loro rapporto come l’indice di efficienza della campagna pubblicitaria:

#Definizione della nuova variabile
real_estate_df$eff_ads <- listings/months_inventory

È interessante provare a vedere se c’è un correlazione tra il numero di annunci attivi e l’efficacia degli annunci stessi:

ggplot(data = real_estate_df, aes(x = listings, y = eff_ads, color = city)) +
  geom_point(size = 3, alpha = 0.7) +
  xlab("Listings") +
  ylab("Effective Ads") +
  ggtitle("Scatterplot di Listings vs Effective Ads") +
  theme_minimal()+
  theme(legend.position = "right")

Il risultato è interessante, effettivamente c’è una certa correlazione tra il numero di annunci attivi e l’efficacia della campagna stessa. Proviamo a vedere invece la correlazione con i mesi necessari per vendere tutti gli annunci.

ggplot(data = real_estate_df, aes(x = months_inventory, y = eff_ads, color = city)) +
  geom_point( size = 3, alpha = 0.7) +
  xlab("Months Inventory") +
  ylab("Effective Ads") +
  ggtitle("Scatterplot di Months Inventory vs Effective Ads") +
  theme_minimal()+
  theme(legend.position = "right")

All’interno della serie di dati relativi ad una stessa città si vede una correlazione inversa tra l’efficacia della campagna è i mesi necessari per vendere tutti gli annunci.

I due scatterplot ci testimoniano che la variabile che abbiamo definito riesce a descrivere l’efficacia di una campagna pubblicitaria.

Creiamo un summary di media e deviazione standard delle nostre nuove variabili in funzione della città.

kable(real_estate_df %>%
  group_by(city) %>%
  summarise(media_price = mean(mean_price),
            dev.st_price = sd(mean_price),
            media_eff_ads = mean(eff_ads),
            dev.st_eff_ads = sd(eff_ads)))
city media_price dev.st_price media_eff_ads dev.st_eff_ads
Beaumont 146640.4 11232.13 172.0570 23.121131
Bryan-College Station 183534.3 15149.35 199.5783 33.744500
Tyler 167676.8 12350.51 260.9829 29.293672
Wichita Falls 119430.0 11398.48 116.7364 6.326538

Il summary prodotto mostra la differenza tra i prezzi delle case nelle quattro città e come anche l’efficacia della campagna pubblicitaria ha una certa dipendenza dalla città stessa.

Per chiarezza possiamo mostrare questo medesimo risultato attraverso dei boxplot:

ggplot() +
    geom_boxplot(aes_string(y = real_estate_df$mean_price, color = city),
                   fill = "lightblue") +
    labs(title = paste("BP per Mean Price e città"),
         y = "Mean Price",
         color = "Città") +
    theme_minimal()+
    theme(legend.position = "bottom")+
    theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
          axis.ticks.x = element_blank())

ggplot() +
    geom_boxplot(aes_string(y = real_estate_df$eff_ads, color = city),
                   fill = "lightblue") +
    labs(title = paste("BP per Eff Ads e città"),
         y = "Eff Ads",
         color="Città") +
    theme_minimal()+
    theme(legend.position = "bottom")+
    theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
          axis.ticks.x = element_blank())

Riproduciamo i boxplot del punto 3, aggiungendo la dipendenza dalla città e dall’anno:

columns <- c("sales", "volume", "median_price", "listings", "months_inventory")

for (col in columns) {
  
  data <- real_estate_df[[col]]
  
  plot <- ggplot(real_estate_df) +
  geom_boxplot(aes_string(y = col, x = "city", fill = "factor(year)")) +
  labs(title = paste("Boxplot per città, anno e", col),
       y = col,
       x = "Città",
       fill = "Anno") +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_fill_brewer(palette = "Set3")
  
  print(plot)
}

Dal boxplot del volume possiamo vedere come la società sia cresciuta negli anni, infatti aumenta la mediana di volume ma ha aumentato anche la variabilità dei suoi ricavi nei diversi anni, infatti il “box” si è allargato costantemente.

Per osservare se c’è una dipendenza dall’anno o dal mese oltre che dalla città, rappresentiamo alcune variabili in funzione del tempo:

# Converte la colonna 'date' in formato Date
real_estate_df$date2 <- as.Date(paste("01", real_estate_df$date, sep = "-"), format = "%d-%b-%Y")

columns <- c("sales", "median_price", "listings", "eff_ads")

for (col in columns) {
  # Crea il grafico
  p <- ggplot(real_estate_df, aes(x = date2, y = !!sym(col), color = city, group = city)) +
    geom_line(linewidth = 1) +
    geom_point(size = 2) +
    xlab("Mese") +
    ylab(col) +
    ggtitle(sprintf("Andamento di %s nel tempo e per Città", col)) +
    scale_x_date(
      breaks = seq(min(real_estate_df$date2), 
                   max(real_estate_df$date2), 
                   by = "4 months"),  # Intervallo ogni 4 mesi
      labels = scales::date_format("%b %Y")  # Formato "Gennaio 2022"
    ) +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "bottom"
    )

  # Stampa il grafico
  print(p)
}

Dal grafico della variabile sales si può osservare come il numero di vendite abbia una certà periodicità infatti raggiungono il picco nel periodo estivo e il minimo durante l’inverno, per osservare meglio questo fenomeno rappresentiamolo attraverso un grafico a barre. Questa periodicità si riflette anche il listing.

La variabile median_price mostra per Beaumomnt e Wichita Falls una certa caoticità mentre per Tyler e Bryan-College station è leggermente crescente.

Mentre la variabile eff_ads per tutta le città tranne Wichita Falls mostra una crescita costante segno che le campagne pubblicitarie hanno avuto un maggior impatto.

# Crea un vettore con i nomi dei mesi
month_names <- c("Jann", "Feb", "Mar", "Apr", "May", "Jun", 
                 "Jul", "Aug", "Sep", "Oct", "Nov", "Dic")

# Converte il numero del mese in nome del mese
real_estate_df$month_name <- factor(real_estate_df$month, levels = 1:12, labels = month_names)

ggplot(data = real_estate_df) +
  geom_col(aes(x = month_name, y = volume, fill = city),
           position = "stack") +
  labs(title = "Distribuzione delle vendite nei mesi",
       x = "Mese",
       y = "Vendite totali") +
  scale_y_continuous(labels = scales::comma) + 
  theme_classic() +
  theme(legend.position = "bottom")

ggplot(data = real_estate_df) +
  geom_col(aes(x = month_name, y = volume, fill = city),
           position = "fill") +
  labs(title = "Distribuzione relativa delle vendite nei mesi",
       x = "Mese",
       y = "Frequenze relative") +
  scale_y_continuous(labels = scales::percent) +  # Mostra le frequenze come percentuali
  theme_classic() +
  theme(legend.position = "bottom")

Questi ultimi due grafici mettono in evidenza in modo chiaro come il mercato immobiliare sia costate durante l’anno nelle città di Wichita Falls e Beaumont mentre per Bryan-College Station, soprattutto, ma in parte anche per la città di Tyler c’è una dipendenza dalla stagionalità.

8. Conclusioni

L’analisi svolta è stata molto ricca e molteplici sono le considerazioni che possono essere fatte:

  1. Il mercato di Wichita Falls è quello meno sviluppato sotto tutti gli indici, potrebbe essere una città in cui fare investimento cercando di aumentare il numero di case in vendita.

  2. La città di Bryan-College Station sembra strategica per il business essendo quella che mostra un valore mediano maggiore.

  3. Nelle città di Wichita Falls e Beaumont ha senso intensificare la pubblicità nel corso di tutto l’anno mentre per le altre due devono essere sfruttati soprattutto i mesi estivi piuttosto che quelli invernali che presentano una maggior difficoltà nel concludere accordi di vendita.