• 論文やレポートに分析結果を綺麗に出力したい

• 複数の回帰結果をまとめて扱いたい

library(tidyverse)
library(gapminder)

gapminder = gapminder %>% 
    mutate(pop = as.double(pop))

gapminder

## # A tibble: 1,704 x 6
##        country continent  year lifeExp      pop gdpPercap
##         <fctr>    <fctr> <int>   <dbl>    <dbl>     <dbl>
##  1 Afghanistan      Asia  1952  28.801  8425333  779.4453
##  2 Afghanistan      Asia  1957  30.332  9240934  820.8530
##  3 Afghanistan      Asia  1962  31.997 10267083  853.1007
##  4 Afghanistan      Asia  1967  34.020 11537966  836.1971
##  5 Afghanistan      Asia  1972  36.088 13079460  739.9811
##  6 Afghanistan      Asia  1977  38.438 14880372  786.1134
##  7 Afghanistan      Asia  1982  39.854 12881816  978.0114
##  8 Afghanistan      Asia  1987  40.822 13867957  852.3959
##  9 Afghanistan      Asia  1992  41.674 16317921  649.3414
## 10 Afghanistan      Asia  1997  41.763 22227415  635.3414
## # ... with 1,694 more rows

• rmdでhtmlとして出力するときはresults='asis'にする
• type = "latex"とするとtex形式で，"text"とするとコンソール上で見やすい形に出力できる

library(stargazer)

gapminder %>% 
    lm(lifeExp ~ year + pop + gdpPercap, data = .) %>% 
    stargazer(type = "html")

gapminder = gapminder %>% 
    mutate_if(is.factor, as.character)

by_country = gapminder %>% 
    group_by(continent, country) %>% 
    nest()

by_country

## # A tibble: 142 x 3
##    continent     country              data
##        <chr>       <chr>            <list>
##  1      Asia Afghanistan <tibble [12 x 4]>
##  2    Europe     Albania <tibble [12 x 4]>
##  3    Africa     Algeria <tibble [12 x 4]>
##  4    Africa      Angola <tibble [12 x 4]>
##  5  Americas   Argentina <tibble [12 x 4]>
##  6   Oceania   Australia <tibble [12 x 4]>
##  7    Europe     Austria <tibble [12 x 4]>
##  8      Asia     Bahrain <tibble [12 x 4]>
##  9      Asia  Bangladesh <tibble [12 x 4]>
## 10    Europe     Belgium <tibble [12 x 4]>
## # ... with 132 more rows

library(sandwich)

lin_reg = function(data) {glm(lifeExp ~ year + pop + gdpPercap, data = data)}
het_se = function(reg) {reg %>% vcovHC() %>% sqrt() %>% diag()}

by_country_with_reg = by_country %>% 
    mutate(reg = map(data, lin_reg), se = map(reg, het_se))

by_country_with_reg

## # A tibble: 142 x 5
##    continent     country              data       reg        se
##        <chr>       <chr>            <list>    <list>    <list>
##  1      Asia Afghanistan <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  2    Europe     Albania <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  3    Africa     Algeria <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  4    Africa      Angola <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  5  Americas   Argentina <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  6   Oceania   Australia <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  7    Europe     Austria <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  8      Asia     Bahrain <tibble [12 x 4]> <S3: glm> <dbl [4]>
##  9      Asia  Bangladesh <tibble [12 x 4]> <S3: glm> <dbl [4]>
## 10    Europe     Belgium <tibble [12 x 4]> <S3: glm> <dbl [4]>
## # ... with 132 more rows

library(magrittr)

by_country_with_reg %>% 
    slice(1:8) %$% 
    stargazer(reg, se = se, type = "html")

by_country_with_reg %>% 
    slice(1:8) %$% 
    stargazer(reg, se = se, type = "html",
        covariate.labels = c("Year", "Population", "GdpPerCap"),
        omit.stat = c("ll", "aic"), omit.table.layout = "bl",
        add.lines = list(c("Continent", continent), c("Country", country)))

by_continent = gapminder %>% 
    group_by(continent) %>% 
    nest() %>% 
    add_column(country = "All")

