Nota: è possibile utilizzare la table of content scorrevole sulla sinistra per andare alle specifiche sezioni.

Analisi descrittiva del mercato immobiliare del Texas

In questo piccolo progetto andrò ad effettuare una semplice analisi descrittiva di un dataset .csv di prova che contiene report dei dati relativi alle vendite immobiliari in Texas.

Approfitterò anche per provare un po’ di differenti funzionalità e pacchetti aggiuntivi.

In uno studio formale sarebbe ovviamente più appropriato mantenere una certa omogeneità, sia per chiarezza sia per un più efficiente utilizzo di risorse.

~In questo scenario cercherò di sperimentare (come chiesto anche nei dettagli dell’esercitazione) nei limiti delle mie 
(molto modeste) capacità di programmazione.~

Nota: tutti i pacchetti caricati con library() vanno ovviamente prima installati, nel caso l’utente non lo abbia già fatto, tramite il comando install.packages(“nomepacchetto”).

1: Individuazione e descrizione delle tipologie di variabili

#Caricamento del dataset

data <- read.csv("realestate_texas.csv", sep = ",")

#Visualizzazione grezza della struttura del dataset:

str(data)
## 'data.frame':    240 obs. of  8 variables:
##  $ city            : chr  "Beaumont" "Beaumont" "Beaumont" "Beaumont" ...
##  $ year            : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
##  $ month           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sales           : int  83 108 182 200 202 189 164 174 124 150 ...
##  $ volume          : num  14.2 17.7 28.7 26.8 28.8 ...
##  $ median_price    : num  163800 138200 122400 123200 123100 ...
##  $ listings        : int  1533 1586 1689 1708 1771 1803 1857 1830 1829 1779 ...
##  $ months_inventory: num  9.5 10 10.6 10.6 10.9 11.1 11.7 11.6 11.7 11.5 ...

Il dataset contiene 240 “elementi” (osservazioni) e 8 diverse variabili. Esaminiamone la tipologia e le modalità (nel caso delle variabili non categoriche ne prenderemo solo alcune casualmente).

  1. City: Formato stringa (chr), variabile qualitativa su scala nominale. Ricaviamo le possibili modalità uniche che essa può assumere nel dataset (tramite una semplice funzione della libreria dplyr):

    #utilizzato l'argomento message=false per evitare di riportare l'output del caricamento della libreria, sarà lo stesso per gli altri blocchi di caricamento librerie
library(dplyr)
    cities <- data %>% group_by(city)

cities %>% group_keys()
    ## # A tibble: 4 × 1
##   city                 
##   <chr>                
## 1 Beaumont             
## 2 Bryan-College Station
## 3 Tyler                
## 4 Wichita Falls

    Queste sono semplicemente le località nelle quali si trovano gli immobili venduti ovvero Beaumont, Bryan-College Station, Tyler e Wichita Falls.

  2. Year: Formato numero intero (int). Quando si tratta di tempo, l’effettiva tipologia della variabile può essere ambigua. Visualizziamo le modalità:

    years <- data %>% group_by(year)

#utilizzo di print.data.frame per rendere l'output più pulito
print.data.frame(years %>% group_keys())  
    ##   year
## 1 2010
## 2 2011
## 3 2012
## 4 2013
## 5 2014

    Formalmente, saremmo in presenza di una variabile quantitativa discreta (solo nel senso che indichiamo gli anni interi, ma di per sé sarebbe continua) su scala ad intervalli. In questo scenario è probabilmente più utile e appropriato trattarla come una variabile qualitativa (o categorica). Sono registrate le vendite su un periodo di 5 anni consecutivi che va da inizio 2010 a fine 2014, è opportuno associare gli anni a “categorie”.

  3. Month: Formato numero intero (int). Abbiamo un’unità normalmente espressa in parole (i nomi dei mesi) ma qui espressa sottoforma di numeri. Si può intuire che abbia la stessa funzione della precedente variabile. Osserviamo, per completezza, le modalità della variabile:

    ##    month
## 1      1
## 2      2
## 3      3
## 4      4
## 5      5
## 6      6
## 7      7
## 8      8
## 9      9
## 10    10
## 11    11
## 12    12

    I mesi dell’anno sono “codificati” con numeri da 1 a 12. La variabile può essere trattata come una variabile qualitativa ordinale, nello specifico ordinale ciclica in quanto tecnicamente non esiste una reale modalità iniziale (o finale), se non nel contesto di un’analisi annuale ordinata. Chiaramente sarebbe possibile trattarla come quantitativa nel caso si dovessero effettuare operazioni numeriche (come calcolare una distanza temporale e trarre eventuali altre informazioni correlate e coerenti nel dataset).

  4. Sales: Formato numero intero (int). Vediamo qualche modalità:

    sample(data$sales, 10)
    ##  [1] 127 202 170 182 388  90 163 111 403 196

    Una semplice variabile quantitativa discreta su scala di rapporti. Indicherà il numero di vendite effettuato nello specifico mese-anno di riferimento.

    Questa variabile, come le altre successive, è legata alle 3 precedenti (nel senso che con opportuni strumenti possibile indicare le prime 3 variabili come “input di ricerca” e ottenere le altre 5 variabili come informazioni).
    Si può anche selezionare solo una o due di loro 3 ed ottenere un dataset ridotto invece che una singola riga.
    Un esempio grezzo per ottenere le informazioni per Beaumont, anno 2011 e mese di maggio:

    c <- "Beaumont"
y <- 2011
m <- 5

info <- data %>%
  filter(city == c, year == y, month == m) %>%
  select(sales, volume, median_price, listings, months_inventory)
#abbiamo "filtrato" per le 3 variabili indicate prima, e selezionato le altre per la visualizzazione

print(info)
    ##   sales volume median_price listings months_inventory
## 1   143 19.684       120700     1832             12.4

  5. Volume: Formato numero decimale (num).

    ##  [1] 25.487 63.759 28.833 29.117 35.395 19.684 16.908 40.306 23.029 49.914

    Variabile quantitativa continua su scala di rapporti. Indica il valore delle vendite in milioni di dollari nel mese-anno di riferimento.

  6. Median_price: Formato numero decimale (num).

    ##  [1] 143100 118200 111100 126200  99300 155200  91700 134200  95000 163700

    Variabile quantitativa continua* su scala di rapporti. Indica il prezzo mediano di vendita nel mese-anno di riferimento.

    *Da notare comunque che si tratta numeri interi arrotondati alle centinaia, seppur inseriti nel dataset in formato decimale.

    In un certo senso apparirebbe come una variabile discreta, ma sarebbe formalmente incorretto.

  7. Listings: Formato numero intero (int).

    ##  [1] 1041 3272 1562 3267 1016 2897 1486 1586  907 2938

    Variabile quantitativa discreta su scala di rapporti. Indica il numero di annunci attivi nel periodo di riferimento.

  8. Months_inventory: Formato numero decimale (num).

    ##  [1] 13.8  9.3  7.9  9.0  8.1  8.5  8.6  4.8  9.0 11.6

    Variabile quantitativa continua su scala di rapporti. Indica il tempo in mesi che sarebbe necessario, col rateo di vendite stimato nel periodo di riferimento, per vendere tutte le inserzioni ancora presenti.

2: Calcolo degli indici e distribuzioni di frequenze

Per le prime 3 variabili, essendo qualitative, calcoliamo le distribuzioni di frequenze. Per le altre calcoleremo i vari indici.

CITY

#Memorizziamo il numero di righe (oggetti) del dataset
n <- dim(data)[1]

#Calcoliamo frequenze assolute, relative, e le loro cumulate
#Tabuliamo i dati della variabile city con la funzione table
city_ni <- table(data$city)    #freq assolute
city_fi <- table(data$city)/n  #freq relative
city_Ni <- cumsum(city_ni)     #somma cumulativa delle freq assolute
city_Fi <- city_Ni/n           #somma cumulativa delle freq relative

#concateniamo i valori ottenuti in un vettore
cbind(city_ni, city_fi, city_Ni, city_Fi)  
##                       city_ni city_fi city_Ni city_Fi
## Beaumont                   60    0.25      60    0.25
## Bryan-College Station      60    0.25     120    0.50
## Tyler                      60    0.25     180    0.75
## Wichita Falls              60    0.25     240    1.00

Notiamo che nel dataset le vendite sono equamente distibuite (city_ni, che rappresenta la frequenza assoluta, è 60 per tutte le modalità). Abbiamo una variabile quadrimodale.

(Ciò ha senso considerando che, come vedremo dopo, per tutte le città sono riportati i dati relativi a un periodo di 60 mesi consecutivi.)

Osserviamo meglio per mezzo di un basilare grafico:

city_distrib <- as.data.frame(cbind(city_ni, city_fi, city_Ni, city_Fi))
#Trasformare i vettori in un dataframe è consigliabile e utile in molti casi,
#ad esempio nel caso si vogliano esportare i risultati in file csv o simili,
#anche se in questo specifico scenario di presentazione non è necessario.

#Per tenere conto di questa operazione, è necessario specificare come primo
#argomento il dataframe in questione (e la colonna da usare, con la notazione $). 
#In questo caso funzionerebbe anche indicando direttamente il vettore city_ni.

par(las=1)    #per rendere i valori dell'asse y orizzontali
barplot(city_distrib$city_ni,
        main = "Distribuzione di frequenza (città)",
        xlab = "Città",
        ylab = "Frequenza assoluta",
        col = "gray",
        names.arg = rownames(city_distrib))

YEAR

Questa volta utilizzeremo la libreria qdap per eseguire le stesse operazioni ma più rapidamente.

library(qdap)

Con il comando dist_tab() di qdap è possibile ottenere direttamente una distribuzione di frequenza (specificando opportunamente la variabile del dataset). Le frequenze relative saranno indicate in formato percentuale.

year_distrib <- dist_tab(data$year)

year_distrib
##   interval freq cum.freq percent cum.percent
## 1     2010   48       48      20          20
## 2     2011   48       96      20          40
## 3     2012   48      144      20          60
## 4     2013   48      192      20          80
## 5     2014   48      240      20         100

Noteremo che anche qui abbiamo una distribuzione equimodale, pentamodale in questo caso. La frequenza sarà di 48 per tutti e 5 gli anni.

Ciò ha senso se consideriamo che, come vedremo, sono riportati i dati per ogni mese di 12 mesi (per ogni anno) per 5 anni.

Se prendiamo un anno, per ogni città ci saranno 12 report (quindi 12 “oggetti” nel dataset).

Se prendiamo tutte e 4 le città, cumulativamente ci saranno 48 report annuali.

Utilizziamo la stessa tipologia di grafico basilare per visualizzare:

par(las=1)  
barplot(year_distrib$freq,
        main = "Distribuzione di frequenza (anni)",
        xlab = "Anno",
        ylab = "Frequenza assoluta",
        ylim = c(0, 50),  #impostiamo manualmente la cima dell'asse y
        col = "gray",
        names.arg = year_distrib$interval)

MONTH

month_distrib <- dist_tab(data$month)

