Današnje predavanje se odnosi na regresijsku analizu, najpopularniji analitički pristup u statistici i data science-u. Modeli koje ćemo spominjati se dominantno koriste u ekonomiji, odnosno ekonometriji no primjenjivost imaju i u drugim društvenim znanostima i analitičko-poslovnim rješenjima. Cilj predavanja je izložiti pregled najvažnijih funkcija i paketa, a ne raspravljati o kakrakteristikama i teoretskim aspektima pojedinih modela.

Software preduvjeti

R paketi

Base R sadržava većinu potrebnih funkcija za osnovnu regresijsku analizu no u današnjem predavanju ćemo koristiti niz dodatnih paketa. Korištenje dodatnih paketa olakšava analizu i omogućuje provođenje sofisticiranijih modela.

  • Novi: broom, estimatr, fixest, sandwich, lmtest, AER, lfe, mfx, margins, modelsummary, vtable
  • Korišteni: tidyverse, hrbrthemes, listviewer

Praktičan način za instalaciju i učitavanje svih potrebnih paketa je izvršavanje donjeg koda. Za pakete broom i modelsummary ćemo koristiti razvojne pakete jer sadržavaju nekoliko funkcionalnosti koje nisu dostupne u trenutnim CRAN verzijama.

## učitaj i instaliraj pakete
if (!require("pacman")) install.packages("pacman")
pacman::p_load(mfx, tidyverse, hrbrthemes, estimatr, fixest, sandwich, lmtest, AER, lfe, margins, vtable)
## razvojne verzije
pacman::p_install_gh("tidymodels/broom") 
pacman::p_install_gh("vincentarelbundock/modelsummary")
##  ggplot2 teme
theme_set(hrbrthemes::theme_ipsum())

Sada kada su potrebni paketi učitani, pogledajmo starwars podatke (već korišteni u prethodnim predavanjima) koje ćemo koristiti u svrhu demonstracije modela:

starwars
## # A tibble: 87 x 14
##    name    height  mass hair_color  skin_color eye_color birth_year sex   gender
##    <chr>    <int> <dbl> <chr>       <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Luke S~    172    77 blond       fair       blue            19   male  mascu~
##  2 C-3PO      167    75 <NA>        gold       yellow         112   none  mascu~
##  3 R2-D2       96    32 <NA>        white, bl~ red             33   none  mascu~
##  4 Darth ~    202   136 none        white      yellow          41.9 male  mascu~
##  5 Leia O~    150    49 brown       light      brown           19   fema~ femin~
##  6 Owen L~    178   120 brown, grey light      blue            52   male  mascu~
##  7 Beru W~    165    75 brown       light      blue            47   fema~ femin~
##  8 R5-D4       97    32 <NA>        white, red red             NA   none  mascu~
##  9 Biggs ~    183    84 black       light      brown           24   male  mascu~
## 10 Obi-Wa~    182    77 auburn, wh~ fair       blue-gray       57   male  mascu~
## # ... with 77 more rows, and 5 more variables: homeworld <chr>, species <chr>,
## #   films <list>, vehicles <list>, starships <list>

Osnove regresijskog modela

lm() funkcija

Glavna naredba za provedbu regresijskog modela u R je lm() funkcija. “lm” je kratica za “linear models”, a funkcijska sintaksa je vrlo intuitivna:

lm(y ~ x1 + x2 + x3 + ..., data = df)

lm() ima izvor podataka kao funkcijski argument (u ovom slučaju hipotetski data frame pod nazivom df). Razlog za to je mogućnost supostojanja različitih objekata u radnom prostoru R pa je potrebno eksplicitno naznačiti koji objekt želimo koristiti za regresiju. To je potrebno čak i ako je df jedini data frame u radnom prostoru R u trenutku izođenja regresije. Alternativna opcija je indeksiranje varijabli:

lm(df$y ~ df$x1 + df$x2 + df$x3 + ...)

Provedimo sada jednostavnu regresiju težine (mass) na visinu (height) koristeći starwars podatke:

ols1 = lm(mass ~ height, data = starwars)
# ols1 = lm(starwars$mass ~ starwars$height) ## alternativno
ols1
## 
## Call:
## lm(formula = mass ~ height, data = starwars)
## 
## Coefficients:
## (Intercept)       height  
##    -13.8103       0.6386

Ovaj regresijski ispis je vrlo jednostavan i sažet zbog toga što se većina korisnih informacija skriva u internoj (list) strukturi “ols1” objekta. U RStudio-u je moguće pregledati tu strukturu pomoću naredbe View(ols1) ili klikom na “ols1” objekt u gornjem desnom panelu. To će otvorriti interaktivni panel u kojem možete detaljnije proučiti ovaj objekt. Taj pristup ne funkcionira u (knit!) R Markdown dokumentu, pa ćemo koristiti listviewer::jsonedit() funkciju za interaktivni pregled:

# View(ols1) ## u uobičajnoj skripti koristi
listviewer::jsonedit(ols1, mode="view") ## za R Markdown

Kao što je vidljivo ols1 objekt sadržava mnoštvo slotova… koji sadržavaju regresijske koeficijente, vektor rezidualnih i fitted (i.e. predicted) vrijednosti, rangove matrice dizajna, input podatke…itd. Za deskriptivni pregled najvažnijih vrijednosti se uobičajno koristi generička summary() funkcija:

summary(ols1) # (slično kao u Stata-i)
## 
## Call:
## lm(formula = mass ~ height, data = starwars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -61.43  -30.03  -21.13  -17.73 1260.06 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.8103   111.1545  -0.124    0.902
## height        0.6386     0.6261   1.020    0.312
## 
## Residual standard error: 169.4 on 57 degrees of freedom
##   (28 observations deleted due to missingness)
## Multiple R-squared:  0.01792,    Adjusted R-squared:  0.0006956 
## F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312

Izvucimo sada regresijske koeficijente:

summary(ols1)$coefficients
##               Estimate  Std. Error    t value  Pr(>|t|)
## (Intercept) -13.810314 111.1545260 -0.1242443 0.9015590
## height        0.638571   0.6260583  1.0199865 0.3120447

Korištenje “tidy” regresijskih koeficijenata iz broom paketa

Iako je najjednostavniji način vađenja koeficijenata summary() funkcija, u praksi je broom (paket) bolji način. broom ima niz korisnih funkcioinalnosti koji regresijske (i druge statističke) objekte pretvaraju u “tidy” data.frame-ove. To je dosta praktično jer se regresijski output često koristi kao input u nekoj drugoj proceduri, npr. za vizualizaciju marginalnih efekata. Sada ćemo pogledati kako pomoću funkcije broom::tidy(..., conf.int = TRUE) možemo prebaciti ols1 regresijski objekt u tidy data.frame koeficijenata i drugih statistika:

# library(broom) ## učitano
tidy(ols1, conf.int = TRUE)
## # A tibble: 2 x 7
##   term        estimate std.error statistic p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)  -13.8     111.       -0.124   0.902 -236.       209.  
## 2 height         0.639     0.626     1.02    0.312   -0.615      1.89

Ove “očišćene” (tidy) koeficijente bismo sada mogli iskoristiti za ggplot2 vizualizaciju, primjerice koristeći geom_pointrange() za prikaz error bar-ova.

broom ima još nekoliko korisnih funkcionalnosti. Npr. broom::glance() daje prikaz modelskih meta podataka (R2, AIC, etc.) u data.frame-u:

