Business Data Analytics

ANALISI ESPLORATIVA DATASET

library(knitr)

Grafico di una variabile

library(ggplot2)

ggplot(db, aes(y=GDPpc)) + geom_boxplot()

Boxplot PIL procapite per ciascuna regione (in ordine di PIL procapite)

ggplot(db, aes(x=reorder(Nazione, -GDPpc, FUN = median), y=GDPpc)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 90, vjust= 0.5, hjust=1))

Boxplot PIL procapite rispetto a dati trattati e non (per vedere se sono omogenei)

ggplot(db, aes(x=Classe, y=GDPpc, fill = Classe)) + geom_boxplot()

Matrice di correlazione tra le variabili

library(corrplot)

col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF","#77AADD", "#4477AA"))

cor_matr = cor(db[,4:9])

quartz()

corrplot(cor_matr, method="color", col=col(200), type="upper",order="hclust", addCoef.col = "black", tl.col="black", tl.srt=45, diag=FALSE)

02 CAUSALITY PSM

PSM:

Propensity Score Matching

La Propensity Score Matching è una tecnica utilizzata quando la regressione lineare semplice di tipo LS ci restituisce delle stime biased nei confronti dei coefficienti beta delle covariate di cui vogliamo studiare l'effetto causale rispetto alla variabile dipendente (di outcome). In particolare, viene utilizzata nel momento in cui abbiamo a che fare con due campioni, uno trattato e uno no, che hanno delle caratteristiche in termini di covariate abbastanza differenti. Quindi il PSM prende il campione dei trattati e cerca all'interno del campione dei non trattati gli individui più simili, in termini di covariate, agli individui più simili nei campioni trattati.

1) Setup

library(stargazer)
library(MatchIt)
library(MASS)
library(mvtnorm)
library(ivreg)
library(rdrobust)
db = read.csv('DB_Causality.csv', sep = ";")

Boxplot del GDP, trattati vs non trattati

library(ggplot2)
ggplot(db, aes(x=Classe, y=GDPpc, fill = Classe)) + geom_boxplot()

2) Analisi

pre-matchingdb_trattati = db[which(db$Trattamento_Fondi_Obiettivo1 == 1), ]
db_non_trattati = db[which(db$Trattamento_Fondi_Obiettivo1 == 0), ]
t.test(db_trattati$GDPpc, db_non_trattati$GDPpc)
t.test(db_trattati$Accumulo_Capitale, db_non_trattati$Accumulo_Capitale)
t.test(db_trattati$Scolarizzazione, db_non_trattati$Scolarizzazione)
t.test(db_trattati$Crescita_Popolazione,db_non_trattati$Crescita_Popolazione)

Se p-value è alto, allora accetto H0, cioè le medie dei due array sono statisticamente molto simili.
Un altro modo per accettare H0 è vedere il valore assoluto della statistica test t (se è alto allora rifiuto H0, cioè le due medie non sono simili).

3) Calcolo Naive Estimator NE
NE è una differenza tra medie, la media dei trattati meno i non trattati.
library(sjPlot)
naive_estimator = mean(db_trattati$Crescita_GDPpc)-mean(db_non_trattati$Crescita_GDPpc)

4) Stimo modello di regressione
Normalizziamo i dati. Se va fatto, viene specificato

all’esame.

Crescita_GDPpc = scale(db$Crescita_GDPpc)

Trattamento_Fondi_Obiettivo1 = db$Trattamento_Fondi_Obiettivo1

GDPpc = scale(db$GDPpc)

Accumulo_Capitale = scale(db$Accumulo_Capitale)

Scolarizzazione = scale(db$Scolarizzazione)

Crescita_Popolazione = scale(db$Crescita_Popolazione)

f_ne = Crescita_GDPpc ~ Trattamento_Fondi_Obiettivo1 + GDPpc +Accumulo_Capitale + Scolarizzazione + Crescita_Popolazione

dataset = data.frame(Crescita_GDPpc,Trattamento_Fondi_Obiettivo1,GDPpc, Accumulo_Capitale, Scolarizzazione, Crescita_Popolazione)

mod_ne = lm(f_ne, data = dataset)

summary(mod_ne)

Call:lm(formula = f_ne, data = dataset)

Residuals:Min 1Q Median 3Q Max-2.6827 -0.3559 0.0291 0.4128 3.2271

Coefficients: Estimate Std. Error t value Pr(>|t|)

(Intercept) -0.250187 0.065738 -3.806 0.000177 ***

Trattamento_Fondi_Obiettivo1 0.942735 0.160096 5.889 1.23e-08 ***

GDPpc -0.085074 0.075871 -1.121 0.263221

Accumulo_Capitale -0.005064 0.050919 -0.099 0.920858

