README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r setup, include = FALSE}
library(knitr)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(gravity)
library(eflm)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# Efficient Fitting of Linear Models

<!-- badges: start -->
[![Project Status: Active – The project has reached a stable, usable
state and is being actively
developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
[![Lifecycle: stable](https://img.shields.io/badge/lifecycle-stable-brightgreen.svg)](https://www.tidyverse.org/lifecycle/#stable)
[![CRAN
status](https://www.r-pkg.org/badges/version/gravity)](https://cran.r-project.org/package=eflm)
[![codecov](https://codecov.io/gh/pachadotdev/eflm/branch/main/graph/badge.svg?token=XI59cmGd15)](https://codecov.io/gh/pachadotdev/eflm)
[![R-CMD-check](https://github.com/pachadotdev/eflm/workflows/R-CMD-check/badge.svg)](https://github.com/pachadotdev/eflm/actions)
<!-- badges: end -->

## Scope

`eflm` package reduces the design matrix from $N\times P$ into 
$P \times P$ for reduced fitting time, and delivers functions that are drop-in 
replacements for `glm` and `lm`, like:
```{r, eval = FALSE}
# just append and 'e' to glm
eglm(mpg ~ wt, data = mtcars)
```

The best computational performance is obtained when R is linked against OpenBLAS,
Intel MKL or other optimized BLAS library. This implementation aims at being
compatible with 'broom' and 'sandwich' packages for summary statistics and
clustering by providing S3 methods.

This package takes ideas from glm2, speedglm, fastglm, speedglm and fixest 
packages, but the implementations here shall keep the functions and outputs as 
closely as possible to the stats package, therefore making the functions 
provided here compatible with packages such as sandwich for robust estimation, 
even if that means to attenuate the speed gains.

The greatest strength of this package is testing. With more than 1600 
(and counting) tests, we try to do exactly the same as lm/glm, even in edge 
cases, but faster.

The ultimate aim of the project is to produce a package that:

* Does exactly the same as lm and glm in less time
* Is equally numerically stable as lm and glm
* Depends only on base R, with no Rcpp or other calls
* Uses R's internal C code such as the `Cdqrls` function that the stats package uses for model fitting
* Can be used in Shiny dashboard and contexts where you need fast model fitting
* Is useful for memory consuming models
* Allows model fitting in cases demanding more memory than free RAM (PENDING)

## Installation

You can install the released version of eflm from CRAN with:
```r
install.packages("eflm")
```

And the development version with:
``` r
remotes::install_github("pachadotdev/eflm")
```

## Progress list

### Stats compatibility

- [x] cooks.distance

### Sandwich compatibility

- [x] estfun
- [x] bread
- [x] vcovCL
- [x] meatCL
- [x] vcovCL
- [x] vcovBS
- [ ] vcovHC
- [ ] meatHC
- [ ] vcovPC
- [ ] meatPC
- [ ] vcovPL
- [ ] meatPL

### Broom compatibility

- [x] augment
- [x] tidy
- [x] glance

### Lmtest compatibility

- [x] resettest

## Benchmarking

The dataset for this benchmark was taken from Yotov et al. (2016) and consists 
in a 28,152 x 8 data frame with 6 numeric and 2 categorical columns of the form:

|Year ($t$) |Trade ($X$) |DIST |CNTG |LANG |CLNY |Exp Year ($\pi$) |Imp Year ($\chi$) |
|-----------|------------|-----|-----|-----|-----|-----------------|------------------|
|1986       |27.8        |12045|0    |0    |0    |ARG1986          |AUS1986           |
|1986       |3.56        |11751|0    |0    |0    |ARG1986          |AUT1986           |
|1986       |96.1        |11305|0    |0    |0    |ARG1986          |BEL1986           |

This data can be found in the `tradepolicy` package.

The variables are:

* `year`: time of export/import flow
* `trade`: bilateral trade
* `log_dist`: log of distance
* `cntg`: contiguity (0/1)
* `lang`: common language (0/1)
* `clny`: colonial relation (0/1)
* `exp_year`/`imp_year`: exporter/importer time fixed effects

For benchmarking I'll fit a PPML model, as it's a computationally expensive model.

```r
ch1_application1 <- tradepolicy::agtpa_applications %>%
  select(exporter, importer, pair_id, year, trade, dist, cntg, lang, clny) %>%
  filter(year %in% seq(1986, 2006, 4))
  
formula <- trade ~ log(dist) + cntg + lang + clny + exp_year + imp_year
eglm(formula, quasipoisson, ch1_application1)
```

To compare `glm`, the proposed `eglm` and Stata's `ppml`, I conducted a test 
with 500 repetitions locally, and reported the median of the realizations as
the fitting time. The plots on the right report the fitting times and used 
memory by running regressions with cumulative subset of the data for 
1986, ..., 2006 (e.g. regress for 1986, then 1986 and 1990, ..., 
then 1986 to 2006), we obtain the next fitting times and memory allocation 
depending on the design matrix dimensions:

```{r echo=FALSE, fig.height=3, fig.width=6, message=FALSE, warning=FALSE, dpi=150, fig.align="center"}
g <- readRDS("benchmarks/06-stata-vs-r.rds")
g2 <- readRDS("benchmarks/05-benchmark-subsets.rds")
g + (g2[[1]] / g2[[2]])
```

Yotov et al. (2016) features complex both partial and general equilibrium 
models. Some partial equilibrium models are particularly slow to fit because of
the allocated memory and the number of fixed effects, such as the Regional Trade 
Agreements (RTAs) model.

In the next table, TG means 'Traditional Gravity' (e.g. vanilla PPML), 
DP means 'Distance Puzzle' and GB stands for 'Globalization', which are 
refinements of the simple PPML model and include dummy variables such as 
specific country pair fixed effects and lagged RTAs.

```{r echo=FALSE}
load("benchmarks/04-glm-vs-eglm.RData")
t1 <- bind_rows(
  model_matrix_traditional_gravity %>% 
    filter(fun %in% c("eglm"), model_matrix_ncol == max(model_matrix_ncol)),
  model_matrix_distance_puzzle %>% 
    filter(fun %in% c("eglm"), model_matrix_ncol == max(model_matrix_ncol)),
  model_matrix_rtas %>% 
    filter(fun %in% c("eglm"), model_matrix_ncol == max(model_matrix_ncol))
) %>% 
  select(-fun)
colnames(t1) <- c("Model", "Rows in design matrix", "Cols in design matrix")
kable(t1)
```

```{r echo=FALSE, fig.height=2, fig.width=6, message=FALSE, warning=FALSE, dpi=150, fig.align="center"}
include <- t1$Model

g2_time <- ggplot(bench_models %>% filter(fun %in% c("eglm", "glm"), model %in% include)) + 
  geom_col(aes(x = model, y = median_seconds, fill = `Function`), position = "dodge2") + 
  scale_fill_manual(values = c("#3B9AB2", "#E1AF00")) + 
  theme_minimal(base_size = 8) + 
  # theme(legend.position = "none") +
  labs(x = "Model", y = "Median Time (Seconds)",
       title = "Time Benchmark")
g2_memory <- ggplot(bench_models %>% filter(fun %in% c("eglm", "glm"), model %in% include)) + 
  geom_col(aes(x = model, y = mem_alloc_mb, fill = `Function`), position = "dodge2") + 
  scale_fill_manual(values = c("#3B9AB2", "#E1AF00")) + 
  theme_minimal(base_size = 8) + 
  # theme(legend.position = "none") +
  labs(x = "Model", y = "Required Memory (MB)",
       title = "Memory Benchmark")
g2_time + g2_memory
```

The results for the RTA model show that the speedups can be scaled, and we can 
show both time reduction and required memory increases.

```{r echo=FALSE}
t2 <- bench_models %>% 
  filter(fun %in% c("eglm", "glm"), model %in% include) %>% 
  select(model, fun, median_seconds) %>% 
  spread(fun, median_seconds) %>% 
  select(model, glm, eglm) %>% 
  mutate(time_gain = paste0(round(100 * (glm - eglm) / glm, 2), "%"))
colnames(t2) <- c("Model", "GLM Time (s)", "EGLM Time (s)", "Time Gain (%)")
kable(t2)
```

Is it important to mention that the increase in memory results in reduced object
size for the stored model.

```{r echo=FALSE}
t4 <- bind_rows(
  object_size_traditional_gravity %>% 
    filter(fun %in% c("glm", "eglm"), model %in% include),
  object_size_distance_puzzle %>% 
    filter(fun %in% c("glm", "eglm"), model %in% include),
  object_size_rtas %>% 
    filter(fun %in% c("glm", "eglm"), model %in% include)
) %>% 
  mutate(fun = ifelse(fun == "glm", "GLM", "EGLM"))
t4 <- t4 %>% 
  select(model, fun, size) %>% 
  mutate(size = round(size, 2)) %>% 
  spread(fun, size) %>% 
  select(model, GLM, EGLM) %>% 
  mutate(object_gain = paste0(round(100 * (GLM - EGLM) / GLM, 2), "%"))
colnames(t4) <- c("Model", "GLM Size (MB)", "EGLM Size (MB)", "Memory Savings (%)")
kable(t4)
```

To conclude my benchmarks, I fitted the PPML model again on DigitalOcean 
droplets, leading to consistent times across scaled hardware. The results can be
seen in the next plot:
```{r echo=FALSE, fig.height=2, fig.width=6, message=FALSE, warning=FALSE, dpi=150, fig.align="center"}
g_do <- readRDS("benchmarks/07-r-digitalocean.rds")
g_do
```

## Edge cases

An elementary example that breaks `eflm` even with QR decomposition can be 
found in Golub et al. (2013), which consists in passing an ill conditioned 
matrix:
```{r include=FALSE}
set.seed(200100); n <- 10000
df <- data.frame(x1 = 1:n, x2 = c(0, 2:n))
df$y <- with(df, 1.999 + 2.000 * x1 + 2.001 * x2 + rnorm(n))
reg1 <- lm(y ~ ., df); reg2 <- elm(y ~ ., df)
```

| Model | (Intercept) | $x_1$ | $x_2$ | 
|-------|-------------|-------|-------|
|REG 1  | 1.98        | 2.98  | 1.02  |
|REG 2  | 1.98        | 4.00  | NA    |

# References

Golub, Gene H, and Charles F Van Loan. 2013. *Matrix Computations*. Vol. 3. JHU press.

Yotov, Yoto V, Roberta Piermartini, José-Antonio Monteiro, and Mario Larch. 2016. 
*An Advanced Guide to Trade Policy Analysis: The Structural Gravity Model*. World Trade
Organization Geneva.