glance(ols1)
## # A tibble: 1 x 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1    0.0179      0.000696  169.      1.04   0.312     1  -386.  777.  783.
## # ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Tehnike izvoza podataka u druge formate (npr. LaTeX tablice) ćemo spomenuti pred kraj predavanja.

Regresija na dijelu podataka

Prethodno procijenjeni regresijski model i nije baš nešto… R2 je 0.018. Nizak R2 baca sumnju na ekstremne vrijednosti u podatcima…

Čini se da bi imalo smisla maknuti ekstremnu vrijednost iz regresijskog modela!? To je moguće napraviti na dva načina: 1) napravi novi data.frame i provedi regresiju i 2) subset-aj podatke (direktno) unutar lm() funkcije.

1) Napravi novi data.frame

R dozvoljava mnoštvo objekata u radnom prostoru pa je moguće napraviti novi data.frame objekt koji isključuje ekstremnu vrijednost Jabba. Za to je moguće koroistiti dplyr (predavanje) ili data.table (predavanje). U ovom slučaju ćemo koristiti dplyr jer je to trenutno kompatiblno sa starwars podatcima:

starwars2 =
  starwars %>% 
  filter(name != "Jabba Desilijic Tiure")
  # filter(!(grepl("Jabba", name))) ## Regex također funkcionira
ols2 = lm(mass ~ height, data = starwars2)
summary(ols2)
## 
## Call:
## lm(formula = mass ~ height, data = starwars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.382  -8.212   0.211   3.846  57.327 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -32.54076   12.56053  -2.591   0.0122 *  
## height        0.62136    0.07073   8.785 4.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.14 on 56 degrees of freedom
##   (28 observations deleted due to missingness)
## Multiple R-squared:  0.5795, Adjusted R-squared:  0.572 
## F-statistic: 77.18 on 1 and 56 DF,  p-value: 4.018e-12

2) Subset unutar lm() funkcije

Provedi regresiju i subset-aj direktno u funkcijskom pozivu.

ols2a = lm(mass ~ height, data = starwars %>% filter(!(grepl("Jabba", name))))
summary(ols2a)
## 
## Call:
## lm(formula = mass ~ height, data = starwars %>% filter(!(grepl("Jabba", 
##     name))))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.382  -8.212   0.211   3.846  57.327 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -32.54076   12.56053  -2.591   0.0122 *  
## height        0.62136    0.07073   8.785 4.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.14 on 56 degrees of freedom
##   (28 observations deleted due to missingness)
## Multiple R-squared:  0.5795, Adjusted R-squared:  0.572 
## F-statistic: 77.18 on 1 and 56 DF,  p-value: 4.018e-12

Kvaliteta modela je sada znatno unaprijeđena jer su podaci očišćeni od outliera pa je R2 porastao na 0.58. To ne znači da možemo olako izbacivati podatke…ponekad je moguće ekstremne vrijednosti uzeti u obzir kroz statistički model…primjerice kroz nestandardnu rezidualnu strukturu!

Nestandardna rezidualna struktura

Statističke neregularnosti (heteroskedastičnost, klasteri itd) su regularna stvar u statistici. Dobra stvar je da postoji mnoštvo načina za uvažavanje nestandardne rezudualne strukture u R, npr. sandwich (paket). Nešto “moderniji” pristup je estimatr (paket) koji daje na brzi i praktičnosti. Slijedi nekoliko oglednih primjera…

Robusne standardne pogreške

Heteroskedstične (HC) “robusne” standardne pogreške modela možemo uvažiti korištenjem estimatr::lm_robust() funkcije. Prikažimo to na primjeru ols1 regresijskog objekta koji smo već koristili. estimatr će prikaz dati u “tidy” formatu, no valja imati na umu da je objekt moguće “ugnijezditi” i u tidy() funkciju!

# library(estimatr) ## učitano
ols1_robust = lm_robust(mass ~ height, data = starwars)
# tidy(ols1_robust, conf.int = TRUE) ## tidy() alternativa
ols1_robust
##               Estimate  Std. Error    t value     Pr(>|t|)    CI Lower
## (Intercept) -13.810314 23.45557632 -0.5887859 5.583311e-01 -60.7792950
## height        0.638571  0.08791977  7.2631109 1.159161e-09   0.4625147
##               CI Upper DF
## (Intercept) 33.1586678 57
## height       0.8146273 57

Ovaj paket koristi Eicker-Huber-White robusnu rezidualnu strukturu kao default, odnosno “HC2” standard errors. Default varijantu je moguće promijeniti sa se_type = argumentom 1. To je korisno ako dolazite iz Stata-e i želite replicirati rezultate. Replikacija rezultata nije uvijek elegantna stvar, a za raspravu na temu replikacije između Stata-e, R i Python-a pogledajte ovdje.

lm_robust(mass ~ height, data = starwars, se_type = "stata")
##               Estimate  Std. Error    t value     Pr(>|t|)    CI Lower
## (Intercept) -13.810314 23.36219608 -0.5911394 5.567641e-01 -60.5923043
## height        0.638571  0.08616105  7.4113649 6.561046e-10   0.4660365
##               CI Upper DF
## (Intercept) 32.9716771 57
## height       0.8111055 57

estimatr podržava i robusnu regresiju u modelu instrumentalnih varijabli (IV), no to ćemo detaljnije objasniti kasnije…

Usputni komentar o HAC (Newey-West) rezidualnoj strukturi

estimatr nema podršku za HAC (i.e. heteroskedasticity and autocorrelation consistent) rezidualnu strukturu a la Newey-West. Inicijativa za dodavanje ove funkcionalnosti postoji na GitHub-u, no za sada je HAC omogućen kroz sandwich paket. Npr. korištenje sandwich::NeweyWest()funkcije na ols1 objektu:

# library(sandwich) ## učitano
# NeweyWest(ols1) ## HAC VCOV
sqrt(diag(NeweyWest(ols1))) ## prikaži HAC SEs
## (Intercept)      height 
##  21.2694130   0.0774265

Ako pak želite koristiti HAC SE rezidualnu strukturu u modelu, preporuča se primjena ols1 objekta u lmtest::coeftest()funkciji. Ova funkcija koristi sandwich paket i omogućava praktičan način testiranja hipoteza u modelu sa raznim specifikacijama. Osnovni primjer:

# library(lmtest) ## već učitano
ols1_hac = lmtest::coeftest(ols1, vcov = NeweyWest)
ols1_hac
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -13.810314  21.269413 -0.6493    0.5187    
## height        0.638571   0.077427  8.2474 2.672e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obratite pozornost i na lakoću kojom se coeftest() uklapa u sintaksu broom paketa:

tidy(ols1_hac, conf.int = TRUE)
## # A tibble: 2 x 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)  -13.8     21.3       -0.649 5.19e- 1  -56.4      28.8  
## 2 height         0.639    0.0774     8.25  2.67e-11    0.484     0.794

Klasteri u rezidualnoj strukturi

Klasteri u rezidualnoj strukturi se najčešće javljaju u panel podatcima. Panel podatke ćemo spomenuti kasnije no ovdje je jedan praktični primjer klastera sa estimatr::lm_robust() funkcijom:

lm_robust(mass ~ height, data = starwars, clusters = homeworld)
##               Estimate  Std. Error    t value     Pr(>|t|)    CI Lower
## (Intercept) -9.3014938 28.84436408 -0.3224718 0.7559158751 -76.6200628
## height       0.6134058  0.09911832  6.1886211 0.0002378887   0.3857824
##               CI Upper       DF
## (Intercept) 58.0170751 7.486034
## height       0.8410291 8.195141

Dummie varijable i interakcije varijabli