Scolarizzazione -0.247953 0.053593 -4.627 5.93e-06

Crescita_Popolazione -0.072538 0.062201 -1.166 0.244635---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8089 on 254 degrees of freedom

Multiple R-squared: 0.3584, Adjusted R-squared: 0.3457

F-statistic: 28.37 on 5 and 254 DF, p-value: < 2.2e-16

La variabile Trattamento_Fondi_Obiettivo1 ha un beta positivo e statisticamente significativo (perché il valore assoluto della statistica test t-value è alto e perché il p-value è molto basso). Ne deduciamo che le regioni trattate (quelle per cui Trattamento_Fondi_Obiettivo1 = 1) hanno un y (Crescita_GDPpc) maggiore. Però, Trattamento_Fondi_Obiettivo1 potrebbe essere biased, a causa dell'assenza di randomizzazione.

5) Propensity Score Matching (PSM)

Il PSM stima il propensity score, sceglie un algoritmo di matching, il quale, dato il PSM, associa ad ogni unità trattata un'unità non trattata.

Dato il campione dei dati trattati e non trattati, matchati rispetto al propensity score, calcoliamo il nuovo modello di regressione. Il propensity score viene stimato mediante regressione logistica. Esistono diversi tipi di algoritmi di matching, ma quello su cui ci focalizziamo è il K-Nearest Neighbour (k-NN). Dato ogni elemento nel gruppo dei trattati, andiamo a cercare i k individui non trattati con il propensity score più vicino a quello dell'unità trattata. Si può usare in due modi:

• Con replacement, ad ogni unità non trattata può essere matchata una o più unità trattate in contemporanea.
• Senza replacement, ad ogni unità non trattata viene matchata una singola unità trattata.

library(MatchIt)
psm_result = matchit(db$Trattamento_Fondi_Obiettivo1 ~ db$GDPpc + db$Accumulo_Capitale + db$Scolarizzazione + db$Crescita_Popolazione, method = "nearest")
method = "nearest" per specificare che vogliamo k-NN.

Di default fa il k-NN senza replacement.

6) Analisi post-matching

Vediamo se abbiamo risolto il problema dell'assenza di randomizzazione del campione.

db_post_psm = match.data(psm_result, data=db)
db_post_psm_trattati = db_post_psm[which(db_post_psm$Trattamento_Fondi_Obiettivo1 == 1), ]
db_post_psm_non_trattati = db_post_psm[which(db_post_psm$Trattamento_Fondi_Obiettivo1 == 0), ]
(db_post_psm_trattati$GDPpc, db_post_psm_non_trattati$GDPpc)
t.test(db_post_psm_trattati$Accumulo_Capitale, db_post_psm_non_trattati$Accumulo_Capitale)
t.test(db_post_psm_trattati$Scolarizzazione, db_post_psm_non_trattati$Scolarizzazione)
t.test(db_post_psm_trattati$Crescita_Popolazione, db_post_psm_non_trattati$Crescita_Popolazione)

Se p-value è alto, accetto H0, cioè le medie tra trattati e non trattati rispetto a quella covariata sono più o meno uguali, quindi ho risolto il problema per quella covariata, altrimenti no.

Graficamente

ggplot(db_post_psm, aes(x=Classe, y=GDPpc, fill = Classe))

+geom_boxplot()
NE post-matching

naive_estimator_post_psm = mean(db_post_psm_trattati$Crescita_GDPpc) - mean(db_post_psm_non_trattati$Crescita_GDPpc)
Stimo il nuovo modello di regressione

Crescita_GDPpc_post_psm = scale(db_post_psm$Crescita_GDPpc)
Trattamento_Fondi_Obiettivo1_post_psm = db_post_psm$Trattamento_Fondi_Obiettivo1
GDPpc_post_psm = scale(db_post_psm$GDPpc)
Accumulo_Capitale_post_psm = scale(db_post_psm$Accumulo_Capitale)
Scolarizzazione_post_psm = scale(db_post_psm$Scolarizzazione)
Crescita_Popolazione_post_psm = scale(db_post_psm$Crescita_Popolazione)
f_ne_post_psm = Crescita_GDPpc_post_psm ~ Trattamento_Fondi_Obiettivo1_post_psm + GDPpc_post_psm + Accumulo_Capitale_post_psm + Scolarizzazione_post_psm + Crescita_Popolazione_post_psm
dataset_post_psm = data.frame(Crescita_GDPpc_post_psm, Trattamento_Fondi_Obiettivo1_post_psm, GDPpc_post_psm, Accumulo_Capitale_post_psm, Scolarizzazione_post_psm, Crescita_Popolazione_post_psm)
mod_ne_post_psm = lm(f_ne_post_psm, data =