month_distrib
##    interval freq cum.freq percent cum.percent
## 1         1   20       20    8.33        8.33
## 2         2   20       40    8.33       16.67
## 3         3   20       60    8.33       25.00
## 4         4   20       80    8.33       33.33
## 5         5   20      100    8.33       41.67
## 6         6   20      120    8.33       50.00
## 7         7   20      140    8.33       58.33
## 8         8   20      160    8.33       66.67
## 9         9   20      180    8.33       75.00
## 10       10   20      200    8.33       83.33
## 11       11   20      220    8.33       91.67
## 12       12   20      240    8.33      100.00

Anche in questo caso, distribuzione equimodale. La frequenza è di 20 per ogni mese.

Essendo il periodo di report di 5 anni, per ogni città ci saranno 5 report di mese corrispondente (ovviamente di anno diverso).

Abbiamo 4 città, dunque 5*4=20.

par(las=1)
barplot(month_distrib$freq,
        main = "Distribuzione di frequenza (mesi)",
        xlab = "Mese",
        ylab = "Frequenza assoluta",
        ylim = c(0, 20),
        col = "gray",
        names.arg = c("gen", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dic")) #sovrascriviamo la notazione in numeri con abbreviazioni dei mesi

—

Da questo momento, passeremo al calcolo degli indici (posizione, variabilità, forma) delle altre 5 variabili.

SALES

Partiamo dagli indici di posizione.

#QUARTILI (comprende minimo, massimo e mediana)
quantile(data$sales)
##    0%   25%   50%   75%  100% 
##  79.0 127.0 175.5 247.0 423.0

0% e 100% corrisponderanno al minimo e al massimo, 50% corrisponderà alla mediana.

Per formalità, verifichiamoli manualmente:

#Utilizzo del comando paste per concatenare in modo semplice
#le stringhe e i risultati dei calcoli nelle funzioni richiamate
#Fun fact: venendo dal corso Python ero abituato ad eseguire
#operazioni di questo tipo con print, ma R apparentemente 
#non permette la stessa sintassi e c'è voluta un'ora solo per 
#capire ci volesse un altro comando
#Non so se ci siano modi migliori, ma questo è adeguato

paste("Il minimo è:", min(data$sales), " la mediana è:", median(data$sales), " il massimo è:", max(data$sales))
## [1] "Il minimo è: 79  la mediana è: 175.5  il massimo è: 423"

Calcoliamo la media (aritmetica):

paste("Media (aritmetica):", mean(data$sales))
## [1] "Media (aritmetica): 192.291666666667"

Non dimentichiamo la moda. Essendoci un numero non basso di osservazioni, utilizzare table non è molto comodo.

Si potrebbe utilizzare una piccola funzione creata manualmente, ad esempio:

mode <- function(x) {
  u <- unique(x)    #tutte le modalità uniche possibili dell'argomento
  tab <- tabulate(match(x, u))   #confronto tra tutte le modalità della variabile e quelle uniche
  u[tab == max(tab)]
}

#Ci ho messo più del dovuto a trovare questa combinazione

Questa funzione crea una tabella di frequenza e poi trova il valore (o i valori) che hanno la frequenza massima.

mode(data$sales)
## [1] 124

Ma visto che siamo sfaticati potremmo voler utilizzare un pacchetto che offre una funzione dedicata, come DescTools.

library(DescTools)
Mode(data$sales)
## [1] 124
## attr(,"freq")
## [1] 5
#Il primo valore riportato in output è la moda, l'ultimo indica quante volte si ripete

Notiamo che la moda è inferiore della mediana, la quale è inferiore della media. Ci ritorneremo a breve.

Passiamo agli indici di variabilità.

IQR (range interquartile, ovvero il distacco tra terzo e primo quartile):

paste("IQR:", IQR(data$sales))
## [1] "IQR: 120"

Varianza (o momento secondo, ovvero il rapporto tra la sommatoria dagli scarti dalla media al quadrato e il numero totale di osservazioni):

var(data$sales)
## [1] 6344.3
#tengo conto del fatto che le funzioni dirette restituiscono un valore leggermente diverso
#da quelli calcolati manualmente, ma ai fini di questo studio non è rilevante

Deviazione standard (radice quadrata della varianza):

sales_sd <- sd(data$sales)

print(sales_sd)
## [1] 79.65111

Coefficiente di variazione (rapporto tra deviazione standard e media):

sales_CV <- sd(data$sales)/mean(data$sales)

#nota, facoltativamente si moltiplica per 100 nel caso si desideri
#utilizzare un valore percentuale

print(sales_CV)
## [1] 0.4142203

Ora passiamo agli indici di forma:

Utilizziamo il pacchetto moments per questa sezione.

library(moments)

(A)Simmetria (skewness) ovvero il rapporto del rapporto tra il cubo del momento terzo della distribuzione (rapporto tra sommatoria degli scarti dalla media al cubo e il numero di osservazioni) e il cubo della deviazione standard:

sales_skew <- skewness(data$sales)

print(sales_skew)
## [1] 0.718104

Essendo maggiore di 0, avremo una distribuzione asimmetrica positiva.

Curtosi (ovvero il rapporto tra momento quarto e deviazione standard alla quarta):

#ricordiamoci di sottrarre 3 al valore ottenuto dalla funzione dedicata,
#essendo essa centrata in 3 e non in 0

sales_kurt <- kurtosis(data$sales)-3

print(sales_kurt)
## [1] -0.3131764

Essendo minore di 0 avremo una distribuzione platicurtica, ovvero più “appiattita” (almeno mediamente) rispetto alla distribuzione normale teorica.

Per farci un’idea, vediamo una rappresentazione grezza della distribuzione (ci interessa principalmente la “forma” del grafico):

plot(density(data$sales), lwd=2)
abline(v=range(data$sales), col=2)

x <- seq(min(data$sales), max(data$sales), length = 200)
y <- dnorm(x, mean = mean(data$sales), sd = sd(data$sales))
lines(x, y, col = "cyan3", lwd = 1)

~(Le due linee rosse sono state inserite per delimitare minimo e massimo, le parti di curva che le superano possono essere ignorate. 
La linea azzurrina rappresenta la distribuzione normale teorica, data stessa media e deviazione standard.)
~
Si può dire che il numero di vendite tipico di una città in un mese è molto più vicino al minimo che al massimo.

Vediamo anche l’andamento cumulativo delle vendite (includendo tutte le città) nell’intero periodo d’analisi per mezzo della funzione aggregate:

#Calcoliamo le vendite cumulative per città e mese (di ogni anno)
cum_sales <- aggregate(data$sales ~ data$month + data$year, data = data, FUN = sum)

#l'output sarà ordinato per mesi ciclicamente, andando avanti di anno in anno
print(cum_sales)
##    data$month data$year data$sales
## 1           1      2010        421
## 2           2      2010        487
## 3           3      2010        755
## 4           4      2010        916
## 5           5      2010        931
## 6           6      2010        866
## 7           7      2010        712
## 8           8      2010        738
## 9           9      2010        598
## 10         10      2010        565
## 11         11      2010        503
## 12         12      2010        604
## 13          1      2011        425
## 14          2      2011        469
## 15          3      2011        668
## 16          4      2011        716
## 17          5      2011        780
## 18          6      2011        885
## 19          7      2011        812
## 20          8      2011        786
## 21          9      2011        629
## 22         10      2011        594
## 23         11      2011        549
## 24         12      2011        565
## 25          1      2012        499
## 26          2      2012        574
## 27          3      2012        711
## 28          4      2012        747
## 29          5      2012        882
## 30          6      2012        898
## 31          7      2012        928
## 32          8      2012        954
## 33          9      2012        707
## 34         10      2012        742
## 35         11      2012        650
## 36         12      2012        643
## 37          1      2013        576
## 38          2      2013        593
## 39          3      2013        814
## 40          4      2013        878
## 41          5      2013       1057
## 42          6      2013       1045
## 43          7      2013       1127
## 44          8      2013       1107
## 45          9      2013        814
## 46         10      2013        738
## 47         11      2013        690
## 48         12      2013        733
## 49          1      2014        627
## 50          2      2014        694
## 51          3      2014        841
## 52          4      2014        977
## 53          5      2014       1127
## 54          6      2014       1177
## 55          7      2014       1136
## 56          8      2014       1044
## 57          9      2014        899
## 58         10      2014        959
## 59         11      2014        745
## 60         12      2014        843
par(las=1)
sales_month <- barplot(cum_sales$`data$sales`,
            main = "Numero di vendite cumulate mese per mese (2010-2014)",
            xlab = "Numero mese (1 anno = 12 mesi)",
            ylab = "Vendite",
            col = "white",
            names.arg = seq(1, 60, 1),
            cex.names = 0.7,
            las = 2)

lines(sales_month, cum_sales$`data$sales`, col=2, lwd=1)

Pare che, complessivamente, ci siano picchi di vendite nei periodi estivi (comprensibile…) e che nel 2013 e 2014 le vendite in generale siano maggiori. Ritorneremo su tendenze come queste in una sezione successiva.

Quantifichiamo meglio annualmente:

cum_sales_y <- aggregate(data$sales ~ data$year, data = data, FUN = sum)

print.data.frame(cum_sales_y)
##   data$year data$sales
## 1      2010       8096
## 2      2011       7878
## 3      2012       8935
## 4      2013      10172
## 5      2014      11069

Nel secondo anno si è registrata una leggera riduzione di vendite totali rispetto al primo, mentre negli anni successivi è notevolmente aumentata.

VOLUME

Questa volta proviamo ad utilizzare il pacchetto summarytools, che promette di calcolare pressoché tutti gli indici essenziali con un singolo comando. Eventuali mancanze le colmeremo manualmente.

library(summarytools)
stats_vol <- descr(data$volume, round.digits = 7, transpose = T)

stats_vol
## Descriptive Statistics  
## data$volume  
## N: 240  
## 
##                      Mean      Std.Dev         Min           Q1       Median           Q3
## ------------ ------------ ------------ ----------- ------------ ------------ ------------
##       volume   31.0051875   16.6514472   8.1660000   17.6290000   27.0625000   40.9030000
## 
## Table: Table continues below
## 
##  
## 
##                       Max          MAD          IQR          CV    Skewness   SE.Skewness
## ------------ ------------ ------------ ------------ ----------- ----------- -------------
##       volume   83.5470000   16.1566335   23.2335000   0.5370536   0.8792182     0.1571376
## 
## Table: Table continues below
## 
##  
## 
##                 Kurtosis       N.Valid     Pct.Valid
## ------------ ----------- ------------- -------------
##       volume   0.1505673   240.0000000   100.0000000

Per referenza,

~Mean = media (aritmetica)
Std.Dev = Deviazione standard
Min e Max = Minimo, massimo
Q1 e Q3 = primo e terzo quartile
Median = mediana
IQR = range interquartile
CV= coefficiente di variazione
Skewness = Asimmetria
Kurtosis = Curtosi~

Sembra esserci tutto a parte moda e varianza, che per certi versi son di minor rilievo in scenari come questo dato che avendo tantissimi dati molto precisi (e a volte non interi) e diversi difficilmente ci saranno dati che si ripeteranno più volte e se ci son più valori di moda questi possono essere molto diversi tra loro, mentre alla varianza viene pressoché sempre preferita la deviazione standard (sua radice quadrata). In ogni caso sono facilmente calcolabili.

paste("Moda:", Mode(data$volume))
## [1] "Moda: 14.003" "Moda: 35.335" "Moda: 42.553"
paste(" Varianza:", var(data$volume))
## [1] " Varianza: 277.270692404027"

Volendo, per curiosità scientifica, si potrebbe cercare la moda dei valori interi approssimati:

Mode(round(data$volume))
## [1] 14
## attr(,"freq")
## [1] 11

O anche dei valori approssimati alla prima cifra decimale, anche se a questo punto comincerebbe a crearsi una certa confusione:

Mode(round(data$volume, digits = 1))
## [1] 13.9 16.1 16.2 17.8 24.9 31.7 32.1 40.9 42.6
## attr(,"freq")
## [1] 3

In generale, spesso ha senso creare delle classi e determinare le classi modali, che è un’operazione molto semplice. Per questo studio non lo faremo (tralasciando successive richieste), ci limiteremo ad osservare le densità di distribuzione.

Nell’output abbiamo anche qualche “nuovo” indice come la MAD (Deviazione Mediana Assoluta) che è in realtà un indice piuttosto utile data la sua “robustezza” (bassa o negligibile sensibilità ad outliers o al peso di distribuzioni fortemente non-normali) ma che per questo studio non approfondiremo.

Abbiamo di nuovo un’asimmetria positiva, ma in questo caso una distribuzione leggermente leptocurtica (distribuzione mediamente più “appuntita” rispetto alla normale teorica).

Come per sales, vediamo qualche grafico:

plot(density(data$volume), lwd=2)
abline(v=range(data$volume), col=2)

x <- seq(min(data$volume), max(data$volume), length = 100)
y <- dnorm(x, mean = mean(data$volume), sd = sd(data$volume))
lines(x, y, col = "cyan3", lwd = 1)

#Guadagni cumulativi per città e mese 
cum_volume <- aggregate(data$volume ~ data$month + data$year, data = data, FUN = sum)

print(cum_volume)
##    data$month data$year data$volume
## 1           1      2010      63.751
## 2           2      2010      76.897
## 3           3      2010     111.876
## 4           4      2010     135.196
## 5           5      2010     143.689
## 6           6      2010     138.191
## 7           7      2010     106.764
## 8           8      2010     114.356
## 9           9      2010      88.826
## 10         10      2010      87.164
## 11         11      2010      73.663
## 12         12      2010      92.070
## 13          1      2011      60.658
## 14          2      2011      69.379
## 15          3      2011     101.566
## 16          4      2011     110.267
## 17          5      2011     122.036
## 18          6      2011     138.986
## 19          7      2011     130.855
## 20          8      2011     120.816
## 21          9      2011      95.103
## 22         10      2011      88.274
## 23         11      2011      86.288
## 24         12      2011      83.347
## 25          1      2012      69.791
## 26          2      2012      82.504
## 27          3      2012     109.129
## 28          4      2012     112.128
## 29          5      2012     147.526
## 30          6      2012     149.335
## 31          7      2012     152.463
## 32          8      2012     154.271
## 33          9      2012     112.759
## 34         10      2012     115.004
## 35         11      2012     102.547
## 36         12      2012      97.386
## 37          1      2013      91.739
## 38          2      2013      89.946
## 39          3      2013     125.502
## 40          4      2013     145.621
## 41          5      2013     184.474
## 42          6      2013     184.315
## 43          7      2013     188.551
## 44          8      2013     185.637
## 45          9      2013     136.781
## 46         10      2013     124.589
## 47         11      2013     110.195
## 48         12      2013     119.965
## 49          1      2014      94.076
## 50          2      2014     114.304
## 51          3      2014     139.621
## 52          4      2014     162.877
## 53          5      2014     196.317
## 54          6      2014     215.236
## 55          7      2014     203.805
## 56          8      2014     185.203
## 57          9      2014     158.514
## 58         10      2014     166.541
## 59         11      2014     123.450
## 60         12      2014     149.125
par(las=1)
volume_month <- barplot(cum_volume$`data$volume`,
            main = "Guadagni cumulati mese per mese (2010-2014)",
            xlab = "Numero mese (1 anno = 12 mesi)",
            ylab = "Guadagni (milioni di dollari)",
            col = "white",
            names.arg = seq(1, 60, 1),
            cex.names = 0.7,
            las = 2)

lines(volume_month, cum_volume$`data$volume`, col=2, lwd=1)

Vediamo le somme cumulate annuali:

cum_volume_y <- aggregate(data$volume ~ data$year, data = data, FUN = sum)

print.data.frame(cum_volume_y)
##   data$year data$volume
## 1      2010    1232.443
## 2      2011    1207.575
## 3      2012    1404.843
## 4      2013    1687.315
## 5      2014    1909.069

Stesse tendenze di quelle intraviste nella variabile sales (ha senso, essendo le due strettamente legate e direttamente proporzionali a parità di prezzo medio).

MEDIAN_PRICE

stats_median_p <- descr(data$median_price, round.digits = 7, transpose = T)

stats_median_p
## Descriptive Statistics  
## data$median_price  
## N: 240  
## 
##                                Mean         Std.Dev             Min               Q1
## ------------------ ---------------- --------------- --------------- ----------------
##       median_price   132665.4166667   22662.1486865   73800.0000000   117100.0000000
## 
## Table: Table continues below
## 
##  
## 
##                              Median               Q3              Max             MAD
## ------------------ ---------------- ---------------- ---------------- ---------------
##       median_price   134500.0000000   150100.0000000   180000.0000000   24092.2500000
## 
## Table: Table continues below
## 
##  
## 
##                                IQR          CV     Skewness   SE.Skewness     Kurtosis       N.Valid
## ------------------ --------------- ----------- ------------ ------------- ------------ -------------
##       median_price   32750.0000000   0.1708218   -0.3622768     0.1571376   -0.6427292   240.0000000
## 
## Table: Table continues below
## 
##  
## 
##                        Pct.Valid
## ------------------ -------------
##       median_price   100.0000000
paste("Moda:", Mode(data$median_price))
## [1] "Moda: 130000"
paste(" Varianza:", var(data$median_price))
## [1] " Varianza: 513572983.089261"

Asimmetria e curtosi sono entrambe negative, quindi avremo una distribuzione asimmetrica negativa e platicurtica.

Vediamo:

plot(density(data$median_price), lwd=2)
abline(v=range(data$median_price), col=2)

x <- seq(min(data$median_price), max(data$median_price), length = 100)
y <- dnorm(x, mean = mean(data$median_price), sd = sd(data$median_price))
lines(x, y, col = "cyan3", lwd = 1)

Proviamo a vedere un “andamento medio” del prezzo mediano (sempre tenendo conto di tutte le città in insieme, per avere un’idea generale).

#andamendo del prezzo mediano medio per città e mese 
mean_median_p_m <- aggregate(data$median_price ~ data$month + data$year, data = data, FUN = mean)

print(mean_median_p_m)
##    data$month data$year data$median_price
## 1           1      2010            135450
## 2           2      2010            131100
## 3           3      2010            122325
## 4           4      2010            128300
## 5           5      2010            130125
## 6           6      2010            133275
## 7           7      2010            127225
## 8           8      2010            131100
## 9           9      2010            125825
## 10         10      2010            134150
## 11         11      2010            131600
## 12         12      2010            131825
## 13          1      2011            122250
## 14          2      2011            121200
## 15          3      2011            126750
## 16          4      2011            132475
## 17          5      2011            128775
## 18          6      2011            130350
## 19          7      2011            130750
## 20          8      2011            135350
## 21          9      2011            129100
## 22         10      2011            124100
## 23         11      2011            127900
## 24         12      2011            125250
## 25          1      2012            114250
## 26          2      2012            127225
## 27          3      2012            129275
## 28          4      2012            126250
## 29          5      2012            130575
## 30          6      2012            136175
## 31          7      2012            137375
## 32          8      2012            135050
## 33          9      2012            135275
## 34         10      2012            125075
## 35         11      2012            136750
## 36         12      2012            127650
## 37          1      2013            128325
## 38          2      2013            130600
## 39          3      2013            128225
## 40          4      2013            131050
## 41          5      2013            141025
## 42          6      2013            139525
## 43          7      2013            136200
## 44          8      2013            139525
## 45          9      2013            141350
## 46         10      2013            140550
## 47         11      2013            134150
## 48         12      2013            138150
## 49          1      2014            120975
## 50          2      2014            140250
## 51          3      2014            130500
## 52          4      2014            139375
## 53          5      2014            141925
## 54          6      2014            148775
## 55          7      2014            142200
## 56          8      2014            142350
## 57          9      2014            138650
## 58         10      2014            143525
## 59         11      2014            141125
## 60         12      2014            144125
par(mar = c(5, 5, 4, 0) + 0.1, las=1, cex.axis= 0.95, mgp = c(4, 1, 0))   
#qualche aggiustamento per fare spazio nel grafico
mean_median_p_month <- barplot(mean_median_p_m$`data$median_price`,
            main = "Media del prezzo mediano mese per mese (2010-2014)",
            xlab = "Numero mese (1 anno = 12 mesi)",
            ylab = "Prezzo mediano medio",
            col = "white",
            names.arg = seq(1, 60, 1),
            cex.names = 0.75,
            las = 2)

lines(mean_median_p_month, mean_median_p_m$`data$median_price`, col=2, lwd=1)

Si potrebbe dire che il prezzo mediano è aumentato in modo timido e graduale dal 2010 a fine 2014. Non sembrano esserci evidenti variazioni stagionali.

Proviamo a raccogliere la media in anni:

mean_median_p_y <- aggregate(data$median_price ~ data$year, data = data, FUN = mean)

print.data.frame(mean_median_p_y)
##   data$year data$median_price
## 1      2010          130191.7
## 2      2011          127854.2
## 3      2012          130077.1
## 4      2013          135722.9
## 5      2014          139481.2

Nel secondo anno del periodo di osservazione il prezzo mediano è generalmente diminuito, dopodiché è progressivamente aumentato, proprio come per le vendite (e guadagni).

LISTINGS

stats_listings <- descr(data$listings, round.digits = 7, transpose = T)

stats_listings
## Descriptive Statistics  
## data$listings  
## N: 240  
## 
##                          Mean       Std.Dev           Min             Q1         Median
## -------------- -------------- ------------- ------------- -------------- --------------
##       listings   1738.0208333   752.7077561   743.0000000   1025.0000000   1618.5000000
## 
## Table: Table continues below
## 
##  
## 
##                            Q3            Max           MAD            IQR          CV    Skewness
## -------------- -------------- -------------- ------------- -------------- ----------- -----------
##       listings   2128.0000000   3296.0000000   879.9231000   1029.5000000   0.4330833   0.6454431
## 
## Table: Table continues below
## 
##  
## 
##                  SE.Skewness     Kurtosis       N.Valid     Pct.Valid
## -------------- ------------- ------------ ------------- -------------
##       listings     0.1571376   -0.8101534   240.0000000   100.0000000
paste("Moda:", Mode(data$listings))
## [1] "Moda: 1581"
paste(" Varianza:", var(data$listings))
## [1] " Varianza: 566568.966091353"

Avremo una distribuzione asimmetrica positiva e platicurtica.

plot(density(data$listings), lwd=2)
abline(v=range(data$listings), col=2)

x <- seq(min(data$listings), max(data$listings), length = 100)
y <- dnorm(x, mean = mean(data$listings), sd = sd(data$listings))
lines(x, y, col = "cyan3", lwd = 1)

Fin’ora questa è sicuramente la distribuzione più irregolare.

Mostriamo l’andamento cumulativo degli annunci attivi:

#Quantità cumulata di annunci attivi ogni mese
cum_listings <- aggregate(data$listings ~ data$month + data$year, data = data, FUN = sum)
par(las=1)
listings_month <- barplot(cum_listings$`data$listings`,
            main = "Totale degli annunci attivi mese per mese (2010-2014)",
            xlab = "Numero mese (1 anno = 12 mesi)",
            ylab = "Annunci attivi",
            col = "white",
            ylim = c(0, 8000),
            names.arg = seq(1, 60, 1),
            cex.names = 0.7,
            las = 2)

lines(listings_month, cum_listings$`data$listings`, col=2, lwd=1)

Come per per il numero di vendite e i guadagni, il numero di annunci attivi pare essere maggiore in estate (soprattutto attorno giugno). La differenza relativa tra periodo invernale ed estivo però è considerevolmente minore in questo caso.

Proviamo a raggruppare per anni la somma degli annunci attivati all’anno:

cum_listings_y <- aggregate(data$listings/12 ~ data$year, data = data, FUN = sum)

print.data.frame(cum_listings_y)
##   data$year data$listings/12
## 1      2010         7304.000
## 2      2011         7398.583
## 3      2012         7107.250
## 4      2013         6710.417
## 5      2014         6240.167

(Per formalità listings dovrebbe essere arrotondato, essendo una variabile intera discreta, ma qui ignoriamo essendo questo un semplice confronto immediato.)

Nel secondo anno il numero totale di annunci attivi ogni mese è stato generalmente maggiore del primo, ma negli anni successivi è costantemente diminuito.

Ad occhio, con i dati così raggruppati, sembrerebbe esserci una correlazione inversa con il prezzo mediano e le precedenti variabili, che avevamo visto diminuire nel secondo anno e aumentare successivamente. Proviamo a vedere quanto è forte la correlazione col prezzo mediano con alcune funzioni.

Nota per la correzione: su questa breve porzione della correlazione non sono per niente sicuro… c’è anche da dire che è un po’ superfluo.
medianP_listings_correlation_y <- c(cor(mean_median_p_y$`data$median_price`, cum_listings_y$`data$listings`), cor(mean_median_p_y$`data$median_price`, cum_listings_y$`data$listings`, method = "kendall"), cor(mean_median_p_y$`data$median_price`, cum_listings_y$`data$listings`, method = "spearman"))


print(medianP_listings_correlation_y)
## [1] -0.9847591 -0.8000000 -0.9000000

Senza entrare particolarmente nei dettagli, il coefficiente di correlazione è un indice che può assumere un valore compreso nel range continuo -1 e 1. Un coefficiente pari a 1 indica una correlazione positiva perfetta (le due variabili aumentano assieme, un coefficiente di -1 indica una correlazione negativa perfetta (se una variabile aumenta, l’altra diminuisce), mentre 0 indicherebbe l’assenza di correlazione. I metodi utilizzati calcolano in modo diverso la correlazione, e hanno ognuno un proprio campo ideale di applicazione, ma per oggi non ci interessa.

Qui il coefficiente è abbastanza vicino a -1 con tutti i metodi, quindi c’è effettivamente una certa correlazione se si considerano le medie annuali.

Vediamo mese per mese:

medianP_listings_correlation_m <- c(cor(mean_median_p_m$`data$median_price`, cum_listings$`data$listings`), cor(mean_median_p_m$`data$median_price`, cum_listings$`data$listings`, method = "kendall"), cor(mean_median_p_m$`data$median_price`, cum_listings$`data$listings`, method = "spearman"))


print(medianP_listings_correlation_m)
## [1] -0.3937336 -0.2619520 -0.3880514

Mensilmente, la correlazione negativa sembra molto debole.

Facciamo la stessa cosa col numero di vendite.

sales_listings_correlation_y <- c(cor(cum_sales_y$`data$sales`, cum_listings_y$`data$listings`), cor(cum_sales_y$`data$sales`, cum_listings_y$`data$listings`, method = "kendall"), cor(cum_sales_y$`data$sales`, cum_listings_y$`data$listings`, method = "spearman"))


print(sales_listings_correlation_y)
## [1] -0.9902385 -1.0000000 -1.0000000

Apparentemente qui la correlazione inversa annuale è pressoché perfetta.

Vediamo mensilmente:

sales_listings_correlation_m <- c(cor(cum_sales$`data$sales`, cum_listings$`data$listings`), cor(cum_sales$`data$sales`, cum_listings$`data$listings`, method = "kendall"), cor(cum_sales$`data$sales`, cum_listings$`data$listings`, method = "spearman"))


print(sales_listings_correlation_m)
## [1] -0.040204981  0.001697313 -0.001111497

Mentre osservando a livello mensile la correlazione pare nulla.

~Immagino che gli andamenti mensili di alcune variabili potrebbero essere “sfasati” tra loro e ciò mascheri la presenza di un forte collegamento, 
oppure osservare le medie annuali comporti una qualche “distorsione” o bias. 
Purtroppo ora come ora non ho le competenze per analizzare la prima ipotesi.~

(Non sono sicuro dell’utilità e accuratezza di queste osservazioni.)

Fine porzione incerta

MONTHS_INVENTORY

stats_months_inv <- descr(data$months_inventory, round.digits = 7, transpose = T)

stats_months_inv
## Descriptive Statistics  
## data$months_inventory  
## N: 240  
## 
##                               Mean     Std.Dev         Min          Q1      Median           Q3
## ---------------------- ----------- ----------- ----------- ----------- ----------- ------------
##       months_inventory   9.1925000   2.3036686   3.4000000   7.8000000   8.9500000   11.0000000
## 
## Table: Table continues below
## 
##  
## 
##                                 Max         MAD         IQR          CV    Skewness   SE.Skewness
## ---------------------- ------------ ----------- ----------- ----------- ----------- -------------
##       months_inventory   14.9000000   2.1497700   3.1500000   0.2506031   0.0407194     0.1571376
## 
## Table: Table continues below
## 
##  
## 
##                            Kurtosis       N.Valid     Pct.Valid
## ---------------------- ------------ ------------- -------------
##       months_inventory   -0.1979448   240.0000000   100.0000000

Distribuzione asimmetrica positiva (ma non troppo distante dalla distribuzione normale, data la bassa skewness) e leggermente platicurtica.

paste("Moda:", Mode(data$months_inventory))
## [1] "Moda: 8.1"
paste(" Varianza:", var(data$months_inventory))
## [1] " Varianza: 5.30688912133891"
plot(density(data$months_inventory), lwd=2)
abline(v=range(data$months_inventory), col=2)

x <- seq(min(data$months_inventory), max(data$months_inventory), length = 100)
y <- dnorm(x, mean = mean(data$months_inventory), sd = sd(data$months_inventory))
lines(x, y, col = "cyan3", lwd = 1)

Facciamo la media mensile (sempre considerando le 4 città contemporaneamente):

months_inventory_m <- aggregate(data$months_inventory ~ data$month + data$year, data = data, FUN = mean)

print(months_inventory_m)
##    data$month data$year data$months_inventory
## 1           1      2010                 8.750
## 2           2      2010                 9.150
## 3           3      2010                 9.475
## 4           4      2010                10.025
## 5           5      2010                 9.700
## 6           6      2010                10.050
## 7           7      2010                10.550
## 8           8      2010                10.625
## 9           9      2010                10.650
## 10         10      2010                10.500
## 11         11      2010                10.400
## 12         12      2010                 9.800
## 13          1      2011                 9.950
## 14          2      2011                10.200
## 15          3      2011                10.900
## 16          4      2011                11.525
## 17          5      2011                12.075
## 18          6      2011                11.925
## 19          7      2011                11.575
## 20          8      2011                11.225
## 21          9      2011                11.025
## 22         10      2011                10.650
## 23         11      2011                10.175
## 24         12      2011                 9.625
## 25          1      2012                 9.950
## 26          2      2012                10.150
## 27          3      2012                10.475
## 28          4      2012                10.625
## 29          5      2012                10.525
## 30          6      2012                10.475
## 31          7      2012                10.350
## 32          8      2012                 9.800
## 33          9      2012                 9.625
## 34         10      2012                 9.225
## 35         11      2012                 8.975
## 36         12      2012                 8.375
## 37          1      2013                 8.550
## 38          2      2013                 8.700
## 39          3      2013                 8.900
## 40          4      2013                 9.050
## 41          5      2013                 8.700
## 42          6      2013                 8.600
## 43          7      2013                 8.225
## 44          8      2013                 8.000
## 45          9      2013                 7.650
## 46         10      2013                 7.500
## 47         11      2013                 7.225
## 48         12      2013                 6.700
## 49          1      2014                 7.000
## 50          2      2014                 7.100
## 51          3      2014                 7.275
## 52          4      2014                 7.375
## 53          5      2014                 7.375
## 54          6      2014                 7.475
## 55          7      2014                 7.375
## 56          8      2014                 7.300
## 57          9      2014                 6.975
## 58         10      2014                 6.825
## 59         11      2014                 6.525
## 60         12      2014                 6.075
par(las=1)
months_inventory_month <- barplot(months_inventory_m$`data$months_inventory`,
            main = "Andamento medio delle stime del tempo per la vendita degli annunci attivi (2010-2014)",
            xlab = "Numero mese (1 anno = 12 mesi)",
            ylab = "Tempo (mesi)",
            col = "white",
            names.arg = seq(1, 60, 1),
            cex.names = 0.7,
            las = 2)

lines(months_inventory_month, months_inventory_m$`data$months_inventory`, col=2, lwd=1)

months_inventory_y <- aggregate(data$months_inventory ~ data$year, data = data, FUN = mean)

print(months_inventory_y)
##   data$year data$months_inventory
## 1      2010              9.972917
## 2      2011             10.904167
## 3      2012              9.879167
## 4      2013              8.150000
## 5      2014              7.056250

Come per Listings, c’è un aumento nel secondo anno e una diminuzione costante negli anni successivi. Intuitivamente ha senso: nel secondo anno, ad esempio, listings è maggiore che negli altri anni (più annunci attivi) e al contempo sales è al minimo (minor numero di vendite). Logicamente, la stima del tempo teoricamente necessario per vendere tutti gli annunci sarà più alta.

-

3: Variabili con maggiore variabilità e con maggior asimmetria

Per identificare la variabile con maggiore variabilità, può essere opportuno confrontare i coefficienti di variazione (CV).

(Le prime 3 variabili, quelle trattate come qualitative, possono essere tranquillamente escluse da questo “esercizio” in funzione della loro equimodalità e del basso numero di modalità.

CVs <- c(sales_CV, stats_vol$CV, stats_median_p$CV, stats_listings$CV, stats_months_inv$CV)
#per sales avevamo calcolato manualmente il CV
colors <- c("green2", "cyan3", "orange2", "pink3", "cyan4")
par(las=1)
barplot(CVs, names.arg = c("Sales","Volume","Median Price","Listings","Months Inv."), col = colors, 
        ylim = c(0, 0.6),
        main = "Confronto della Variabilità Relativa (CV)",
        xlab = "Variabile quantitativa", ylab = "Coefficiente di Variazione")

Appare evidente che Volume sia la variabile più… variabile.
In ordine decrescente dunque abbiamo:

  1. Volume

  2. Listings

  3. Sales

  4. Months_Inventory

  5. Median_Price

Un’operazione analoga può essere eseguita per visualizzare la variabile più asimmetrica.

skew_vals <- c(sales_skew, stats_vol$Skewness, stats_median_p$Skewness, stats_listings$Skewness, stats_months_inv$Skewness)
par(las=1)
barplot(skew_vals, names.arg = c("Sales","Volume","Median Price","Listings","Months Inv."), col = colors,
        ylim = c(-0.4, 1),
        main = "Confronto dell'asimmetria (skewness)",
        xlab = "Variabile quantitativa", ylab = "Asimmetria")

abline(h = 0, col = "black", lwd = 2) #un semplice asse orizzontale

Anche qui Volume è la più asimmetrica, quantitativamente. Abbiamo un’unica variabile asimmetrica negativamente ovvero median price, che in senso assoluto è la quarta più asimmetrica.

In ordine decrescente (assoluto):

  1. Volume

  2. Sales

  3. Listings

  4. Median_Price*

  5. Months_Inventory

4: Divisione di una delle variabili quantitative in classi con distribuzione di frequenza, grafico a barre e indice di Gini

Proviamo con la variabile median price.

#Dividiamo in 9 classi di "larghezza" 12000
medianP_cl <- cut(data$median_price, c(72000, 84000, 96000, 108000, 120000, 132000, 144000, 156000, 168000, 180000), dig.lab = 6)
medianP_cl_dist <- as.data.frame(dist_tab(medianP_cl))

print(medianP_cl_dist)
##   dataframe.interval dataframe.freq dataframe.cum.freq dataframe.percent dataframe.cum.percent
## 1      (72000,84000]              2                  2              0.83                  0.83
## 2      (84000,96000]             17                 19              7.08                  7.92
## 3     (96000,108000]             25                 44             10.42                 18.33
## 4    (108000,120000]             22                 66              9.17                 27.50
## 5    (120000,132000]             36                102             15.00                 42.50
## 6    (132000,144000]             50                152             20.83                 63.33
## 7    (144000,156000]             55                207             22.92                 86.25
## 8    (156000,168000]             23                230              9.58                 95.83
## 9    (168000,180000]             10                240              4.17                100.00

~(Nota: per convenzione, il limite superiore di una classe è escluso, il limite inferiore è incluso. In questo caso però è l’opposto.
Ai fini pratici non cambia nulla in questo specifico scenario, ma è giusto specificarlo. 
Per creare classi con il criterio convenzionale è sufficiente aggiungere il parametro right = F nella funzione cut)~

par(mar = c(5, 3.7, 3.5, 0) + 0, las=1, cex.axis= 1, mgp = c(2.5, 1, 0))
barplot(medianP_cl_dist$dataframe.freq, names.arg = medianP_cl_dist$dataframe.interval,
        main = "Distribuzione di frequenza prezzo mediano (classi)",
        cex.names = 0.7,
        ylim = c(0, 60),
        xlab = "Classi", ylab = "Frequenza (assoluta)")

Passiamo all’indice di Gini, che ricordiamo di base misura la “propensione” di una variabile qualitativa ad assumere le sue diverse modalità, andando a considerare la sua distribuzione di frequenze.

Nota: Altra parte un po’ incerta… mi sembra di capire che c’è una differenza formale sostanziale in come si calcola ed esprime l’indice quando si ha a che fare con classi/variabili quantitative piuttosto che semplici variabili qualitative, in virtù di valori e “pesi” differenti. Probabilmente è un campo un po’ differente. Io considererò le classi come semplici “modalità” qualitative, come lo sarebbero dei colori o dei nomi.

Creiamo una semplice funzione per il calcolo dell’indice (compresa normalizzazione):

gini.index<- function(x){
  ni=table(x)
  fi=ni/length(x)
  fi2=fi^2
  J = length(table(x))
  
  gini = 1-sum(fi2)
  #normalizzazione
  gini.norm = gini/((J-1)/J)
  
  return(gini.norm)
}
gini.index(medianP_cl)
## [1] 0.9521094

Se i calcoli sono corretti, abbiamo un’eterogeneità del 95,2%.

Fine parte incerta

5: “Indovinare” l’indice di Gini per la variabile city.

Riguardiamo un momento la distribuzione per city:

print.data.frame(city_distrib)
##                       city_ni city_fi city_Ni city_Fi
## Beaumont                   60    0.25      60    0.25
## Bryan-College Station      60    0.25     120    0.50
## Tyler                      60    0.25     180    0.75
## Wichita Falls              60    0.25     240    1.00

Avevamo individuato una distribuzione equimodale (quadrimodale), dunque un’equidistribuzione, il ché significa massima eterogeneità, e quindi indice di Gini pari a 1.

Confermiamo, per formalità:

gini.index(data$city)
## [1] 1

6: Probabilità di ottenere specifiche modalità nel leggere una riga “a caso”.

Iniziamo dalla probabilità di prendere una riga che riporti “Beaumont” (città).

Già intutivamente potremmo dire che sarà del 25% perché sappiamo che la distribuzione di city è quadrimodale, e ci sono 60 righe per ogni città, per un totale di 240. 60/240 = 0.25 -> 25%.

Verifichiamo come se non lo sapessimo:

Beaum_freq <- sum(data$city == "Beaumont")
Beaum_prob <- Beaum_freq/nrow(data)   #casi favorevoli/casi possibili

paste0(Beaum_prob*100,"%")
## [1] "25%"

Passiamo alla probabilità di avere “luglio” (che sarà indicato come 7 nel dataset per la variabile month).

Jul_freq <- sum(data$month == "7")   #sarà 20 per ogni mese
Jul_prob <- Jul_freq/nrow(data)

paste0(Jul_prob*100,"%")
## [1] "8.33333333333333%"

(Una possibilità su 12)

Infine vediamo la probabilità di avere “dicembre 2012”:

Dec_2012_freq <- sum(data$month == "12" & data$year == "2012") 
#4, uno per città

Dec_2012_prob <- Dec_2012_freq/nrow(data)

paste0(Dec_2012_prob*100, "%")
## [1] "1.66666666666667%"

(Una possibilità su 60)

7: Aggiunta di una colonna col prezzo medio tramite le altre variabili a disposizione

In questo scenario, il prezzo medio è facilmente ricavabile dal rapporto tra i guadagni (volume) e il numero di vendite (sales).
Lo ricaviamo per ogni riga e creeremo una colonna con i risultati.

data$mean_price <- (data$volume*1000000)/data$sales

Verifichiamo se il dataset è stato correttamente aggiornato:

head(data, 10)
##        city year month sales volume median_price listings months_inventory mean_price
## 1  Beaumont 2010     1    83 14.162       163800     1533              9.5   170626.5
## 2  Beaumont 2010     2   108 17.690       138200     1586             10.0   163796.3
## 3  Beaumont 2010     3   182 28.701       122400     1689             10.6   157697.8
## 4  Beaumont 2010     4   200 26.819       123200     1708             10.6   134095.0
## 5  Beaumont 2010     5   202 28.833       123100     1771             10.9   142737.6
## 6  Beaumont 2010     6   189 27.219       122800     1803             11.1   144015.9
## 7  Beaumont 2010     7   164 22.706       124300     1857             11.7   138451.2
## 8  Beaumont 2010     8   174 25.237       136800     1830             11.6   145040.2
## 9  Beaumont 2010     9   124 17.233       121100     1829             11.7   138975.8
## 10 Beaumont 2010    10   150 23.904       138500     1779             11.5   159360.0
#***Questo blocco non verrà eseguito***
#Sarebbe possibile scegliere una posizione specifica per la colonna
#Ad esempio, se avessimo voluto metterla dopo median_price:

pos <- which(names(data) == "median_price") + 1


data <- cbind(dat[, 1:pos], mean_price = data$mean_price, data[, (pos+1):ncol(data)])

#***Fine blocco non eseguibile***

8: Aggiunta di una colonna che dia un’idea di “efficacia” degli annunci di vendita

La quantità di dati disponibile è in realtà relativamente limitata e c’è un certo margine di “soggettività” nella definizione di efficacia, ma si possono comunque fare delle considerazioni.
La più ovvia di gran lunga sarebbe quella di osservare il rapporto tra vendite e annunci attivi: un rapporto maggiore implicherebbe una maggiore “efficacia”.

Un parametro minore da considerare potrebbe essere il rapporto tra prezzo medio e prezzo mediano. Un prezzo medio maggiore può essere interpretato come un maggior successo delle inserzioni e una maggiore capacità di generare profitto con gli annunci più di “fascia alta” o comunque di rendere “attraenti” gli immobili, spingendo a valutazioni maggiori.

La variazione stagionale (col picco estivo e il fondo invernale) è probabilmente dovuta al semplice fatto che le persone siano più invogliate/facilitate a fare questo tipo di grandi acquisti in determinati periodi dell’anno, soprattutto considerando che la tendenza è perfettamente ciclica. Difficile immaginare abbia a che fare con un’effettiva “efficacia” degli annunci, tantomeno con variazioni così ampie.

Altre metriche sarebbero in qualche modo variazioni delle precedenti. Ad esempio, potremmo considerare il tempo stimato di vendita degli annunci attivi (months_inventory), ma ciò sarebbe comunque correlato alle vendite e agli annunci attivi (anche se non calcolato in modo diretto da un rapporto vendite/annunci.)

I guadagni per annuncio attivo nemmeno fornirebbero informazioni che non sarebbero in qualche modo ridondanti (o persino fuorvianti).

Potremmo quindi ricavare una combinazione ponderata delle prime due metriche citate, dando alla seconda un “peso” minore, e riscalare secondo alcuni fattori. Faremo in modo che ad un indice di 1 corrisponda (idealmente) un rapporto vendite/annunci “tipico” (per esempio il rapporto tra la mediana di sales e la mediana di listings) e a un rapporto prezzo medio/prezzo mediano pari a 1.

sales_to_list_ratio <- (data$sales)/(data$listings)
mean_to_median_ratio <- (data$mean_price)/(data$median_price)

sales_to_list_typ <- median(data$sales)/median(data$listings)

#normalizziamo rispetto a 1
sales_to_list_norm <- (sales_to_list_ratio)/(sales_to_list_typ)

#Definiamo un indice AEI (Ad Efficacy Index) e creiamo la colonna
data$AEI <- (0.8*sales_to_list_norm)+(0.2*mean_to_median_ratio)

-

Vediamo rapidamente le caratteristiche di questa nuova variabile:

stats_AEI <- descr(data$AEI, round.digits = 7, transpose = T)

stats_AEI
## Descriptive Statistics  
## data$AEI  
## N: 240  
## 
##                  Mean     Std.Dev         Min          Q1      Median          Q3         Max
## --------- ----------- ----------- ----------- ----------- ----------- ----------- -----------
##       AEI   1.1090300   0.3473404   0.6077843   0.8905339   1.0454133   1.2350649   3.0963658
## 
## Table: Table continues below
## 
##  
## 
##                   MAD         IQR          CV    Skewness   SE.Skewness    Kurtosis       N.Valid
## --------- ----------- ----------- ----------- ----------- ------------- ----------- -------------
##       AEI   0.2500839   0.3396027   0.3131930   2.0813694     0.1571376   6.8318433   240.0000000
## 
## Table: Table continues below
## 
##  
## 
##               Pct.Valid
## --------- -------------
##       AEI   100.0000000
plot(density(data$AEI), lwd=2)
abline(v=range(data$AEI), col=2)

x <- seq(min(data$AEI), max(data$AEI), length = 100)
y <- dnorm(x, mean = mean(data$AEI), sd = sd(data$AEI))
lines(x, y, col = "cyan3", lwd = 1)

-

Infine diamo un’occhiata all’andamento di questa nuova variabile per ogni città separatamente, in un unico grafico. Essendo un’operazione un po’ meno semplice ci serviremo della libreria ggplot2 (che servirà anche per i quesiti della prossima macro-sezione).

library(ggplot2)
ggplot(data, aes(x = (as.Date(paste(year, month, "01", sep = "-"))), y = AEI, color = city)) +
    geom_line() +
    geom_point() +
    scale_color_brewer(palette = "Set1") +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
    scale_y_continuous(limits = c(0, 3.5))+
    labs(title = "Andamento dell'AEI per Città", x = "Data", y = "AEI", color = "Città") +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
          axis.title = element_text(size = 14),
          legend.position = "bottom")

Vediamo anche con 4 grafici separati:

ggplot(data, aes(x = (as.Date(paste(year, month, "01", sep = "-"))), y = AEI, color = city)) +
    geom_line() +
    geom_point() +
    scale_color_brewer(palette = "Set1") +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
    scale_y_continuous(limits = c(0, 3.5))+
    labs(title = "Andamento dell'AEI per Città", x = "Data", y = "AEI", color = "Città") +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
          axis.title = element_text(size = 14),
          legend.position = "bottom") +
    facet_wrap(~ city, ncol = 2)

Pare che l’andamento in Bryan-College Station sia alquanto irregolare se confrontato con le altre città, con picchi estivi assai più evidenti. Al contrario, sembra esserci la maggiore costanza per Tyler.
Probabilmente, una buona porzione degli andamenti stagionali che abbiamo notato in precedenza per variabili come sales, volume, listings, etc etc… è dovuta al mercato di Bryan-College Station, soprattutto negli ultimi 2 anni.


Vediamo l’indice anche numericamente e graficamente con un confronto medio:

per_city_AEI <- aggregate(data$AEI ~ data$city, data = data, FUN = mean)

print.data.frame(per_city_AEI)
##               data$city  data$AEI
## 1              Beaumont 1.0091317
## 2 Bryan-College Station 1.3199782
## 3                 Tyler 0.9269012
## 4         Wichita Falls 1.1801089
col = c("red3", "blue3", "green4", "purple3" )

par(las=1)
barplot(per_city_AEI$`data$AEI`, names.arg = per_city_AEI$`data$city`, col = col,
        ylim = c(0, 1.5),
        main = "Indice di efficacia medio per città nel periodo di riferimento",
        xlab = "Città", ylab = "AEI")

abline(h = 0, col = "black", lwd = 2)

Bryan-College è la città dove, in media, l’indice che abbiamo calcolato risulta più alto, Tyler quella dove risulta più basso.

-

9: Summary per alcune variabili condizionatamente a città, anno, mese

Nota: Da non confondere con alcune operazioni (vagamente) simili che ho eseguito nel task 2

(Commento personale secondario: il cheatsheet di dplyr è oggettivamente utile ma personalmente non sono riuscito a trovarlo super intuitivo/esplicativo.
Mi sono trovato maggiormente ad utilizzare altri manuali/guide.)

Iniziamo creando un sommario di alcuni indici principali (in ordine: minimo, massimo, deviazione standard, media, mediana) di sales, volume, median price e listings, raggruppando città e anno:

summary_city_year <- data %>%
    group_by(city, year) %>%
    summarise(
        Sales_min = min(sales),
        Sales_max = max(sales),
        Sales_sd = sd(sales),
        Sales_mean = round(mean(sales)),
        Sales_mdn = round(median(sales)),
        Vol_min = min(volume),
        Vol_max = max(volume),
        Vol_sd = sd(volume),
        Vol_mean = mean(volume),
        Vol_mdn = median(volume),
        MdnP_min = min(median_price),
        MdnP_max = max(median_price),
        MdnP_sd = sd(median_price),
        MdnP_mean = round(mean(median_price)),
        MdnP_mdn = round(median(median_price)),
        List_min = min(listings),
        List_max = max(listings),
        List_sd = sd(listings),
        List_mean = round(mean(listings)),
        List_mdn = round(median(listings))
    )

print.data.frame(summary_city_year, digits = 4)
##                     city year Sales_min Sales_max Sales_sd Sales_mean Sales_mdn Vol_min Vol_max Vol_sd Vol_mean Vol_mdn MdnP_min MdnP_max MdnP_sd MdnP_mean MdnP_mdn List_min List_max List_sd List_mean List_mdn
## 1               Beaumont 2010        83       202    36.92        156       157  14.162   28.83  4.954    22.65   23.30   121100   163800   13354    133117   128400     1533     1857  101.91      1731     1756
## 2               Beaumont 2011       108       177    22.66        144       148  15.489   28.48  4.301    21.10   21.38   111100   144600    9603    125642   126900     1596     1845   74.67      1748     1764
## 3               Beaumont 2012       110       218    28.39        172       176  13.496   31.15  4.922    24.47   25.32   110000   134700    7973    126533   128300     1570     1765   54.56      1691     1694
## 4               Beaumont 2013       140       273    37.73        201       202  20.342   42.03  6.437    30.31   30.13   121100   147000    7785    132400   132750     1534     1708   54.92      1640     1657
## 5               Beaumont 2014       148       262    36.49        214       210  18.077   41.19  7.050    32.13   33.26   106700   142400    9835    132250   134150     1500     1672   57.37      1587     1590
## 6  Bryan-College Station 2010        89       286    70.75        168       153  15.151   47.45 10.818    28.73   26.29   148300   165300    5474    153533   151800     1298     1984  165.16      1562     1563
## 7  Bryan-College Station 2011        94       284    62.19        167       148  15.242   47.78 10.313    28.93   25.88   146700   157300    3709    151417   150300     1362     1840  156.09      1606     1572
## 8  Bryan-College Station 2012       115       296    74.28        197       161  19.789   55.45 13.490    35.36   30.07   140700   170000    7096    153567   153700     1442     1834  152.98      1610     1576
## 9  Bryan-College Station 2013       125       402    95.85        238       188  18.977   76.12 19.541    45.12   35.19   146900   167300    5429    159392   159900     1057     1750  228.36      1406     1424
## 10 Bryan-College Station 2014       152       403    86.69        260       246  29.457   83.55 17.971    52.81   48.44   155300   180000    7776    169533   170900      882     1271  127.44      1106     1096
## 11                 Tyler 2010       155       316    48.98        228       229  24.411   49.91  8.394    36.35   34.27   129400   143100    4782    135175   134450     2727     3296  226.43      3051     3099
## 12                 Tyler 2011       143       313    49.62        239       247  21.050   52.32  9.406    38.55   39.77   120600   152600    8505    136217   134650     2720     3266  184.60      3070     3098
## 13                 Tyler 2012       169       322    46.40        264       276  25.386   57.39 10.230    44.01   45.36   124200   152100    7983    139250   138800     2633     3072  122.95      2910     2918
## 14                 Tyler 2013       197       369    53.05        287       288  32.083   63.05 10.326    50.32   51.35   132400   155600    6726    146100   146100     2500     2998  159.72      2824     2884
## 15                 Tyler 2014       238       423    56.85        332       340  36.916   80.81 12.761    59.60   63.17   130700   161600    8543    150467   151900     2272     2875  172.17      2670     2716
## 16         Wichita Falls 2010        89       167    26.62        123       123   8.951   20.88  4.066    14.97   15.01    86400   119200   10361     98942    98350      904     1028   45.17       959      957
## 17         Wichita Falls 2011        79       135    19.76        106       111   8.166   15.28  2.515    12.05   13.21    73800   113600   10632     98142   101500      844     1052   58.59       975      982
## 18         Wichita Falls 2012        90       132    14.25        112       116   9.695   17.79  2.662    13.23   12.78    82100   118800   12347    100958    99900      801      941   41.22       896      907
## 19         Wichita Falls 2013        79       159    26.00        121       122   9.666   19.06  3.106    14.85   15.01    85900   121300   10383    105000   103600      743      923   52.91       841      847
## 20         Wichita Falls 2014        89       150    21.09        117       111   9.626   18.67  3.131    14.54   13.90    90000   135300   12444    105675   104250      746      973   73.05       877      884

-

A questo punto potremmo visualizzare ogni elemento del sommario in un grafico (in questo caso usiamo quello a barre), come ad esempio la media delle vendite:

ggplot(summary_city_year, aes(x=year, y=Sales_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media vendite per città ogni anno", x= "Anno", y="Vendite medie", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_continuous(limits = c(0, 350)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")

-

O la mediana del numero di annunci attivi:

ggplot(summary_city_year, aes(x=year, y=List_mdn, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Mediana del numero di annunci attivi per città ogni anno", x= "Anno", y="Mediana annunci attivi", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_continuous(limits = c(0, 3100)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")

-

O la media del prezzo mediano:

ggplot(summary_city_year, aes(x=year, y=MdnP_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media del prezzo mediano per città ogni anno", x= "Anno", y="Media prezzo mediano", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_continuous(limits = c(0, 200000)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")

-

O la deviazione standard degli incassi:

ggplot(summary_city_year, aes(x=year, y=Vol_sd, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Deviazione standard degli incassi per città ogni anno", x= "Anno", y="Deviazione standard incassi", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")

-

Possiamo provare a fare la stessa cosa ma raggruppando città e mese.

summary_city_month <- data %>%
    group_by(city, month) %>%
    summarise(
        Sales_min = min(sales),
        Sales_max = max(sales),
        Sales_sd = sd(sales),
        Sales_mean = round(mean(sales)),
        Sales_mdn = round(median(sales)),
        Vol_min = min(volume),
        Vol_max = max(volume),
        Vol_sd = sd(volume),
        Vol_mean = mean(volume),
        Vol_mdn = median(volume),
        MdnP_min = min(median_price),
        MdnP_max = max(median_price),
        MdnP_sd = sd(median_price),
        MdnP_mean = round(mean(median_price)),
        MdnP_mdn = round(median(median_price)),
        List_min = min(listings),
        List_max = max(listings),
        List_sd = sd(listings),
        List_mean = round(mean(listings)),
        List_mdn = round(median(listings))
    )

-

Vediamo di nuovo qualche grafico. Media delle vendite:

ggplot(summary_city_month, aes(x=month, y=Sales_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media vendite per città nei singoli mesi", x= "Mese", y="Vendite medie", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12, labels = month.name) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")

Proviamo anche a visualizzare separatamente per le 4 città:

ggplot(summary_city_month, aes(x=month, y=Sales_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media vendite per città nei singoli mesi", x= "Mese", y="Vendite medie", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")+
  facet_wrap(~ city, ncol=2)

-

Vediamo anche la media degli annunci per mese, direttamente separata per città:

ggplot(summary_city_month, aes(x=month, y=List_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media del numero di annunci attivi per città nei singoli mesi", x= "Anno", y="Media annunci attivi", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = 1:12) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")+
  facet_wrap(~ city, ncol=2)

Da notare che per Bryan-College Station sembra sia la primavara il periodo con più annunci attivi, in controtendenza con le altre città dove è l’estate ad avere il picco.

Vediamo infine con la media del prezzo mediano, sempre già separata per città:

ggplot(summary_city_month, aes(x=month, y=MdnP_mean, fill=city)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Media del prezzo mediano per città nei singoli mesi", x= "Mese", y="Media prezzo mediano", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 12),
        legend.position = "bottom")+
  facet_wrap(~ city, ncol=2)

Potremmo continuare per eternità con tante combinazioni, ma fermiamoci qui.

-

Ancora più grafici con ggplot2

1: Distribuzione del prezzo mediano delle case tra le varie città con i boxplot.

Creiamo subito i boxplot con ggplot:

ggplot(data, aes(x=city, y=median_price))+
    geom_boxplot(fill="lightblue1",
                 outlier.color = "red",
                 outlier.size = 3)+
    geom_jitter(width = 0.27, size = 1.1, 
                aes(color = median_price))+
    scale_y_continuous(limits = c(0, NA))+
    stat_summary(fun = mean, geom = "point", shape = 16, 
                 size = 2.5, colour = "orange")+
    labs(x= "Città", y="Prezzo mediano")

Nota: outlier(s) indicati con puntino rosso, jitter points in gradazioni di blu aggiunti per avere ulteriore rappresentazione della distribuzione, punto di media indicato con puntino arancione.

Beaumont, Bryan-College Station e Wichita Falls presentano un singolo report con valore considerato outlier, con un dato del prezzo mediano significativamente superiore a tutti gli altri report. Tyler non sembra presentare outliers.

Appare evidente che Bryan-College Station abbia i prezzi più alti di tutte le città, seguita da Tyler, Beaumont e infine Wichita Falls con i prezzi di gran lunga più bassi, come avevamo già intravisto nelle sezioni precedenti.

Wichita Falls in generale presenta relativamente poche vendite, pochi annunci attivi e prezzi bassi.
Abbiamo però visto un AEI in seconda posizione (sezione 8 del “capitolo” precedente), dunque non c’è necessariamente motivo di immaginare la presenza di una qualche “crisi”, più realisticamente si tratta semplicemente di un mercato più stazionario e in “piccola scala”.

Bryan-College Station sembra avere il mercato “più attraente”, in ogni caso. Pur non presentando il maggior numero di annunci attivi o di vendite assolute, le vendite percentuali sono le maggiori (Tyler ad esempio porta molti più annunci a fronte di un numero di vendite non di molto superiore) e i prezzi sono i più alti. L’indice AEI che avevamo calcolato era il più alto.
Buona parte di questi risultati però è probabilmente dovuta agli ultimi 2 anni di misurazione, con un rapporto vendite/annunci attivi notevolmente elevato rispetto agli altri anni e in generale alle altre città.
In particolare si è potuto notare che è il numero di annunci attivi ad essere calato in modo più evidente specificamente per questa città, a differenza del numero di vendite che è aumentato in modo abbastanza coerente tra le città (escludendo Whichita Falls, che come già accennato è rimasta prossoché stazionaria).

In termini di guadagno assoluto, comunque, Tyler è al primo posto. Lo vedremo meglio nella prossima sezione.

-

2: Distribuzione del valore totale delle vendite tra le varie città e i vari anni con i boxplot (o varianti)

Proviamo a rappresentare con qualche grafico “ibrido”. In particolare, utilizziamo una combinazione di boxplot azzurro e violinplot dimezzato rosa, con l’aggiunta di jitter points arancioni e un trattino di media viola. Per il momento limitiamoci alle città:

library(PupillometryR)
#utilizziamo questo pacchetto per aggiungere una piccola funzione specifica
#geom_flat_violin
ggplot(data, aes(x=city, y=volume))+
    geom_flat_violin(position = position_nudge(x = .2, y = 0), 
                     adjust = .5, trim = T, fill="pink2") +
    geom_jitter(width = 0.18, size = 1.15, aes(color = volume)) +
    scale_color_gradient(low = "orange4", high = "orange")+
    geom_boxplot(width = .15, outlier.color = "red", fill="lightblue2")+
    stat_summary(fun = mean, geom = "crossbar", 
                 width=0.06, colour = "purple")+
    labs(title = "Distribuzione dei guadagni per città",
         x= "Città", y="Guadagni (milioni di dollari)")+
    theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 13))+
    scale_y_continuous(limits = c(0, NA))
## Warning: Using the `size` aesthetic with geom_polygon was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Proviamo ora ad integrare gli anni.

ggplot(data, aes(x=city, y=volume))+
    geom_flat_violin(position = position_nudge(x = .2, y = 0), 
                     adjust = .5, trim = T, fill="pink2") +
    geom_jitter(width = 0.18, size = 1.15, aes(color = volume)) +
    scale_color_gradient(low = "orange4", high = "orange")+
    geom_boxplot(width = .15, outlier.color = "red", fill="lightblue2")+
    stat_summary(fun = mean, geom = "crossbar", 
                 width=0.06, colour = "purple")+
    labs(title = "Distribuzione dei guadagni per città negli anni",
         x= "Città", y="Guadagni (milioni di dollari)")+
    theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 13))+
    scale_y_continuous(limits = c(0, NA))+
    facet_wrap(~ year, ncol=2, scales = "free")

—

Bonus: un altro modo per ottenere una simile rappresentazione con 5 grafici creati in modo indipendente
library(ggpubr)  #utilizzo della funzione ggarrange per unire grafici
year_boxplot_vol <- function(data_subset, year) {
  ggplot(data_subset, aes(x = city, y = volume)) +
    geom_flat_violin(position = position_nudge(x = .2, y = 0), 
                     adjust = .5, trim = TRUE, fill = "pink2") +
    geom_jitter(width = 0.18, size = 1.15, aes(color = volume)) +
    scale_color_gradient(low = "orange4", high = "orange") +
    geom_boxplot(width = .15, outlier.color = "red", fill = "lightblue2") +
    stat_summary(fun = mean, geom = "crossbar", 
                 width = 0.06, colour = "purple") +
    labs(title = paste("Distribuzione dei guadagni per città -", year),
         x = "Città", y = "Guadagni (milioni di dollari)") +
    theme(plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
          axis.title = element_text(size = 10)) +
    scale_y_continuous(limits = c(0, NA))
}


#Per chiarezza, la funzione accetta un subset dei dati (in questo caso
#relativo a un anno) e restituisce il grafico
#Sarà richiamata per tutti gli anni fino ad ottenere i differenti 
#grafici tutti in un blocco
p1 <- year_boxplot_vol(subset(data, year == 2010), "2010")
p2 <- year_boxplot_vol(subset(data, year == 2011), "2011")
p3 <- year_boxplot_vol(subset(data, year == 2012), "2012")
p4 <- year_boxplot_vol(subset(data, year == 2013), "2013")
p5 <- year_boxplot_vol(subset(data, year == 2014), "2014")

#applicazione della funzione con ogni subset
ggarrange(p1, p2, p3, p4, p5, ncol = 2, nrow = 3)

#combinazione dei 5 grafici

In questo scenario però il risultato non è necessariamente migliore.

—

In ogni caso, si nota che Tyler abbia in genere i guadagni maggiori. Bryan-College Station ha dei guadagni elevati ma più “instabili”, con una fascia di distribuzione che si estende generalmente di più, con un massimo superiore a quello di Tyler in alcuni anni. Questo è coerente con la forte variabilità stagionale in termini di vendite che abbiamo visto in sezioni precedenti, superiore per B-C.S. rispetto alle altre città.

Wichita Falls ha invece una fascia molto ristretta in termini di guadagni, anche questo è coerente con il mercato “stazionario” e in piccola scala visto in precedenza, con prezzi modesti e numeri di vendite poco variabili tra stagioni e anni.

Non c’è molto da segnalare per Beaumont, se non nel 2012 dove la fascia è estremamente concentrata (alcuni valori vengono considerati addirittura outliers).

C’è da dire, comunque, che questo tipo di dataset presenta un numero di “osservazioni” relativamente limitato.

-

3: Confronto del totale delle vendite tra i vari mesi (e città) per mezzo di grafici a barre sovrapposte e normalizzati (bonus: aggiungere l’anno)

#Arrivato a questo punto mi accorgo che in precedenza ho inconsciamente utilizzato geom_bar impostando manualmente l'argomento stat="identity" invece che usare direttamente geom_col che lo usa di default.
#Pur funzionando correttamente, capisco che formalmente sia un po' improprio.
ggplot(data, aes(x=month, y=sales, fill=city)) +
  geom_col(position = "stack") +
  labs(title = "Totale delle vendite per città nei singoli mesi", x= "Mese", y="Vendite (unità)", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12, labels = month.name) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 17, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom")

-

Ricreiamo lo stesso grafico ma normalizzato (sufficiente cambiare l’argomento di geom_col):

ggplot(data, aes(x=month, y=sales, fill=city)) +
  geom_col(position = "fill") +
  labs(title = "Totale delle vendite per città nei singoli mesi", x= "Mese", y="Vendite (frazione relativa)", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12, labels = month.name) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 17, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom")

-

Adesso invece proviamo ad inserire anche l’anno come variabile di rappresentazione, usando nuovamente ggarrange.

year_barplot_mSal <- function(data_subset, year) {
  ggplot(data_subset, aes(x=month, y=sales, fill=city)) +
  geom_col(position = "stack") +
  labs(title = paste("Totale delle vendite per città nei singoli mesi -", 
                     year), 
       x= "Mese",
       y="Vendite (unità)", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12, labels = month.name) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        axis.title = element_text(size = 13),
        legend.position = "bottom")
  
}
bar_p1 <- year_barplot_mSal(subset(data, year == 2010), "2010")
bar_p2 <- year_barplot_mSal(subset(data, year == 2011), "2011")
bar_p3 <- year_barplot_mSal(subset(data, year == 2012), "2012")
bar_p4 <- year_barplot_mSal(subset(data, year == 2013), "2013")
bar_p5 <- year_barplot_mSal(subset(data, year == 2014), "2014")
ggarrange(bar_p1, bar_p2, bar_p3, bar_p4, bar_p5, ncol = 1, nrow = 5)

-

Proviamo infine a rifare lo stesso ma con normalizzazione. Questa volta proviamo più semplicemente con facet_grid:

ggplot(data, aes(x=month, y=sales, fill=city)) +
  geom_col(position = "fill") +
  labs(title = "Totale delle vendite per città nei singoli mesi", x= "Mese", y="Vendite (frazione relativa)", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12, labels = month.name) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 17, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom")+
  facet_grid(year ~ .)

Per curiosità, proviamo a fare lo stesso con gli anni come colonne:

ggplot(data, aes(x=month, y=sales, fill=city)) +
  geom_col(position = "fill") +
  labs(title = "Totale delle vendite per città nei singoli mesi", x= "Mese", y="Vendite (frazione relativa)", fill="Città") +
  scale_fill_brewer(palette = "Set1") +
  scale_x_continuous(breaks = 1:12) +
  scale_y_continuous(limits = c(0, NA)) + 
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 11),
        legend.position = "bottom")+
  facet_grid(~ year)

Non è grandiosamente elegante… ma ce lo faremo bastare.

-

4: Line chart di una variabile per confronti commentati fra città e periodi storici

Per questa operazione ci concentreremo sulla variabile sales. Oltre alla linechart aggiungeremo delle linee tratteggiate che definiscano vagamente la “tendenza”.
Osserviamo prima tutto assieme (anche se sarà molto confusionario):

ggplot(data, aes(x = (as.Date(paste(year, month, "01", sep = "-"))), 
                 y = sales, color = city))+
  geom_line()+
  geom_point(size=1.5)+
  geom_text(aes(label = sales), vjust = -0.7, size = 2.5)+
  geom_smooth(method = "loess", span = 0.3, se = F, 
              linetype = "dotdash", linewidth = 0.5)+
  labs(title = "Vendite nel tempo tra le città",
       x = "Anno",
       y = "Vendite (unità)",
       color = "Città") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits= c(0, NA))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom",
        panel.grid.major = element_line(linewidth=0.30, linetype = 'twodash',
                                        color = "grey60"),
        panel.grid.minor = element_line(linewidth = 0.4, linetype = 'dotted',
                                        color = "grey70"))+
  scale_color_brewer(palette = "Set1")