Za demonstraciju principa u ovom poglavlju je potrebno subset-ati starwars podatke tako da uključuju samo ljudsku vrstu. Zbog toga ćemo napraviti novi podatkovni skup, a dodatno i novu (faktorsku) varijablu “spol” (“gender”) jer je važno pokazati kako R tretira faktorske varijable:

humans = 
  starwars %>% 
  filter(species=="Human") %>%
  mutate(gender_factored = as.factor(gender)) %>% ## faktorska verzija varijable "gender"
  select(contains("gender"), everything())
humans
## # A tibble: 35 x 15
##    gender    gender_factored name   height  mass hair_color skin_color eye_color
##    <chr>     <fct>           <chr>   <int> <dbl> <chr>      <chr>      <chr>    
##  1 masculine masculine       Luke ~    172    77 blond      fair       blue     
##  2 masculine masculine       Darth~    202   136 none       white      yellow   
##  3 feminine  feminine        Leia ~    150    49 brown      light      brown    
##  4 masculine masculine       Owen ~    178   120 brown, gr~ light      blue     
##  5 feminine  feminine        Beru ~    165    75 brown      light      blue     
##  6 masculine masculine       Biggs~    183    84 black      light      brown    
##  7 masculine masculine       Obi-W~    182    77 auburn, w~ fair       blue-gray
##  8 masculine masculine       Anaki~    188    84 blond      fair       blue     
##  9 masculine masculine       Wilhu~    180    NA auburn, g~ fair       blue     
## 10 masculine masculine       Han S~    180    80 brown      fair       brown    
## # ... with 25 more rows, and 7 more variables: birth_year <dbl>, sex <chr>,
## #   homeworld <chr>, species <chr>, films <list>, vehicles <list>,
## #   starships <list>

Dummie varijable kao faktori

Dummie varijable su ključna komponenta velikog broja regresijskih modela, a podrška za rukovanje tim varijablama je dosta loša u nekim (i.e. većini) statističkim programima (npr. potrebna je tabulacija nove matrice binarnih varijabli nakon čega slijedi pripisivanje orginalnim podtacima…). R ima dobar pristup za stvaranje i evaluacijju dummie varijabli: jednostavno specificirajte varijablu kao faktor.2

Ovdje je primjer regresije sa “gendered_factored” varijablom koju smo maločas stvorili:

summary(lm(mass ~ height + gender_factored, data = humans))
## 
## Call:
## lm(formula = mass ~ height + gender_factored, data = humans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.068  -8.130  -3.660   0.702  37.112 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              -84.2520    65.7856  -1.281   0.2157  
## height                     0.8787     0.4075   2.156   0.0441 *
## gender_factoredmasculine  10.7391    13.1968   0.814   0.4259  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.19 on 19 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.444,  Adjusted R-squared:  0.3855 
## F-statistic: 7.587 on 2 and 19 DF,  p-value: 0.003784

Cijela stvar sa prebacivanjem string (i.e. character) varijable u faktorsku nije bila nužna u regresijskom pozivu jer će R to učiniti automatski…no dobro je znati što se događa u pozadini. Pogledajte:

## ne faktorska verzija "gender" varijable; R zna!!!
summary(lm(mass ~ height + gender, data = humans))
## 
## Call:
## lm(formula = mass ~ height + gender, data = humans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.068  -8.130  -3.660   0.702  37.112 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     -84.2520    65.7856  -1.281   0.2157  
## height            0.8787     0.4075   2.156   0.0441 *
## gendermasculine  10.7391    13.1968   0.814   0.4259  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.19 on 19 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.444,  Adjusted R-squared:  0.3855 
## F-statistic: 7.587 on 2 and 19 DF,  p-value: 0.003784

Interakcijski efekti

Kao i kod dummie varijabli, R ima praktičnu sintaksu za specifikaciju interakcijskih varijabli direktno u regresijskom modelu.3 Ovo je standardna sintaksa:

  • x1:x2 “križanje” varijabli (ekvivalentno uključivanju x1 × x2 interakcije)
  • x1/x2 “ugnježđivanje” druge varijable u prvu (ekvivalentno x1 + x1:x2)
  • x1*x2uključuje sve glavne i interakcijske varijable (ekvivalentno x1 + x2 + x1:x2)

Općenito je preporučljivo uključiti sve glavne (parent) odnose uz njihove interakcije pa je * default opcija.

Primjerice, može nas zanimati da li je odnos između težine i visine posredovan spolom osobe!? Tada bismo proveli regresiju u ovakvom obliku:

\[Mass = \beta_0 + \beta_1 D_{Male} + \beta_2 Height + \beta_3 D_{Male} \times Height\]

Za implementaciju u R:

ols_ie = lm(mass ~ gender * height, data = humans)
summary(ols_ie)
## 
## Call:
## lm(formula = mass ~ gender * height, data = humans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.250  -8.158  -3.684  -0.107  37.193 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)            -61.0000   204.0565  -0.299    0.768
## gendermasculine        -15.7224   219.5440  -0.072    0.944
## height                   0.7333     1.2741   0.576    0.572
## gendermasculine:height   0.1629     1.3489   0.121    0.905
## 
## Residual standard error: 15.6 on 18 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.4445, Adjusted R-squared:  0.3519 
## F-statistic: 4.801 on 3 and 18 DF,  p-value: 0.01254

Analiza panel podataka

Fiksni efekti sa fixest paketom

Najjednostavniji način za uključivanje fiksnih efekata u regresijski model je korištenje dummie varijabli, no to je vrlo neučinkovito…Uostalom koji je smisao svih teoretskih aspekata within-group transformacija ako ih ne možemo praktično primijeniti u statističkom programu?! U R-u postoji više opcija za analizu fiksnih efekata: lfe (paket) koji je dosta sličan Stata-inom reghdfe (pristupu) i fixest (paket) koji ćemo koristiti u ovom predavanju.

fixest je relativno novi paket koji ima podršku za analizu nelinearnih modela, visko-dimenzionalnih fiksnih efekata, višestrukih klastera itd., a procedure u paketu su jako brze! (Pogledaj za više o odnosu na lfe i reghdfe pakete.) Ovdje ćemo razmotriti samo najosnovnije principe rada sa fixest paketom, a za detalje pogledajte.

Jednostavni FE model

Glavna funkcija u ovom paketu je fixest::feols() i koristi se za procjenu fiksnih efekata u linearnim modelima. Sintaksa zahtjeva da se prvo specificira “normalni” regresijski model, a nakon | lista fiksnih efekata. Pogledajmo jedan primjer regresije težine (mass) na visinu (height) pri čemu ćemo kontrolirati za fiksne efekte na razini vrste4 ( koristimo starwars podatke!):

# library(fixest) ## učitano
ols_fe = feols(mass ~ height | species, data = starwars) ## fiksni efekti nakon "|"
ols_fe
## OLS estimation, Dep. Var.: mass
## Observations: 58 
## Fixed-effects: species: 31
## Standard-errors: Clustered (species) 
##        Estimate Std. Error t value  Pr(>|t|)    
## height 0.974876   0.044291 22.0105 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 9.69063     Adj. R2: 0.99282 
##                 Within R2: 0.662493

Modelski objekt je automatski klasterirao standardne pogreške po varijabli fiksnih efekata (i.e. species). Za obične (vanilla) reziduale je potrebno specificirati se argument u summary.fixest() funkciji na sljedeći način:

summary(ols_fe, se = 'standard')
## OLS estimation, Dep. Var.: mass
## Observations: 58 
## Fixed-effects: species: 31
## Standard-errors: IID 
##        Estimate Std. Error t value   Pr(>|t|)    
## height 0.974876   0.136463  7.1439 1.3797e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 9.69063     Adj. R2: 0.99282 
##                 Within R2: 0.662493

