Laboratorio R per le scienze sociali – Giorno 9
Federico Vegetti
Regressione con R
In questa ultima lezione vedremo come fare una analisi di regressione con R. La regressione è uno strumento molto comune per analizzare l’associazione tra uno o più variabili indipendenti e una variabile dipendente. L’obiettivo di questa sessione non è quello di insegnarvi la teoria statistica dietro la regressione. Tuttavia abbiamo bisogno di capire abbastanza teoria per essere sicuri che quello che stiamo facendo sia corretto, e per essere in grado di interpretare i risultati della nostra analisi. Se lo facciamo, la statistica può davvero aiutarci a estrarre conoscenza dai dati.
Useremo sia i dati V-Dem che una nuova versione dei dati CSES con più variabili.
Prima di tutto, come al solito, carichiamo le librerie e i dati.
library(rio)
library(tidyverse)
setwd("/folder/where/you/keep/the/data")
vdem <- import("vdem_small.dta")
cses <- import("cses2018new.xlsx")Come promemoria, il dataset V-Dem contiene le seguenti variabili:
country: il nome del paese (tutti i paesi dell’UE-28)year: l’anno di osservazione (dal 1990 al 2018)liberal: l’indice di democrazia liberale pubblicato da V-Dem, che cattura quanto è importante in un dato paese e anno “proteggere i diritti individuali e delle minoranze contro la tirannia dello Stato e la tirannia della maggioranza”. Varia da 0 (bassa protezione) a 1 (alta protezione)free_exp: l’indice di libertà di espressione pubblicato da V-Dem, che cattura quanto il governo di un paese in un determinato anno rispetta “la libertà di stampa e dei media, la libertà della gente comune di discutere questioni politiche a casa e nella sfera pubblica, così come la libertà di espressione accademica e culturale”. Varia da 0 (bassa libertà) a 1 (alta libertà)gdp_pc_log: il PIL pro capite (logaritmo)postcom: una variabile dummy che ha valore 1 per i paesi dell’ex blocco comunista, e 0 per gli altri
La nuova versione del dataset CSES contiene le seguenti variabili:
id: codice intervistatoage: età dell’intervistatoage_gr: età dell’intervistato, raggruppata in 3 categoriefemale: variabile dummy che indica il sesso dell’intervistato. 1 = femmine, 0 = maschiedu_low: variabile dummy che indica che l’intervistato ha un basso titolo di studio (fino alla scuola media)edu_high: variabile dummy che indica che l’intervistato ha un basso titolo di studio (laurea, dottorato)mip1, mip2: le due risposte alla domanda aperta “Secondo lei, quali sono i due problemi più importanti che l’Italia deve affrontare in questo momento?”eco_eval: la valutazione di come lo stato dell’economia del paese è cambiato nel corso dell’anno precedente, codificato in 3 classi (1 = è migliorato; 0 = è rimasto uguale; -1 = è peggiorato)area: l’area geografica di residenzavoted: variabile dummy che indica se l’intervistato ha votato alle elezioni 2018 (1) o no (0)ptv_PD: una variabile ordinale che misura la probabilità espressa che l’intervistato voterà mai per il PD. Va da 0 (“Non voterò mai per il partito”) a 10 (“Sicuramente un giorno voterò per il partito”)ptv_FI: una variabile ordinale che misura la probabilità espressa che l’intervistato voterà mai per Forza Italia. Va da 0 (“Non voterò mai per il partito”) a 10 (“Sicuramente un giorno voterò per il partito”)ptv_M5S: una variabile ordinale che misura la probabilità espressa che l’intervistato voterà mai per il M5S. Va da 0 (“Non voterò mai per il partito”) a 10 (“Sicuramente un giorno voterò per il partito”)ptv_LEGA: una variabile ordinale che misura la probabilità espressa che l’intervistato voterà mai per la Lega. Va da 0 (“Non voterò mai per il partito”) a 10 (“Sicuramente un giorno voterò per il partito”)
Stimare una regressione lineare
Iniziamo con i dati V-Dem. Vogliamo capire quali sono le variabili (tra quelle che abbiamo nel nostro dataset) che influenzano il grado di democrazia liberale in un paese, misurato con l’indice liberal: questa sarà la nostra variabile dipendente. (Nota: quando parlo di “influenzare” o “effetto” intendo solo in senso statistico. Facendo una semplice regressione non possiamo dire nulla sulla direzione della causalità tra le variabili dipendenti e quelle indipendenti).
Procediamo per gradi. Per iniziare, vediamo come liberal cambia nel tempo per tutti i paesi, guardando l’effetto della variabile year. Prima di tutto, possiamo vedere se c’è qualche associazione lineare facendo un grafico delle due variabili, aggiungendo il layer geom_smooth() di ggplot2. Specificando method = "lm", chiediamo a ggplot2 di darci una linea retta. In questo modo ggplot2 stimerà un modello di regressione lineare con i dati, e traccerà la linea di regressione risultante insieme a un intervallo di confidenza.
Quindi ci aspettiamo che year abbia un effetto positivo e significativo su liberal. Nel 1990, il livello medio di liberal nel nostro campione era un po’ più basso di \(0,88\), mentre nel 2018 era un po’ più basso di \(0,91\). Un modello di regressione dovrebbe mostrare esattamente lo stesso risultato, anche se in modo diverso.
In R, possiamo stimare una regressione lineare usando la funzione lm() (che sta per “linear model”), che fa già parte di R base tramite la libreria stats. La funzione prende diversi argomenti, ma i due più importanti sono:
- La formula che specifica l’equazione da calcolare.
- I dati a cui il modello deve essere applicato.
La formula è l’elemento più importante: è qui che specifichiamo la forma funzionale della relazione tra la nostra variabile dipendente, \(Y\), e la variabile indipendente (o le variabili indipendenti), \(X\). Ricordate: la regressione lineare applica una funzione lineare ai dati, e trova i parametri della funzione che minimizzano la somma dell’errore quadratico. Una funzione lineare ha la seguente forma:
\[ Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_n X_n \]
Dove:
- \(\alpha\) è la costante, o intercetta, il valore medio di \(Y\) quando tutti gli \(X\) hanno valore \(0\).
- \(\beta_n\) è il coefficiente della variabile \(X_n\), e descrive quanto \(Y\) cambia in media per ogni aumento unitario di \(X_n\)
In R, le formule sono scritte nella forma Y ~ X1 + X2 + Xn, usando la tilde ~ per specificare che la variabile dipendente \(Y\) è una funzione delle variabili \(X\), e il segno più + per elencare tutte le variabili indipendenti. Nel nostro caso, \(Y\) è liberal e l’unica \(X\) che abbiamo è year, quindi la formula sarebbe liberal ~ year.
Stimiamo il modello di regressione e mettiamo i risultati in un oggetto chiamato m.v:
Ora possiamo vedere i risultati applicando la funzione summary() all’oggetto m.v:
##
## Call:
## lm(formula = liberal ~ year, data = vdem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50058 -0.01542 0.01984 0.04916 0.10342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.4982644 0.7158325 -2.093 0.036658 *
## year 0.0011929 0.0003572 3.340 0.000877 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08447 on 805 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.01367, Adjusted R-squared: 0.01244
## F-statistic: 11.15 on 1 and 805 DF, p-value: 0.0008773Due cose:
- Non è necessario mettere il risultato del modello in un oggetto, potremmo mettere la funzione
lm()direttamente all’interno della funzionesummary(), tuttavia alcuni modelli sono molto lenti da stimare, per cui è sempre meglio salvare i risultati in un oggetto che possiamo richiamare in futuro. - Tecnicamente, la formula dovrebbe includere un
1all’inizio, per specificare che vogliamo stimare anche il termine costante \(\alpha\), tuttavia questo viene fatto in automatico. Al contrario, se vogliamo adattare un modello senza intercetta, dovremo specificarlo nella formulaY ~ -1 + X1 + X2 + Xn.
Tornando ai risultati, la costante ha un valore che è al di fuori del range della variabile liberal, che va da 0 a 1. Questo succede perché l’intercetta \(\alpha\) rappresenta la media di \(Y\) quando tutti gli \(X\) sono 0, e nel nostro caso, la variabile year parte dal 1990. Quindi il valore di -1.498 che abbiamo ottenuto è il livello medio di liberal nell’anno 0. Ecco perché è utile centrare le variabili indipendenti (come abbiamo visto qualche settimana fa): per ottenere un coefficiente dell’intercetta che abbia un senso.
In questo caso creiamo un’altra variabile chiamata year_r che va da 0 (1990) a 28 (2018), così il valore dell’intercetta rifletterà la media di liberal nell’anno 1990:
# Creiamo year_r
vdem$year_r <- vdem$year - 1990
# Stimiamo la regressione con "year_r" al posto di "year"
m.v <- lm(liberal ~ year_r,
data = vdem)
summary(m.v)##
## Call:
## lm(formula = liberal ~ year_r, data = vdem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50058 -0.01542 0.01984 0.04916 0.10342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8755834 0.0058434 149.84 < 0.0000000000000002 ***
## year_r 0.0011929 0.0003572 3.34 0.000877 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08447 on 805 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.01367, Adjusted R-squared: 0.01244
## F-statistic: 11.15 on 1 and 805 DF, p-value: 0.0008773Notare due cose:
- Ora l’intercetta è
0,8755, che corrisponde a quanto osservato nel grafico sopra - Il coefficiente di
year_rrimane lo stesso di quello diyearnel modello precedente
Regressione multivariata
Il vero punto di forza della regressione è che possiamo andare oltre le relazioni bivariate (che potremmo valutare con un grafico, o con un coefficiente di correlazione). Per farlo dobbiamo solo aggiornare la formula con tutte le variabili indipendenti che vogliamo includere nel modello. Aggiungiamo anche la variabile dummy postcom, per tenere conto del fatto che nel nostro campione ci sono paesi con un passato recente di regime autoritario, e la variabile gdp_pc_log, per tenere conto delle differenze nel PIL pro capite tra paesi.
##
## Call:
## lm(formula = liberal ~ gdp_pc_log + postcom + year_r, data = vdem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42636 -0.02254 0.01583 0.04274 0.14165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1749160 0.1048117 -1.669 0.09556 .
## gdp_pc_log 0.1072681 0.0105022 10.214 < 0.0000000000000002 ***
## postcom 0.0194981 0.0095351 2.045 0.04122 *
## year_r -0.0014385 0.0004656 -3.090 0.00208 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07426 on 747 degrees of freedom
## (57 observations deleted due to missingness)
## Multiple R-squared: 0.2444, Adjusted R-squared: 0.2414
## F-statistic: 80.55 on 3 and 747 DF, p-value: < 0.00000000000000022Variabili indipendenti categoriali
R può includere variabili indipendenti categoriali nei modelli di regressione con facilità. Come con Stata o SPSS, ogni volta che si includono variabili categoriche in un modello di regressione, dovrebbero essere trattate come variabili factor. In questo modo, il modello tratterà ogni livello della variabile come se fosse una variabile dummy (con valore 1 se quel livello è presente e 0 se non è presente), e li includerà tutti * tranne uno*. Quest’ultimo sarà usato come categoria di riferimento.
Questo è più facile a farsi che a dirsi: vogliamo aggiungere la variabile country al nostro modello di regressione, perché vogliamo controllare le caratteristiche specifiche dei diversi paesi che non sono catturate dalle altre variabili (e quindi rimuoviamo postcom perché è solo una versione più “grossolana” di una variabile che identifica il paese). La variabile country è una variabile stringa, il che significa che non ha alcun valore numerico. Quindi non abbiamo alternative a trattarla come un insieme di variabili dummy, una per paese, che saranno incluse una per una nel modello. La buona notizia è che questo viene fatto da R in automatico quando specifichiamo che una variabile indipendente nel modello di regressione è una variabile factor, utilizzando la funzione factor():
##
## Call:
## lm(formula = liberal ~ gdp_pc_log + year_r + factor(country),
## data = vdem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.300166 -0.009607 0.002297 0.017050 0.156007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.35488984 0.13009961 2.728 0.006530 **
## gdp_pc_log 0.05379419 0.01286640 4.181 0.00003259552773324 ***
## year_r 0.00017073 0.00042709 0.400 0.689463
## factor(country)Belgium 0.02265242 0.01181470 1.917 0.055594 .
## factor(country)Bulgaria 0.01769709 0.01867206 0.948 0.343557
## factor(country)Croatia -0.18226488 0.01599683 -11.394 < 0.0000000000000002 ***
## factor(country)Cyprus -0.05397978 0.01272520 -4.242 0.00002504577382082 ***
## factor(country)Czech Republic 0.01209336 0.01406125 0.860 0.390048
## factor(country)Denmark 0.05799294 0.01177612 4.925 0.00000104878035823 ***
## factor(country)Estonia 0.04361005 0.01643941 2.653 0.008159 **
## factor(country)Finland 0.05163341 0.01184831 4.358 0.00001504161428939 ***
## factor(country)France -0.01444402 0.01186475 -1.217 0.223854
## factor(country)Germany 0.05226527 0.01177571 4.438 0.00001048151100835 ***
## factor(country)Greece -0.00883044 0.01299922 -0.679 0.497162
## factor(country)Hungary -0.01746448 0.01536637 -1.137 0.256109
## factor(country)Ireland -0.02606138 0.01180527 -2.208 0.027586 *
## factor(country)Italy -0.00559723 0.01185006 -0.472 0.636829
## factor(country)Latvia -0.00002962 0.01714797 -0.002 0.998622
## factor(country)Lithuania 0.06691258 0.01690128 3.959 0.00008271495256470 ***
## factor(country)Luxembourg -0.05706389 0.01288943 -4.427 0.00001102582121599 ***
## factor(country)Malta -0.14666791 0.01385883 -10.583 < 0.0000000000000002 ***
## factor(country)Netherlands 0.04268775 0.01181529 3.613 0.000324 ***
## factor(country)Poland 0.04154527 0.01647503 2.522 0.011893 *
## factor(country)Portugal 0.04704948 0.01348288 3.490 0.000513 ***
## factor(country)Romania -0.15659043 0.01969436 -7.951 0.00000000000000715 ***
## factor(country)Slovakia -0.03571675 0.01545041 -2.312 0.021075 *
## factor(country)Slovenia 0.01221482 0.01320984 0.925 0.355444
## factor(country)Spain 0.01735027 0.01237202 1.402 0.161232
## factor(country)Sweden 0.04546530 0.01177546 3.861 0.000123 ***
## factor(country)United Kingdom 0.01789274 0.01181578 1.514 0.130386
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04327 on 721 degrees of freedom
## (57 observations deleted due to missingness)
## Multiple R-squared: 0.7524, Adjusted R-squared: 0.7425
## F-statistic: 75.56 on 29 and 721 DF, p-value: < 0.00000000000000022Notare che in questo caso la categoria di base (o categoria “di riferimento”) è l’Austria. Questo significa che il coefficiente che osserviamo per ogni paese è la differenza nel livello medio di liberal tra quel paese e l’Austria. R usa come categoria di riferimento la prima categoria nell’ordine dei valori, che nel caso delle variabili stringa è la prima lettera dell’ordine alfabetico. Questo comportamento può essere cambiato usando la funzione relevel(). Per esempio, usiamo l’Italia come categoria di base:
m.v <- lm(liberal ~ gdp_pc_log + year_r + relevel(factor(country), ref = "Italy"),
data = vdem)
summary(m.v)##
## Call:
## lm(formula = liberal ~ gdp_pc_log + year_r + relevel(factor(country),
## ref = "Italy"), data = vdem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.300166 -0.009607 0.002297 0.017050 0.156007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3492926 0.1287750 2.712 0.006839 **
## gdp_pc_log 0.0537942 0.0128664 4.181 0.00003259552773324 ***
## year_r 0.0001707 0.0004271 0.400 0.689463
## relevel(factor(country), ref = "Italy")Austria 0.0055972 0.0118501 0.472 0.636829
## relevel(factor(country), ref = "Italy")Belgium 0.0282497 0.0117811 2.398 0.016744 *
## relevel(factor(country), ref = "Italy")Bulgaria 0.0232943 0.0176616 1.319 0.187612
## relevel(factor(country), ref = "Italy")Croatia -0.1766676 0.0151407 -11.668 < 0.0000000000000002 ***
## relevel(factor(country), ref = "Italy")Cyprus -0.0483826 0.0122835 -3.939 0.00008983574127156 ***
## relevel(factor(country), ref = "Italy")Czech Republic 0.0176906 0.0133819 1.322 0.186595
## relevel(factor(country), ref = "Italy")Denmark 0.0635902 0.0118648 5.360 0.00000011238326119 ***
## relevel(factor(country), ref = "Italy")Estonia 0.0492073 0.0155521 3.164 0.001621 **
## relevel(factor(country), ref = "Italy")Finland 0.0572306 0.0117755 4.860 0.00000144006774509 ***
## relevel(factor(country), ref = "Italy")France -0.0088468 0.0117761 -0.751 0.452748
## relevel(factor(country), ref = "Italy")Germany 0.0578625 0.0118591 4.879 0.00000131201724705 ***
## relevel(factor(country), ref = "Italy")Greece -0.0032332 0.0124949 -0.259 0.795892
## relevel(factor(country), ref = "Italy")Hungary -0.0118672 0.0145490 -0.816 0.414956
## relevel(factor(country), ref = "Italy")Ireland -0.0204642 0.0119730 -1.709 0.087847 .
## relevel(factor(country), ref = "Italy")Latvia 0.0055676 0.0162085 0.343 0.731323
## relevel(factor(country), ref = "Italy")Lithuania 0.0725098 0.0159757 4.539 0.00000662883905038 ***
## relevel(factor(country), ref = "Italy")Luxembourg -0.0514667 0.0134840 -3.817 0.000147 ***
## relevel(factor(country), ref = "Italy")Malta -0.1410707 0.0132070 -10.682 < 0.0000000000000002 ***
## relevel(factor(country), ref = "Italy")Netherlands 0.0482850 0.0119974 4.025 0.00006310168760566 ***
## relevel(factor(country), ref = "Italy")Poland 0.0471425 0.0155754 3.027 0.002560 **
## relevel(factor(country), ref = "Italy")Portugal 0.0526467 0.0128885 4.085 0.00004906397073946 ***
## relevel(factor(country), ref = "Italy")Romania -0.1509932 0.0186471 -8.097 0.00000000000000239 ***
## relevel(factor(country), ref = "Italy")Slovakia -0.0301195 0.0146665 -2.054 0.040372 *
## relevel(factor(country), ref = "Italy")Slovenia 0.0178121 0.0126636 1.407 0.159990
## relevel(factor(country), ref = "Italy")Spain 0.0229475 0.0120313 1.907 0.056876 .
## relevel(factor(country), ref = "Italy")Sweden 0.0510625 0.0118515 4.309 0.00001871475705537 ***
## relevel(factor(country), ref = "Italy")United Kingdom 0.0234900 0.0117807 1.994 0.046536 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04327 on 721 degrees of freedom
## (57 observations deleted due to missingness)
## Multiple R-squared: 0.7524, Adjusted R-squared: 0.7425
## F-statistic: 75.56 on 29 and 721 DF, p-value: < 0.00000000000000022Esportare i risultati della regressione: texreg
Usare la funzione summary() sull’oggetto prodotto dalla regressione ci permette di vedere i risultati del nostro modello. Tuttavia, come esportiamo i risultati? Abbiamo davvero bisogno di copiare/incollare ogni singolo numero in una tabella di Excel o Word? E se vogliamo esportare i nostri risultati in PDF o HTML usando R Markdown? Abbiamo bisogno di mostrare il risultato come output di R come abbiamo fatto finora?
Naturalmente no. R ha un sacco di librerie per esportare i risultati di una regressione in diversi formati. È una questione di gusti quale usare, qui vedremo quello che trovo più flessibile, la libreria texreg. texreg vi permette di produrre codice HTML o LaTeX per produrre una tabella ben formattata basata sul risultato della regressione (ad esempio come l’oggetto m.v che abbiamo usato finora).
Cominciamo stimando un modello diverso (e migliore) usando i dati CSES. Vogliamo vedere quali sono i fattori che influenzano la preferenza degli elettori per il Movimento 5 Stelle (M5S). Misuriamo questa preferenza usando la variabile “probabilità di voto” (PTV) per il M5S, che nel nostro dataset si chiama ptv_M5S. Come promemoria, questa è una variabile che va da 0 a 10, ed è la risposta alla seguente domanda: “Quanto è probabile che lei voterà mai per il M5S?”. L’opzione di risposta più bassa è \(0\), e significa “Non voterei mai per il M5S”. La più alta è 10, e significa “Voterò sicuramente per il M5S prima o poi”.
Vogliamo vedere se l’età, il sesso, il livello di istruzione (basso o alto), la valutazione dell’economia e la zona geografica degli intervistati sono associati a una PTV più alta o più bassa per il M5S. Poiché il PTV va da 0 a 10, usiamo di nuovo una regressione lineare.
m.c <- lm(ptv_M5S ~ age + female + edu_low + edu_high +
eco_eval + factor(area),
data = cses)
summary(m.c)##
## Call:
## lm(formula = ptv_M5S ~ age + female + edu_low + edu_high + eco_eval +
## factor(area), data = cses)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5989 -3.0861 0.3361 2.6017 7.2649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.579138 0.319963 17.437 < 0.0000000000000002 ***
## age -0.021860 0.004871 -4.488 0.0000076424042 ***
## female -0.437209 0.159397 -2.743 0.006151 **
## edu_low 0.008780 0.230631 0.038 0.969635
## edu_high -0.639614 0.179547 -3.562 0.000377 ***
## eco_eval -0.723595 0.105931 -6.831 0.0000000000115 ***
## factor(area)Nord-Est 0.140031 0.254405 0.550 0.582095
## factor(area)Nord-Ovest -0.501856 0.225743 -2.223 0.026331 *
## factor(area)Sud 0.995694 0.218328 4.561 0.0000054471305 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.364 on 1797 degrees of freedom
## (195 observations deleted due to missingness)
## Multiple R-squared: 0.08234, Adjusted R-squared: 0.07825
## F-statistic: 20.16 on 8 and 1797 DF, p-value: < 0.00000000000000022Come possiamo esportare questa tabella in modo che possa essere resa in un bel formato in questo documento? Cominciamo installando la libreria texreg.
texreg ha diverse funzioni che sono pensate per esportare l’output di una varietà di modelli di regressione in un modo che appare più bello da vedere (e leggibile) di quello che si ottiene con la funzione summary(). Per esempio, un equivalente funzionale di summary(), che mostra l’output sulla console ma non è destinato all’esportazione, è la funzione screenreg():
##
## ===================================
## Model 1
## -----------------------------------
## (Intercept) 5.58 ***
## (0.32)
## age -0.02 ***
## (0.00)
## female -0.44 **
## (0.16)
## edu_low 0.01
## (0.23)
## edu_high -0.64 ***
## (0.18)
## eco_eval -0.72 ***
## (0.11)
## factor(area)Nord-Est 0.14
## (0.25)
## factor(area)Nord-Ovest -0.50 *
## (0.23)
## factor(area)Sud 1.00 ***
## (0.22)
## -----------------------------------
## R^2 0.08
## Adj. R^2 0.08
## Num. obs. 1806
## ===================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Tuttavia, siamo interessati a due funzioni in particolare: htmlreg(), per esportare tabelle in HTML, e texreg(), per esportare tabelle in LaTeX (e quindi in PDF).
Nota: In texreg le funzioni prendono come input un oggetto prodotto da una regressione, e stampano nella consolle il codice per esportare la tabella formattata, sia in HTML che in LaTeX. Qui vediamo come funziona quando è integrato con R Markdown. Tuttavia, potete anche usare texreg in una normale sessione di R, e copiare/incollare il codice prodotto dalle funzioni in un file HTML o LaTeX.
htmlreg()
htmlreg() produrrà del codice HTML per fare una tabella. Può essere usato in R Markdown quando creiamo un documento HTML o Word. Inoltre, permette di salvare una tabella in un file formato .doc, che può essere aperto in Word. Come con summary() e screenreg(), l’input è il risultato di un modello di regressione:
| Model 1 | |
|---|---|
| (Intercept) | 5.58*** |
| (0.32) | |
| age | -0.02*** |
| (0.00) | |
| female | -0.44** |
| (0.16) | |
| edu_low | 0.01 |
| (0.23) | |
| edu_high | -0.64*** |
| (0.18) | |
| eco_eval | -0.72*** |
| (0.11) | |
| factor(area)Nord-Est | 0.14 |
| (0.25) | |
| factor(area)Nord-Ovest | -0.50* |
| (0.23) | |
| factor(area)Sud | 1.00*** |
| (0.22) | |
| R2 | 0.08 |
| Adj. R2 | 0.08 |
| Num. obs. | 1806 |
| p < 0.001; p < 0.01; p < 0.05 | |
Importante: se volete che l’output sia renderizzato in HTML, dovete aggiungere l’opzione results = 'asis' al code snippet. Per esempio, il codice per produrre l’output di cui sopra è il seguente:
La funzione htmlreg() ha molte opzioni. Un paio di cose che miglioreranno molto la nostra tabella sono la posizione degli errori standard (attualmente posizionati sotto i coefficienti) e i nomi delle variabili. Possiamo avere gli errori standard di fianco ai coefficienti specificando single.row = TRUE. Se vogliamo che nomi delle variabili siano più interpretabili, dovremo scrivere i nomi che vogliamo in un vettore, e mettere questo vettore nell’opzione custom.coef.names. Inoltre, possiamo decidere cosa scrivere nella didascalia:
# Vettore con i nomi delle variabili
names <- c("Intercetta", "Età", "Sesso (F)", "Bassa istruzione", "Alta istruzione",
"Valutazione dell'economia", "Regione: Nord-Est", "Regione: Nord-Ovest", "Regione: Sud")
# Tabella
htmlreg(m.c, single.row = T,
custom.coef.names = names,
caption = "Modello di regressione per probabilità di voto al M5S")| Model 1 | |
|---|---|
| Intercetta | 5.58 (0.32)*** |
| Età | -0.02 (0.00)*** |
| Sesso (F) | -0.44 (0.16)** |
| Bassa istruzione | 0.01 (0.23) |
| Alta istruzione | -0.64 (0.18)*** |
| Valutazione dell’economia | -0.72 (0.11)*** |
| Regione: Nord-Est | 0.14 (0.25) |
| Regione: Nord-Ovest | -0.50 (0.23)* |
| Regione: Sud | 1.00 (0.22)*** |
| R2 | 0.08 |
| Adj. R2 | 0.08 |
| Num. obs. | 1806 |
| p < 0.001; p < 0.01; p < 0.05 | |
Un punto di forza di texreg è che possiamo anche fare una singola tabella per diversi modelli di regressione. Questo è utile quando vogliamo analizzare più di una variabile dipendente usando le stesse variabili indipendenti. L’unica differenza è che invece di mettere il nome del modello come primo argomento nella funzione htmlreg(), mettiamo una lista con tutti i modelli che vogliamo includere nella tabella.
Per esempio, dato che nei dati CSES che stiamo usando abbiamo le variabili PTV per quattro partiti (M5S, PD, FI e Lega), possiamo stimare quattro modelli di regressione ed esportare i risultati in un’unica tabella (possiamo specificare quali coefficienti si riferiscono a quale partito tramite l’opzione custom.model.names):
m.c.1 <- lm(ptv_M5S ~ age + female + edu_low + edu_high +
eco_eval + factor(area),
data = cses)
m.c.2 <- lm(ptv_PD ~ age + female + edu_low + edu_high +
eco_eval + factor(area),
data = cses)
m.c.3 <- lm(ptv_FI ~ age + female + edu_low + edu_high +
eco_eval + factor(area),
data = cses)
m.c.4 <- lm(ptv_LEGA ~ age + female + edu_low + edu_high +
eco_eval + factor(area),
data = cses)
# Tabella
htmlreg(list(m.c.1, m.c.2, m.c.3, m.c.4),
custom.coef.names = names,
custom.model.names = c("M5S", "PD", "FI", "Lega"),
caption = "Modello di regressione per probabilità di voto a 4 partiti italiani")| M5S | PD | FI | Lega | |
|---|---|---|---|---|
| Intercetta | 5.58*** | 2.38*** | 2.31*** | 3.37*** |
| (0.32) | (0.27) | (0.29) | (0.32) | |
| Età | -0.02*** | 0.02*** | -0.00 | 0.00 |
| (0.00) | (0.00) | (0.00) | (0.00) | |
| Sesso (F) | -0.44** | 0.39** | -0.13 | -0.39* |
| (0.16) | (0.14) | (0.14) | (0.16) | |
| Bassa istruzione | 0.01 | 0.02 | 0.32 | 0.13 |
| (0.23) | (0.19) | (0.21) | (0.23) | |
| Alta istruzione | -0.64*** | 0.41** | -0.57*** | -0.84*** |
| (0.18) | (0.15) | (0.16) | (0.18) | |
| Valutazione dell’economia | -0.72*** | 2.06*** | -0.54*** | -1.14*** |
| (0.11) | (0.09) | (0.10) | (0.10) | |
| Regione: Nord-Est | 0.14 | 0.06 | 0.47* | 0.71** |
| (0.25) | (0.22) | (0.23) | (0.25) | |
| Regione: Nord-Ovest | -0.50* | 0.04 | 0.54** | 0.57* |
| (0.23) | (0.19) | (0.20) | (0.22) | |
| Regione: Sud | 1.00*** | 0.37* | 0.45* | -0.41 |
| (0.22) | (0.19) | (0.20) | (0.22) | |
| R2 | 0.08 | 0.24 | 0.04 | 0.09 |
| Adj. R2 | 0.08 | 0.24 | 0.03 | 0.09 |
| Num. obs. | 1806 | 1813 | 1818 | 1811 |
| p < 0.001; p < 0.01; p < 0.05 | ||||
Salvare la tabella in un file Word
texreg permette anche di esportare la tabella in un file .doc, che può essere aperto con MS Word e condiviso con altre persone senza esportare l’intero documento. Questo viene fatto in due modi. Il primo è utilizzando sempre la funzione htmlreg(), con l’importante differenza che ora specifichiamo anche il nome del file in cui vogliamo salvare la tabella (con l’opzione file) e una serie di altre opzioni che servono ad adattare il codice HTML a Word.
Importante: il nome del file deve avere l’estensione .doc, altrimenti dovremo cambiare estensione per aprirlo con Word. Notate anche che questo non produrrà alcun output nel presente documento:
htmlreg(list(m.c.1, m.c.2, m.c.3, m.c.4),
custom.coef.names = names,
custom.model.names = c("M5S", "PD", "FI", "Lega"),
caption = "Modello di regressione per probabilità di voto a 4 partiti italiani",
file = "PTV_models.doc",
inline.css = FALSE, doctype = TRUE, html.tag = TRUE, head.tag = TRUE, body.tag = TRUE)Un’alternativa per esportare le tabelle in Word è usare la nuova funzione wordreg(). Questa funzione esporta la tabella usando RMarkdown. Tuttavia, la tabella è molto più “basica” rispetto a quella fatta con htmlreg():
texreg()
Se lo script su cui state lavorando deve produrre un file PDF (invece di un file HTML o Word) dovrete usare una funzione diversa per esportare la tabella, la funzione texreg(). La sintassi è la seguente (e produce un codice LaTeX che non viene renderizzato da R Markdown):
texreg(list(m.c.1, m.c.2, m.c.3, m.c.4),
custom.coef.names = names,
custom.model.names = c("M5S", "PD", "FI", "Lega"),
caption = "Modello di regressione per probabilità di voto a 4 partiti italiani")##
## \begin{table}
## \begin{center}
## \begin{tabular}{l c c c c}
## \hline
## & M5S & PD & FI & Lega \\
## \hline
## Intercetta & $5.58^{***}$ & $2.38^{***}$ & $2.31^{***}$ & $3.37^{***}$ \\
## & $(0.32)$ & $(0.27)$ & $(0.29)$ & $(0.32)$ \\
## Età & $-0.02^{***}$ & $0.02^{***}$ & $-0.00$ & $0.00$ \\
## & $(0.00)$ & $(0.00)$ & $(0.00)$ & $(0.00)$ \\
## Sesso (F) & $-0.44^{**}$ & $0.39^{**}$ & $-0.13$ & $-0.39^{*}$ \\
## & $(0.16)$ & $(0.14)$ & $(0.14)$ & $(0.16)$ \\
## Bassa istruzione & $0.01$ & $0.02$ & $0.32$ & $0.13$ \\
## & $(0.23)$ & $(0.19)$ & $(0.21)$ & $(0.23)$ \\
## Alta istruzione & $-0.64^{***}$ & $0.41^{**}$ & $-0.57^{***}$ & $-0.84^{***}$ \\
## & $(0.18)$ & $(0.15)$ & $(0.16)$ & $(0.18)$ \\
## Valutazione dell'economia & $-0.72^{***}$ & $2.06^{***}$ & $-0.54^{***}$ & $-1.14^{***}$ \\
## & $(0.11)$ & $(0.09)$ & $(0.10)$ & $(0.10)$ \\
## Regione: Nord-Est & $0.14$ & $0.06$ & $0.47^{*}$ & $0.71^{**}$ \\
## & $(0.25)$ & $(0.22)$ & $(0.23)$ & $(0.25)$ \\
## Regione: Nord-Ovest & $-0.50^{*}$ & $0.04$ & $0.54^{**}$ & $0.57^{*}$ \\
## & $(0.23)$ & $(0.19)$ & $(0.20)$ & $(0.22)$ \\
## Regione: Sud & $1.00^{***}$ & $0.37^{*}$ & $0.45^{*}$ & $-0.41$ \\
## & $(0.22)$ & $(0.19)$ & $(0.20)$ & $(0.22)$ \\
## \hline
## R$^2$ & $0.08$ & $0.24$ & $0.04$ & $0.09$ \\
## Adj. R$^2$ & $0.08$ & $0.24$ & $0.03$ & $0.09$ \\
## Num. obs. & $1806$ & $1813$ & $1818$ & $1811$ \\
## \hline
## \multicolumn{5}{l}{\scriptsize{$^{***}p<0.001$; $^{**}p<0.01$; $^{*}p<0.05$}}
## \end{tabular}
## \caption{Modello di regressione per probabilità di voto a 4 partiti italiani}
## \label{table:coefficients}
## \end{center}
## \end{table}Se state lavorando su una semplice sessione R (quindi non una sessione R Markdown) potete semplicemente copiare il codice qui sopra e incollarlo in un documento LaTeX. Nel file R_day9_lab_texreg.Rmd c’è un esempio di come la sintassi sopra viene renderizzata da LaTeX in un file PDF.