-

Focus sulle linee di tendenza, di cui un set che segue con maggiore precisione l’andamento dei singoli dati e un altro meno visibile che lo definisce più genericamente e linearmente:

ggplot(data, aes(x = (as.Date(paste(year, month, "01", sep = "-"))), 
                 y = sales, color = city))+
  geom_smooth(method = "loess", span = 0.15, se = F, 
              linetype = "twodash", linewidth = 1)+
  geom_smooth(method = "loess", span = 0.8, se = F, 
              linetype = "dotted", linewidth = 0.8)+
  labs(title = "Vendite nel tempo tra le città",
       x = "Anno",
       y = "Vendite (unità)",
       color = "Città") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits= c(0, NA))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom",
        panel.grid.major = element_line(linewidth=0.30, linetype = 'twodash',
                                        color = "grey60"),
        panel.grid.minor = element_line(linewidth = 0.4, linetype = 'dotted',
                                        color = "grey70"))+
  scale_color_brewer(palette = "Set1")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Si nota come Tyler si confermi la città con le maggiori vendite (superata solo per brevi periodi da Bryan-College Station).
Bryan-College Station ha un andamento particolarmente variabile, portandolo talvolta al di sotto di Beaumont.
Whichita Falls è la più stabile e l’unica che non sembra avere una crescita sostanziale negli anni.
Si nota anche che i picchi (sia massimi che minimi) non sono sempre perfettamente coincidenti tra città e anni.