Podatkovni skup sa koeficijentima ćemo pospremiti u zasebni objekt jer će nam trebati kasnije, a pri tome koristimo vanilla reziduale. Usput primjetite na koji načinbroom::tidy() metoda za fixest objekte prihvaća se argument. Ova procedura je ujedno i jako praktična za provjeru više različitih modela:

# coefs_fe = tidy(summary(ols_fe, se = 'standard'), conf.int = TRUE) ## isto kao dolje
coefs_fe = tidy(ols_fe, se = 'standard', conf.int = TRUE)

Visko-dimenzionalni FE i višetruki klasteri

Kako što smo već naveli, fixest omogućava proizvoljan broj varijabli za fiksne efekte i višestruke klastere (do 4). Pogledajmo kako to izgleda u praksi (uz dodavanje “homeworld” varijable fiksnim efektima):

## Wdvostruki fiksni efekti: species i homeworld
ols_hdfe = feols(mass ~ height |  species + homeworld, data = starwars)
ols_hdfe
## OLS estimation, Dep. Var.: mass
## Observations: 55 
## Fixed-effects: species: 30,  homeworld: 38
## Standard-errors: Clustered (species) 
##        Estimate Std. Error t value Pr(>|t|)    
## height 0.755844   0.332888 2.27057  0.03078 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 7.45791     Adj. R2: 1.00768 
##                 Within R2: 0.487231

Rezidualna struktura prethodnog modela je automatski klasterirana po vrsti (species), odnosno prvoj varijabli iza |. Ukoliko želimo klaster po “species” i “homeworld” varijablama5 možemo koristiti funkcijske argumente se ili cluster u summary.fixest() funkciji. Pripisati ćemo model ols_hdfe objektu:

## klaster po species i homeworld
# ols_hdfe = summary(ols_hdfe, se = 'twoway') ## Isto kao niže
ols_hdfe = summary(ols_hdfe, cluster = c('species', 'homeworld'))
ols_hdfe
## OLS estimation, Dep. Var.: mass
## Observations: 55 
## Fixed-effects: species: 30,  homeworld: 38
## Standard-errors: Clustered (species & homeworld) 
##        Estimate Std. Error t value   Pr(>|t|)    
## height 0.755844   0.116416 6.49263 4.1625e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 7.45791     Adj. R2: 1.00768 
##                 Within R2: 0.487231

Usporedba modelskih koeficijenata

fixest ima ugrađenu coefplot() funkciju za prikaz rezulatata procjene. Ovo je korisno za pregled efekata kroz vrijeme iako (pogledaj!) se procjena koeficijenata za različite modele se najčešće radi sa ggplot2 paketom. Sljedeći primjer se zasniva na funkcionalnosti koju omogućjea spremanje objekta kao data frame-a pomoćubroom::tidy() funkcije. Upravo to olakšava vizualizaciju u ovom slučaju:

# library(ggplot2) ## učitano
## "tidy" output od ols_hdfe objekta
coefs_hdfe = tidy(ols_hdfe, conf.int = TRUE)
bind_rows(
  coefs_fe %>% mutate(reg = "Model 1\nFE bez klastera"),
  coefs_hdfe %>% mutate(reg = "Model 2\nHDFE dvostrani klaster")
  ) %>%
  ggplot(aes(x=reg, y=estimate, ymin=conf.low, ymax=conf.high)) +
  geom_pointrange() +
  labs(Title = "Marginalni efekt visine na težinu") +
  geom_hline(yintercept = 0, col = "orange") +
  ylim(-0.5, NA) + ## Added a bit more bottom space to emphasize the zero line
  labs(
    title = "'Utjecaj' visine na težinu",
    caption = "Podatci: Likovi iz Star Wars univerzuma"
    ) +
  theme(axis.title.x = element_blank())

Normalno bismo očekivali veće standardne pogreške sa klasterskom rezidualnom strukturom, no ovdje je taj efekt poništen kroz povećanu preciznost koju daju fiksni efekti. Ipak ne treba zboraviti da je cjelokupni primjer zasnovan na imaginarnim podatcima pa nema smisla previše razmišljati o substantivnim implikacijama ovih rezultata. Bitna je sintaksa ;-)

Komentar o standardnim pogreškama

Upravo smo vidjeli koje opcije ima fixest pri specifikaciji različitih rezidualnih struktura. Ukratko provedite model sa se ili cluster argumentima u summary.fixest() (ili broom::tidy()) ako niste zaovoljni sa default varijantama. Tu postoje još dvije stvari na koje valja skrenuti pozornost!

Prvo, ako dolazite iz drugog statističkog jezika (Stata!?), prilagodba rezidualne strukture nakon što je proveden model može izgledati neobično…no taj način ima znatne prednosti. Primjerice, on nam omogućava analizu različitih specifikacija po on-the-fly principu bez da ponovno provodimo model. fixest je uistinu brz no vremenski gubitci zbog ponovnog provođenja modela kod većih modela mogu biti znatni.

Drugo, usporedba standardnih pogrešaka u različitim programima je komplicirana stvar. Osim mnoštva nerješenih teoretskih aspekata (posebno kod višestrukih klastera) tu je i problem specifičnosti svakog pojedinog statističkog paketa (Pogledaj za diskusiju!) Usporedba u slučaju fixest paketa je je detaljno opisana u vignette o replikaciji rezidualne strukture kod drugih paketa. 6

Slučajni (random) i mješoviti (mixed) efekti

Fiksni efekti se znatno češće sreću u praktičnim ekonometrijskim analizama nego što je to slučaj sa random ili mixed efektima. Ovdje isto valja spomenuti bayesijanske hijerarhisje modele koji se znatno rijeđe pojavljuju u udžbenicima, a vrlo su zanimljivi. U svakom slučaju, R ima podršku za analizu slučajnih efekata kroz plm (paket) i nlme (paket). 7

Instrumentalne varijable

Kao i kod drugih modela, R pruža više mogućnosti za provedbu IV regresijske analize. Ovdje ćemo razmotriti AER::ivreg(), estimatr::iv_robust(), i lfe::felm() funkcije. Sve ove funkcije imaju sličnu sintaksu, pri čemu je prvi regresijski stupanj specificiran nakon | , a poslije primarnog regresijskog modela. Postoje tu i neke druge suptilne razlike pa odlučite sami koji pristup vam više odgovara. Za demonstraciju ćemo koristiti panel podatke o pušenju po SAD federalnim zemljama iz AER paketa. Prvo ćemo učitati podatke, potom napraviti neke potrebne varijable i onda kratko pogledati kako ti podatci izgledaju. Primjer je ograničen na 1995. godinu jer je cilj razumjeti IV sintaksu, a ne nužno kako ova metodologija funkcionira u kontekstu panel podataka.

## učitaj podatke
data("CigarettesSW", package = "AER")
## stvori novi df sa modificiranim varijablama
cigs =
  CigarettesSW %>%
  mutate(
    rprice = price/cpi,
    rincome = income/population/cpi,
    rtax = tax/cpi,
    tdiff = (taxs - tax)/cpi
    ) %>%
  as_tibble()