by_country_with_continent = by_continent %>% 
    bind_rows(by_country) %>% 
    mutate(reg = map(data, lin_reg), se = map(reg, het_se))

by_country_with_continent

## # A tibble: 147 x 5
##    continent               data     country       reg        se
##        <chr>             <list>       <chr>    <list>    <list>
##  1      Asia <tibble [396 x 5]>         All <S3: glm> <dbl [4]>
##  2    Europe <tibble [360 x 5]>         All <S3: glm> <dbl [4]>
##  3    Africa <tibble [624 x 5]>         All <S3: glm> <dbl [4]>
##  4  Americas <tibble [300 x 5]>         All <S3: glm> <dbl [4]>
##  5   Oceania  <tibble [24 x 5]>         All <S3: glm> <dbl [4]>
##  6      Asia  <tibble [12 x 4]> Afghanistan <S3: glm> <dbl [4]>
##  7    Europe  <tibble [12 x 4]>     Albania <S3: glm> <dbl [4]>
##  8    Africa  <tibble [12 x 4]>     Algeria <S3: glm> <dbl [4]>
##  9    Africa  <tibble [12 x 4]>      Angola <S3: glm> <dbl [4]>
## 10  Americas  <tibble [12 x 4]>   Argentina <S3: glm> <dbl [4]>
## # ... with 137 more rows

by_country_with_continent %>% 
    slice(1:8) %$% 
    stargazer(reg, se = se, type = "html",
        covariate.labels = c("Year", "Population", "GdpPerCap"),
        omit.stat = c("ll", "aic"), omit.table.layout = "bl",
        add.lines = list(c("Continent", continent), c("Country", country)))

lin_reg2 = function(data) {glm(lifeExp ~ year*pop*gdpPercap, data = data)}

by_continent_with_reg = by_continent %>% 
    mutate(reg1 = map(data, lin_reg), reg2 = map(data, lin_reg2)) %>% 
    gather("model", "reg", starts_with("reg")) %>% 
    mutate(se = map(reg, het_se)) %>% 
    arrange(continent, model)

by_continent_with_reg

## # A tibble: 10 x 6
##    continent               data country model       reg        se
##        <chr>             <list>   <chr> <chr>    <list>    <list>
##  1    Africa <tibble [624 x 5]>     All  reg1 <S3: glm> <dbl [4]>
##  2    Africa <tibble [624 x 5]>     All  reg2 <S3: glm> <dbl [8]>
##  3  Americas <tibble [300 x 5]>     All  reg1 <S3: glm> <dbl [4]>
##  4  Americas <tibble [300 x 5]>     All  reg2 <S3: glm> <dbl [8]>
##  5      Asia <tibble [396 x 5]>     All  reg1 <S3: glm> <dbl [4]>
##  6      Asia <tibble [396 x 5]>     All  reg2 <S3: glm> <dbl [8]>
##  7    Europe <tibble [360 x 5]>     All  reg1 <S3: glm> <dbl [4]>
##  8    Europe <tibble [360 x 5]>     All  reg2 <S3: glm> <dbl [8]>
##  9   Oceania  <tibble [24 x 5]>     All  reg1 <S3: glm> <dbl [4]>
## 10   Oceania  <tibble [24 x 5]>     All  reg2 <S3: glm> <dbl [8]>

by_continent_with_reg %>% slice(1:8) %$% 
    stargazer(reg, se = se, type = "html",
        covariate.labels = c("Year", "Population", "GdpPerCap",
            "Year * Pop", "Year * GdpPerCap", "Pop * GdpPerCap",
            "Year * Pop * GdpPerCap"),
        omit.stat = c("ll", "aic"), omit.table.layout = "bl",
        add.lines = list(c("Continent", continent), c("Model", model)))

library(desctable)

stats_list = list(
  "Mean" = is.factor ~ percent | (is.double ~ mean),
  "SD"  = is.factor ~ NA | (is.double ~ sd)
)

gapminder %>%
    select(-c(country, year)) %>%
    filter(continent %in% c("Asia", "Europe")) %>%
    group_by(continent) %>%
    desctable(stats = stats_list) %>%
    datatable()