-

Aggiungiamo la funzione geom_text_repel del pacchetto ggrepel per mitigare il problema della sovrapposizione delle etichette numeriche nella visualizzazione dell’intero periodo.

library(ggrepel)

Vediamo anche i grafici separatamente, limitandoci ad usare facet_wrap:

ggplot(data, aes(x = (as.Date(paste(year, month, "01", sep = "-"))), 
                 y = sales, color = city))+
  geom_line(linewidth=0.35)+
  geom_point(size=1.2)+
  geom_text_repel(aes(label = sales), size = 2.0, nudge_y = 2,
                  max.overlaps = 30, direction = "both",
                  max.iter = 50000, max.time = 1.5)+
  geom_smooth(method = "loess", span = 0.3, se = F, 
              linetype = "dotted", linewidth = 0.3)+
  labs(title = "Vendite nel tempo tra le città",
       x = "Anno",
       y = "Vendite (unità)",
       color = "Città") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits= c(0, NA))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title = element_text(size = 14),
        legend.position = "bottom",
        panel.grid.major = element_line(linewidth=0.30, linetype = 'twodash',
                                        color = "grey60"),
        panel.grid.minor = element_line(linewidth = 0.4, linetype = 'dotted',
                                        color = "grey70"))+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~ city)

-

Creiamo grafici differenti per i 5 anni, per visualizzare con maggiore dettaglio:

year_lineplot_mSal <- function(data_subset, year) {
  ggplot(data_subset, aes(x=month, y=sales, color=city)) +
  geom_line(linewidth = 0.5)+
  geom_point(size=1.3)+
  geom_text_repel(aes(label = sales), size = 2.5,
                  max.overlaps = 20, direction = "y",
                  max.iter = 30000, max.time = 1.2)+
  geom_smooth(method = "loess", span = 0.7, se = F, 
              linetype = "dotted", linewidth = 0.55)+
  labs(title = paste("Andamento delle vendite -", 
                     year),
       x = "Mese",
       y = "Vendite (unità)",
       color = "Città") +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  scale_y_continuous(limits = c(0, NA))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 10),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8),
        legend.key.size = unit(0.4, 'cm'),
        legend.key.height = unit(0.4, 'cm'),
        legend.key.width = unit(0.4, 'cm'),
        legend.position = "bottom",
        panel.grid.major = element_line(linewidth=0.30, linetype = 'twodash',
                                        color = "grey60"),
        panel.grid.minor = element_line(linewidth = 0.4, linetype = 'dotted',
                                        color = "grey70"))+
  scale_color_brewer(palette = "Set1")
    
}
line_p1 <- year_lineplot_mSal(subset(data, year == 2010), "2010")
line_p2 <- year_lineplot_mSal(subset(data, year == 2011), "2011")
line_p3 <- year_lineplot_mSal(subset(data, year == 2012), "2012")
line_p4 <- year_lineplot_mSal(subset(data, year == 2013), "2013")
line_p5 <- year_lineplot_mSal(subset(data, year == 2014), "2014")