## napravi podskup podataka za 1995
cigs95 = cigs %>% filter(year==1995)
cigs95
## # A tibble: 48 x 13
##    state year    cpi population packs    income   tax price  taxs rprice rincome
##    <fct> <fct> <dbl>      <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>
##  1 AL    1995   1.52    4262731 101.   83903280  40.5  158.  41.9   104.    12.9
##  2 AR    1995   1.52    2480121 111.   45995496  55.5  176.  63.9   115.    12.2
##  3 AZ    1995   1.52    4306908  72.0  88870496  65.3  199.  74.8   130.    13.5
##  4 CA    1995   1.52   31493524  56.9 771470144  61    211.  74.8   138.    16.1
##  5 CO    1995   1.52    3738061  82.6  92946544  44    167.  44     110.    16.3
##  6 CT    1995   1.52    3265293  79.5 104315120  74    218.  86.4   143.    21.0
##  7 DE    1995   1.52     718265 124.   18237436  48    166.  48     109.    16.7
##  8 FL    1995   1.52   14185403  93.1 333525344  57.9  188.  68.5   123.    15.4
##  9 GA    1995   1.52    7188538  97.5 159800448  36    157.  37.4   103.    14.6
## 10 IA    1995   1.52    2840860  92.4  60170928  60    191.  69.1   125.    13.9
## # ... with 38 more rows, and 2 more variables: rtax <dbl>, tdiff <dbl>

Pretpostavimo da nas zanima regresijski odnos broja (paketa) cigareta kozumiranih po glavi stanovnika i njihove cijene te osobnog dohotka. Empirijski problem se odnosi na to da je broj konzumiranih paketa endogen, odnosno simultano određen na srani ponude i potražnje. Zbog toga treba koristiti instrumentalnu varijablu koja se u ovom slučaju odnosi na porezne varijable. Ovo je regresijski model koji nas zanima:

\[price_i = \pi_0 + \pi_1 tdiff_i + + \pi_2 rtax_i + v_i \hspace{1cm} \text{(Prva razina)}\] \[packs_i = \beta_0 + \beta_2\widehat{price_i} + \beta_1 rincome_i + u_i \hspace{1cm} \text{(Druga razina)}\]

Opcija 1: AER::ivreg()

Započnimo sa AER::ivreg() funkcijom jer i podatci dolaze iz tog paketa. Prva regresijska razina je specificirana nakon | i uključuje sve egzogene varijable.

# library(AER) ## učitano
## IV regresija 
iv_reg = 
  ivreg(
    log(packs) ~ log(rprice) + log(rincome) | ## glavna regresija; "rprice" je endogena
      log(rincome) + tdiff + rtax, ## lista svih *egzogenih* varijabli uključujući i  "rincome"
    data = cigs95
    )
summary(iv_reg, diagnostics = TRUE)
## 
## Call:
## ivreg(formula = log(packs) ~ log(rprice) + log(rincome) | log(rincome) + 
##     tdiff + rtax, data = cigs95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.6006931 -0.0862222 -0.0009999  0.1164699  0.3734227 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.8950     1.0586   9.348 4.12e-12 ***
## log(rprice)   -1.2774     0.2632  -4.853 1.50e-05 ***
## log(rincome)   0.2804     0.2386   1.175    0.246    
## 
## Diagnostic tests:
##                  df1 df2 statistic p-value    
## Weak instruments   2  44   244.734  <2e-16 ***
## Wu-Hausman         1  44     3.068  0.0868 .  
## Sargan             1  NA     0.333  0.5641    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1879 on 45 degrees of freedom
## Multiple R-Squared: 0.4294,  Adjusted R-squared: 0.4041 
## Wald test: 13.28 on 2 and 45 DF,  p-value: 2.931e-05

Za ekonomiste naviknute na Stata-u će ovo biti u najboljem slučaju neintuitivno.8 U ovom slučaju nismo specificirali endogene varijable (i.e. “rplice”) direktno nego smo rekli R-u koje su egzogene varijable. R je nakon toga zaključio gdje endogene varijable trebaju biti instrumentalizirane i proveo prvu regresijsku razinu u pozadini. Ovaj modelski set-up dobiva još smisla ako promislite o teoretskim osnovama IV pristupa!

AER::ivreg() također omogućava alternativni način specifikacije prve regresijske razine. Sada ćemo označiti endogenu “rprice” varijablu sa. -price i uključiti samo instrumentalne varijable nakon |. Output je jednak kao i u prethodnom slučaju:

## eksplicitna specifikacija instrumenata
ivreg(
  log(packs) ~ log(rprice) + log(rincome) | 
    . -log(rprice) + tdiff + rtax, ## alternativni način specifikacije prve regresijske razine
  data = cigs95
  )

Opcija 2: estimatr::iv_robust()

Drugi način se odnosi na estimatr koji smo prethodno spomenuli. Default postavka u ovom pristupu je HC2 rezidualna struktura, a naravno da su dopuštene i druge opcije (i klasteri također). Sintaksa je skoro ista kao u prethodnom primjeru, a sve što treba promijeniti je funkcijski poziv iz AER::ivreg() u estimatr::iv_robust().

# library(estimatr) ## učitano
## IV regresija sa robusnim SE
iv_reg_robust = 
  iv_robust( ## sve je isto osim funkcijskog poziva
    log(packs) ~ log(rprice) + log(rincome) | 
      log(rincome) + tdiff + rtax,
    data = cigs95
    )
summary(iv_reg_robust, diagnostics = TRUE)
## 
## Call:
## iv_robust(formula = log(packs) ~ log(rprice) + log(rincome) | 
##     log(rincome) + tdiff + rtax, data = cigs95)
## 
## Standard error type:  HC2 
## 
## Coefficients:
##              Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
## (Intercept)    9.8950     0.9777  10.120 3.569e-13   7.9257  11.8642 45
## log(rprice)   -1.2774     0.2547  -5.015 8.739e-06  -1.7904  -0.7644 45
## log(rincome)   0.2804     0.2547   1.101 2.768e-01  -0.2326   0.7934 45
## 
## Multiple R-squared:  0.4294 ,    Adjusted R-squared:  0.4041 
## F-statistic:  15.5 on 2 and 45 DF,  p-value: 7.55e-06

Opcija 3: felm::lfe()

felm() funkcija iz lfe paketa je vjerojatno najelegantnija opcija u usporedbi sa pretthodne dvije jer ima najintuitivniju sintaksu.9 Njena sintaksa je jako slična Stata-inom načinu specifikacije prve regresijske razine, gdje se navode samo endogene varijable i instrumenti.

# library(lfe) ## učitano
iv_felm = 
  felm(
    log(packs) ~ log(rincome) |
      0 | ## bez fiksnih efekata 
      (log(rprice) ~ tdiff + rtax), ## prva razina; obrati pažnju na zagrade
    data = cigs95
  )
summary(iv_felm)
## 
## Call:
##    felm(formula = log(packs) ~ log(rincome) | 0 | (log(rprice) ~      tdiff + rtax), data = cigs95) 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60069 -0.08622 -0.00100  0.11647  0.37342 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.8950     1.0586   9.348 4.12e-12 ***
## log(rincome)         0.2804     0.2386   1.175    0.246    
## `log(rprice)(fit)`  -1.2774     0.2632  -4.853 1.50e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1879 on 45 degrees of freedom
## Multiple R-squared(full model): 0.4294   Adjusted R-squared: 0.4041 
## Multiple R-squared(proj model): 0.4294   Adjusted R-squared: 0.4041 
## F-statistic(full model):13.28 on 2 and 45 DF, p-value: 2.931e-05 
## F-statistic(proj model): 13.28 on 2 and 45 DF, p-value: 2.931e-05 
## F-statistic(endog. vars):23.56 on 1 and 45 DF, p-value: 1.496e-05

U gornjem primjeru smo unijeli “0” na mjesto fiksnih efekata jer koristimo samo subset podataka. Sljedeći primjer se odnosi na IV regresiju sa felm()funkcijom pri čemu su uključeni svi podatci (i.e. puni panel) te “year” and “state” fiksni efekti za kontrolu panel strukture.

iv_felm_all = 
  felm(
    log(packs) ~ log(rincome) |
      year + state | ## uključi fiksne efekte
      (log(rprice) ~ tdiff + rtax), 
    data = cigs ## puni panel
  )
summary(iv_felm_all)
## 
## Call:
##    felm(formula = log(packs) ~ log(rincome) | year + state | (log(rprice) ~      tdiff + rtax), data = cigs) 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.08393 -0.03851  0.00000  0.03851  0.08393 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## log(rincome)         0.4620     0.3081   1.500    0.141    
## `log(rprice)(fit)`  -1.2024     0.1712  -7.024  9.4e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06453 on 45 degrees of freedom
## Multiple R-squared(full model): 0.9668   Adjusted R-squared: 0.9299 
## Multiple R-squared(proj model): 0.5466   Adjusted R-squared: 0.04281 
## F-statistic(full model):26.21 on 50 and 45 DF, p-value: < 2.2e-16 
## F-statistic(proj model): 27.71 on 2 and 45 DF, p-value: 1.436e-08 
## F-statistic(endog. vars):49.33 on 1 and 45 DF, p-value: 9.399e-09

Drugi modeli

Generalizirani linearni modeli (logit, etc.)

Za provedbu generaliziranih linearnih modela (GLM) R ima ugrađenu (base) glm() funkciju u kojoj je potrebno specificirati neku od opcija koje opisuju rezidualnu strukturu i željeni model. Ovo je primjer za logit model:

glm_logit = glm(am ~ cyl + hp + wt, data = mtcars, family = binomial)
tidy(glm_logit, conf.int = TRUE)
## # A tibble: 4 x 7
##   term        estimate std.error statistic p.value  conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>     <dbl>     <dbl>
## 1 (Intercept)  19.7       8.12       2.43   0.0152   8.56      44.3   
## 2 cyl           0.488     1.07       0.455  0.649   -1.53       3.12  
## 3 hp            0.0326    0.0189     1.73   0.0840   0.00332    0.0884
## 4 wt           -9.15      4.15      -2.20   0.0276 -21.4       -3.48

Prije nego što pogledamo kako izvaditi marginalne efekte fiz nelinearnih modela…obratite pozornost na mfx paket pomoću kojeg je moguće dobiti marginalne efekte iz niza GLM modela. Npr.:

# library(mfx) ## učitano
## Oprez: mfx učitava MASS paket pa se javlja "namespace conflict"
## sa dplyr-ovom select() funkcijom. Eksplicitno definirajte što želite koristiti:
## e.g. `select = dplyr::select`
## uzmi marginalne efekte
glm_logitmfx = logitmfx(glm_logit, atmean = TRUE, data = mtcars)
## moguć je i unos u funkciju direktno
# glm_logitmfx = logitmfx(am ~ cyl + hp + wt, atmean = TRUE, data = mtcars)
tidy(glm_logitmfx, conf.int = TRUE)
## # A tibble: 3 x 8
##   term  atmean estimate std.error statistic p.value conf.low conf.high
##   <chr> <lgl>     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 cyl   TRUE    0.0538    0.113       0.475   0.635 -0.178     0.286  
## 2 hp    TRUE    0.00359   0.00290     1.24    0.216 -0.00236   0.00954
## 3 wt    TRUE   -1.01      0.668      -1.51    0.131 -2.38      0.359

Bayes-ova regresija

Bayes-ijanski pristup statistici je vrlo opširna tema, a ovo je jako kratki prikaz mogućnosti za provedu te vrste analize u R. Sučelja za provedbu bayes-ijanskih modela su implementirana kroz MCMC i Bayesian software engines: Stan, JAGS, TensorFlow (via Greta), itd. Ovo je samo jedan jednostavni primjer bayes-ijanske analize sa rstanarm paketom. Primijetite da ovaj paket nismo instalirali na početku jer instalacija često zna izazvati svakakve probleme.10

# install.packages("rstanarm") ## izvršite ovo za početak
library(rstanarm)
bayes_reg = 
  stan_glm(
    mass ~ gender * height,
    data = humans, 
    family = gaussian(), prior = cauchy(), prior_intercept = cauchy()
    )
summary(bayes_reg)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      mass ~ gender * height
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 22
##  predictors:   4
## 
## Estimates:
##                          mean   sd     10%    50%    90% 
## (Intercept)             -65.7   76.1 -161.2  -65.6   27.7
## gendermasculine          -0.5    9.1   -7.3   -0.1    5.9
## height                    0.8    0.5    0.2    0.8    1.4
## gendermasculine:height    0.1    0.1    0.0    0.1    0.2
## sigma                    15.9    2.7   12.7   15.6   19.4
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 82.5    4.9 76.3  82.5  88.8 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                        mcse Rhat n_eff
## (Intercept)            1.6  1.0  2183 
## gendermasculine        0.3  1.0  1227 
## height                 0.0  1.0  2084 
## gendermasculine:height 0.0  1.0  1584 
## sigma                  0.1  1.0  2447 
## mean_PPD               0.1  1.0  3334 
## log-posterior          0.0  1.0  1376 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Još neki modeli

Postoji previše modela i empirijskih procedura da bismo ih sve pokrili u jednom predavanju. Veliki dio tih modela dolazi po default-u u osnovnoj R instalaciji. Ovdje su istaknuti neki novi paketi za specifične modele:

  • Difference-in-differences (sa varijabilnim T (vremenska dimenzija) itd): did (link) i DRDID (link)
  • Sintetička kontrola: gsynth (link) i scul (link)
  • Count modeli (hurdle modeli, itd.): pscl (link)
  • Lasso: biglasso (link)
  • Causal forests: grf (link)
  • itd

Na poslijetku vrijedi pogledati i neke korisne resurse (knjige i tutoriale) i linkove za ekonometriju na kraju ovog predavanja.

Marginalni efekti

Izračun marginalnih efekata u regresiji je jednostavan u slučaju kada nema nelinearnosti u modelu…dovoljno je pogledati vrijednosti koeficijenata u regresijskom ispisu. Za nelinearne modele poput logit, probit i sl., marginalne efekte je potrebno izračunati…a za to postoje funkcije u R. Već smo spominjali mfx paket za izračun marginalnih efekata u GLM modelima no ovdje ćemo istaknuti dvije metode za izračun marginalnih efekata kod više različitih vrsta modela: 1) margins paket 2) pristup za analizu svih modela sa interakcijskim varijablama.

margins paket

margins (paket) je izvrstan način za rad sa marginalnim efektima u različitim modelima.11 Za detalje vrijedi pročitati vignette, a sada ćemo prikazati ogledni primjer.

U prethodnom regresijskom primjeru smo razmatrali odnos težine vs. visine i spola ljudi…Da bismo vidjeli prosječne marginalne efekte (average marginal effect (AME)) ovih zavisnih varijabli, moguće je koristiti margins::margins() funkciju:

# library(margins) ## učitano
ols_ie_marg = margins(ols_ie)

Kao i kod običnog regresijskog modela, prikaz procjene je moguće pregledati na sljedeći način:

# summary(ols_ie_marg) ## Same effect
tidy(ols_ie_marg, conf.int = TRUE)

Za usporedbu marginalnih efekata na specifičnim vrijednostima (npr. odnos AME visine i težine po spolu) možemo napraviti sljedeće:

ols_ie %>% 
  margins(
    variables = "height", ## glavna varijabla od interesa
    at = list(gender = c("masculine", "feminine")) ## na kojim specifičnim vrijednostima
    ) #%>% 
## Average marginal effects at specified values
## lm(formula = mass ~ gender * height, data = humans)
##  at(gender) height
##   masculine 0.8962
##    feminine 0.7333
#  tidy(conf.int = TRUE) ## očisti

Za vizualizaciju je moguće koristiti margins::cplot() funkciju koja omogućava prikaz uvjetnih marginalnih efekata:

cplot(ols_ie, x = "gender", dx = "height", what = "effect")

Ovo je samo demonstrativni prikaz bez puno substantivnog sadržaja!

cplot() funkciju je moguće koristiti za prikaz predviđenih (i.e. predicted) vrijednosti zavisne varijable (ovdje: “mass”) uvjetno po nekoj od nezavisnih varijabli:

par(mfrow=c(1, 2)) ## definiraj grid za prikaz grafika
cplot(ols_ie, x = "gender", what = "prediction")
cplot(ols_ie, x = "height", what = "prediction")

par(mfrow=c(1, 1)) ## resetiraj grid

cplot() koristi base R sustav za izradu grafikona pa zbog toga je potrebno definirati grid za grafikone pomoću par() funkcije. ggplot2 podrška za marginsplot paket je opisana ovdje.

Valja još spomenti i emmeans paket,koji je dosta sličan margins paketu.

Specijalni slučaj: / kratica za interakcije među varijablama

Korištenjem / operatora je moguće specificirati puni skup marginalnih efekata za interakcijske varijable. Princip je da standardnu interakcijsku specifikaciju f1 * x2 možemo zamijeniti sa f1 / x2 i automatski ćemo dobiti puni skup marginalnih efekata. Formalnije rečeno, varijable u modelu su ugniježdene (engl. nested).

Ovo je jedan brzi primjer opisanog principa:.

# ols_ie = lm(mass ~ gender * height, data = humans) ## originalni model
ols_ie_marg2 = lm(mass ~ gender / height, data = humans)
tidy(ols_ie_marg2, conf.int = TRUE)
## # A tibble: 4 x 7
##   term                   estimate std.error statistic p.value conf.low conf.high
##   <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)             -61.0     204.      -0.299   0.768  -4.90e+2    368.  
## 2 gendermasculine         -15.7     220.      -0.0716  0.944  -4.77e+2    446.  
## 3 genderfeminine:height     0.733     1.27     0.576   0.572  -1.94e+0      3.41
## 4 gendermasculine:height    0.896     0.443    2.02    0.0582 -3.46e-2      1.83

Mrginalni efekti gender × height interakcija (i.e. 0.733 i 0.896) su isti kao kada smo koristili margins::margins() funbkciju u prethodnom slučaju.

Ovaj je pristup posebno koristan kod velilkih modela sa puno vraijabli. margins paket koristi numeričke (delta) metode koje su memorijski zahtjevne, a korištenje / ne opterećuje dodatno memoriju.

Prezentacija rezultata

Tablice

Regresijske tablice

U R postoji uistinu mnogo različitih opcija za izradu regresijskih tablica.12 Izvrsna opcija za tablice je modelsummary paket za izradu i izvoz regresijskih tablica. Ovaj paket je fleksibilan i podražva čitav niz modelskih outputa. modelsummary također podržava vizualizaciju koeficijenata i deskriptivnu statistiku podataka. Dokumentacija za paket je također izvrsna!

# library(modelsummary) ## Učitano
## Hint: msummary() je kratica za modelsummary()
msummary(list(ols1, ols_ie, ols_fe, ols_hdfe))
Model 1 Model 2 Model 3 Model 4
(Intercept) -13.810 -61.000
(111.155) (204.057)
height 0.639 0.733 0.975 0.756
(0.626) (1.274) (0.044) (0.116)
gendermasculine -15.722
(219.544)
gendermasculine × height 0.163
(1.349)
Num.Obs. 59 22 58 55
R2 0.018 0.444 0.997 0.998
R2 Adj. 0.001 0.352 0.993 1.008
R2 Within 0.662 0.487
R2 Pseudo
AIC 777.0 188.9 492.1 513.1
BIC 783.2 194.4 558.0 649.6
Log.Lik. -385.503 -89.456 -214.026 -188.552
F 1.040 4.801
Std.Errors by: species by: species & homeworld
FE: species X X
FE: homeworld X


Još jedna sjajna stvar vezano uz modelsummary je kompatibilnost sa R Markdown-om. Funkcija će automatski prilagoditi tablice željenom output-u: HTML, LaTeX/PDF, RTF, itd. Vrstu ouputa je moguće i specificirati prema potrebi ako ne koristite R Markdown, a želite izvesti tablicu u neki specifičan format.

Tablice za deskriptivnu statistiku

modelsummary::datasummary*() funkcija/e omogućavaju mnoštvo tablica za deskriptivnu statistiku,a za puni popis funkcionalnosti pogledajte dokumentaciju. Ovdje je jednostavan primjer na subset-u prethodno korištenog “humans” data frame-a.

datasummary_balance(~ gender,
                    data = select(humans, gender, height, mass, birth_year, eye_color))
feminine (N=9)
masculine (N=26)
Mean Std. Dev. Mean Std. Dev. Diff. in Means Std. Error
height 160.2 7.0 182.3 8.2 22.1 3.0
mass 56.3 16.3 87.0 16.5 30.6 10.1
birth_year 46.4 18.8 55.2 26.0 8.8 10.2
N Pct. N Pct.
eye_color blue 3 33.3 9 34.6
blue-gray 0 0.0 1 3.8
brown 5 55.6 12 46.2
dark 0 0.0 1 3.8
hazel 1 11.1 1 3.8
yellow 0 0.0 2 7.7


Slične mogućnosti nudi i vtable paket. Ovaj paket ima mnoštvo sličnosti sa Stata-om pa će biti posebno pogodan za tablice deskriptivne statsitike u ekonomiji. Ovo je ekvivalent prethodne tablice:

# library(vtable) ## učitano
## st() je kratica za sumtable()
st(select(humans, gender, height, mass, birth_year, eye_color), 
   group = 'gender')
Summary Statistics
gender
feminine
masculine
Variable N Mean SD N Mean SD
height 8 160.25 6.985 23 182.348 8.189
mass 3 56.333 16.289 19 86.958 16.549
birth_year 5 46.4 18.77 20 55.165 26.02
eye_color 9 26
… blue 3 33.3% 9 34.6%
… blue-gray 0 0% 1 3.8%
… brown 5 55.6% 12 46.2%
… dark 0 0% 1 3.8%
… hazel 1 11.1% 1 3.8%
… yellow 0 0% 2 7.7%


vtable::st() funkcija je pametna utoliko što automatski izbire bitne varijable i izbacije nebitne (npr. faktori sa puno razina). Ovdje je primjer za “starwars” data frame koji smo već koristili:

st(starwars)
Summary Statistics
Variable N Mean Std. Dev. Min Pctl. 25 Pctl. 75 Max
height 81 174.358 34.77 66 167 191 264
mass 59 97.312 169.457 15 55.6 84.5 1358
birth_year 43 87.565 154.691 8 35 72 896
sex 83
… female 16 19.3%
… hermaphroditic 1 1.2%
… male 60 72.3%
… none 6 7.2%
gender 83
… feminine 17 20.5%
… masculine 66 79.5%