-

-

Visualizziamo, a due a due, confronti tra un anno e il successivo con ggarrange:

-

library(RColorBrewer)

-

Concludiamo visualizzando grafici per città invece che per anno, facendo anche in modo di indicare i valori massimi per ogni anno (nota: idea iniziale suggerita dall’AI, poi opportunamente personalizzata):

city_lineplot_mSal <- function(data_subset, city) {
  max_points <- data_subset %>%
        group_by(year) %>%
        summarize(max_sales = max(sales), .groups = 'drop') %>%
        left_join(data_subset, by = c("year", "max_sales" = "sales"))
  
  ggplot(data_subset, aes(x=month, y=sales, color=factor(year))) +
  geom_line(linewidth = 0.45)+
  geom_point(size=1.3)+
  geom_point(data = max_points, aes(x = month, y = max_sales), 
             shape = 17, size = 3) +
  geom_text_repel(data = max_points, 
                  aes(x = month, y = max_sales, 
                      label = max_sales),
                  max.overlaps = 20, direction = "both",
                  max.iter = 30000, max.time = 1.2,
                  size = 3.8, fontface= "bold")+
  geom_smooth(method = "loess", span = 0.7, se = F, 
              linetype = "dotted", linewidth = 0.3)+
  labs(title = paste("Andamento delle vendite -", 
                     city),
       x = "Mese",
       y = "Vendite (unità)",
       color = "Anno") +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  scale_y_continuous(limits = c(0, NA))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title = element_text(size = 10),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8),
        legend.key.size = unit(0.4, 'cm'),
        legend.key.height = unit(0.4, 'cm'),
        legend.key.width = unit(0.4, 'cm'),
        legend.position = "bottom",
        panel.grid.major = element_line(linewidth=0.30, linetype = 'twodash',
                                        color = "grey60"),
        panel.grid.minor = element_line(linewidth = 0.4, linetype = 'dotted',
                                        color = "grey70"))+
  scale_color_manual(values = brewer.pal(7, "YlOrRd")[3:7])
    
}
c_line_p1 <- city_lineplot_mSal(subset(data, city == "Beaumont"), "Beaumont")
c_line_p2 <- city_lineplot_mSal(subset(data, city == "Bryan-College Station"), "Bryan-College Station")
c_line_p3 <- city_lineplot_mSal(subset(data, city == "Tyler"), "Tyler")
c_line_p4 <- city_lineplot_mSal(subset(data, city == "Wichita Falls"), "Wichita Falls")

E con questo siamo giunti al termine di quest’analisi descrittiva.

Tra i punti salienti abbiamo:

-Mercato in leggera regressione nel 2011, ma in discreta crescita negli anni successivi.

-Tyler è la città col maggior mercato, maggiori vendite e guadagni ma anche tantissimi annunci attivi, Whichita Falls quella col mercato più stabile e in piccola scala.

-Bryan-College Station è la città col maggior rapporto percentuale di vendite/annunci, i prezzi generalmente più alti (di poco più di Tyler) e un mercato assai variabile.

-I prezzi mediani sono sensibilmente aumentati negli anni, generalmente.

-Il numero di annunci attivi ha avuto picco nel 2011 ed è sceso rapidamente negli anni successivi,
mentre le vendite sono aumentate.

-R Markdown è un po’ buggato (ok questo non è rilevante).

The end.