Grafikoni

Prikazi koefficijenata

Već smo vidjeli kako izraditi i prikazati koeficijente procijenjenih modela. To je (u principu) izvrstan način za brzi pregled koeficijenata modela. Ovdje valja istaknuti još neke funkcionalnosti koje modelsummary paket pruža, primjerice modelplot() funkciju. Posebno je korisno što modelplot() prihvaća i ggplot2 sintaksu.

# library(modelsummary) ## učitano
mods = list('FE, bez klastera' = summary(ols_fe, se = 'standard'),  # bez SE klastera 
            'HDFE, dvostrani klaster' = ols_hdfe)
modelplot(mods) +
  ## radi i sa drugim ggplot2 naredbama
  coord_flip() + 
  labs(
    title = "'Utjecaj' visine na težinu",
    subtitle = "Usporedba modela fiksnih efekata"
    )

Ovdje je još jedan primjer za Masculine × Height koeficijente iz prethodnog interakcijskog modela sa punim marginalnim efektima:

ie_mods = list('Djelomični efekti' = ols_ie, 'Marginalni efekti' = ols_ie_marg2)
modelplot(ie_mods, coef_map = c("gendermasculine:height" = "Masculine × Height")) +
  coord_flip() + 
  labs(
    title = "Utjecaj' visine na težinu",
    subtitle = "Djelomični vs marginalni efekti"
    )

Predikcija i validacija modela

Najjednostavniji način za vizualnu provjeru modela (i.e. validacija i predikcija) je ggplot2. Tu je posebno pogodna geom_smooth() funkcija koju smo već spominjali, a koja omogućuje mapiranje modela na grafikon. Na osnovistarwars2 data.frame-a, bez outlier-a, bi to izgledalo otprilike ovako :

ggplot(starwars2, aes(x = height, y = mass)) + 
    geom_point(alpha = 0.7) +
    geom_smooth(method = "lm") ## ?geom_smooth za druge opcije
## `geom_smooth()` using formula 'y ~ x'

Ovdje je važno spomenuti predict() funkciju iz base R. Zamislimo situaciju u kojoj procjenjujemo bivarijatni regresijski odnos visine i težine, kao u prethodim slučajevima…13 Ovaj put ćemo procijeniti model na training podatcima koji se satoje od 30 likova rangiranih po visini. Procedura izgleda ovako:

## procijeni model na training uzorkua (30 najnižih)
ols1_train = lm(mass ~ height, data = starwars %>% filter(rank(height) <=30))
## predvidi težinu svih 30 starwars likova 
## uključi 95% intervale pouzdanosti;vidi ?predict.lm za druge intervale
##i opcije
predict(ols1_train, newdata = starwars2, interval = "prediction") %>%
  head(5) ## prvih nekoliko redova
##        fit       lwr       upr
## 1 68.00019 46.307267  89.69311
## 2 65.55178 43.966301  87.13725
## 3 30.78434  8.791601  52.77708
## 4 82.69065 60.001764 105.37954
## 5 57.22718 35.874679  78.57968

Prthodno korišteni data.frame se može jednostavno kombinirati sa originalnim podatcima u ggplot2 vizualizaciji. Za tu svrhu je posebno koristan broom paket, odnosno njegova augment() funkcija koja omogućava dodavanje modelskih predviđanja podatkovnom skupu. augment() funkcija prihvaća iste argumente kao i predict() funkcija, naravno, uz neke manje izmjene.14

## alternativa za predict(): koristi augment() za dodavanje .fitted i .resid, kao i za 
## .conf.low i .conf.high predikcijskih intervala podatcima.
starwars2 = augment(ols1_train, newdata = starwars2, interval = "prediction")
## prikaži nove varijable (sve imaju a "." prefiksx)
starwars2 %>% select(contains("."), everything()) %>% head()
## # A tibble: 6 x 18
##   .fitted .lower .upper .resid name           height  mass hair_color skin_color
##     <dbl>  <dbl>  <dbl>  <dbl> <chr>           <int> <dbl> <chr>      <chr>     
## 1    68.0  46.3    89.7   9.00 Luke Skywalker    172    77 blond      fair      
## 2    65.6  44.0    87.1   9.45 C-3PO             167    75 <NA>       gold      
## 3    30.8   8.79   52.8   1.22 R2-D2              96    32 <NA>       white, bl~
## 4    82.7  60.0   105.   53.3  Darth Vader       202   136 none       white     
## 5    57.2  35.9    78.6  -8.23 Leia Organa       150    49 brown      light     
## 6    70.9  49.1    92.8  49.1  Owen Lars         178   120 brown, gr~ light     
## # ... with 9 more variables: eye_color <chr>, birth_year <dbl>, sex <chr>,
## #   gender <chr>, homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Sada možemo vidjeti rezulatate modela na 30 najnižih starwaars likova u usporedbi sa punim podatkovnim skupom.

starwars2 %>%
  ggplot(aes(x = height, y = mass, col = rank(height)<=30, fill = rank(height)<=30)) +
  geom_point(alpha = 0.7) +
  geom_line(aes(y = .fitted)) +
  scale_color_discrete(name = "Training uzorak?", aesthetics = c("colour", "fill")) +
  labs(
    title = "Predviđanje težina na osnovi visine",
    caption = "Regresijska linija 95% intervali pouzdanosti."
    ) +
  xlab("Visina") +
  ylab("Težina")

Dodatni resursi (uglavnom za ekonometriju)


  1. Pogledaj za detalje!↩︎

  2. Faktori su varijable sa jedinstvenim kvalitativnim razinama, npr. “muški”, “ženski”, “trans”, itd.↩︎

  3. Iako postoji mnoštvo razloga da to napravite u zasebnom koraku (npr. standardizacija).↩︎

  4. Pošto smo specificirali “species” u slotu za fiksne efekte, feols() funkcija će automatski pretvoriti potrebne varijable u faktorske.↩︎

  5. Ovo je samo demonstrativni primjer.↩︎

  6. Za detalje vrijedi pročitati!↩︎

  7. plm paket također podržava FE (i pooling) modele. Ipak, čini mi se da su fixest i lfe praktičniji.↩︎

  8. Uz pretpostavku da ste veće napravili logaritamske varijable i subset-irali podatke, naredba u Stata-i bi izgledala otprilike ovako: ivreg log_packs = log_rincome (log_rprice = tdiff rtax).↩︎

  9. fixest paket koji smo razmatrali maločas ne podržava IV regresiju. Sa druge strane, lfe::felm()funkcija ima skore svu funkcinalnost tog paketa, jednio što je malo sporija.↩︎

  10. Primjerice na ovom računalu je bilo potrebno instalirati stan i rstanarm pri čemu je R stalno ispadao kroz instalacijski proces koji je trebalo ponavljati više puta.↩︎

  11. Ovaj pristup ipak ne podržava fixest (niti lfe) modele. Neka alternativna rješenja ipak postoje.↩︎

  12. fixest paket ima etable() funkciju. Ovo je optimizirana funkcija za proizvodnju sjajnih tablica uz minimalni trud no ograničena je na objekte koje proizvodifixest funkcija. Za detalje pogledajte ovdje i ovdje.↩︎

  13. Bivarijatna regresija omogućava izradu 2D vizualizacije.↩︎

  14. Ovdje dodajemo “.fitted”, “.resid”, “.conf.low”, i “.conf.high” varijable (i.e. kolone) data.frame-u. Konvencija u augment() funkciji je prefiksiranje dodanih varijabli sa “.” kako bi se izbjegli konflikti sa postojećim varijablama.↩︎