leeper / reggie Goto Github PK

Stata-like Regression Functionality for R

License: Other

Makefile 5.94% R 94.06%

reggie's Introduction

Stata-like Regression Functionality

This is a work-in-progress to explore how to design Stata-like regression modelling tools for R, namely those that allow plug-and-play variance-covariance estimation procedures and also to provide arguments to modelling functions in data-formula order (rather than the traditional formula-data order) thus enabling easy use in data analysis pipelines via %>%.

Contributions and feedback are welcome on GitHub.

Code Examples

In addition to plug-and-play variance-covariance procedures, the reg() function also provides pretty print methods.

library("reggie")

# reg
reg(ChickWeight, weight ~ Time + Diet)

## Generalized Linear Model
- Model:  weight ~ Time + Diet
- Family: gaussian (link: identity)
- Data (n=578): ChickWeight

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)    10.92       3.36     3.3    0.001
Time            8.75       0.22    39.5   <2e-16
Diet2          16.17       4.09     4.0    8e-05
Diet3          36.50       4.09     8.9   <2e-16
Diet4          30.23       4.11     7.4    2e-13

# reg
reg(ChickWeight, weight ~ Time + Diet, vcov_type = "const")

## Generalized Linear Model
- Model:  weight ~ Time + Diet
- Family: gaussian (link: identity)
- Data (n=578): ChickWeight

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)    10.92       3.36     3.3    0.001
Time            8.75       0.22    39.5   <2e-16
Diet2          16.17       4.09     4.0    8e-05
Diet3          36.50       4.09     8.9   <2e-16
Diet4          30.23       4.11     7.4    2e-13

# reg, vce(robust)
reg(ChickWeight, weight ~ Time + Diet, vcov_type = "HC0")

## Generalized Linear Model
- Model:  weight ~ Time + Diet
- Family: gaussian (link: identity)
- Data (n=578): ChickWeight

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)    10.92       2.82     3.9    1e-04
Time            8.75       0.26    33.6   <2e-16
Diet2          16.17       4.41     3.7    2e-04
Diet3          36.50       4.49     8.1    4e-16
Diet4          30.23       3.13     9.7   <2e-16

# reg, vce(boot)
reg(ChickWeight, weight ~ Time + Diet, vcov_type = "boot")

## Generalized Linear Model
- Model:  weight ~ Time + Diet
- Family: gaussian (link: identity)
- Data (n=578): ChickWeight

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)    10.92       2.80     3.9    1e-04
Time            8.75       0.26    33.3   <2e-16
Diet2          16.17       4.54     3.6    4e-04
Diet3          36.50       4.46     8.2    3e-16
Diet4          30.23       3.13     9.6   <2e-16

# reg, vce(cluster Chick)
reg(ChickWeight, weight ~ Time + Diet, vcov_cluster = ~Chick)

## Generalized Linear Model
- Model:  weight ~ Time + Diet
- Family: gaussian (link: identity)
- Data (n=578): ChickWeight

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)    10.92       5.39     2.0     0.04
Time            8.75       0.53    16.7   <2e-16
Diet2          16.17      10.91     1.5     0.14
Diet3          36.50       9.86     3.7    2e-04
Diet4          30.23       6.67     4.5    6e-06

# bootstrap, cluster(Chick) reps(5000): reg

# reg(ChickWeight, weight ~ Time + Diet, vcov_cluster = ~ Chick, vcov_type = 'boot')

# DOESN'T CURRENTLY WORK, BUT WHY?

# svy: reg
library("survey")
data(api)
dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, fpc = ~fpc)
reg(dstrat, api00 ~ ell + meals + mobility)

## Generalized Linear Model
- Model:  api00 ~ ell + meals + mobility
- Family: gaussian (link: identity)
- Data (n=200):
Stratified Independent Sampling design
svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)   820.89      10.08    81.5   <2e-16
ell            -0.48       0.39    -1.2      0.2
meals          -3.14       0.28   -11.1   <2e-16
mobility        0.23       0.39     0.6      0.6

The "model" object class contains the underlying model object as its model argument, and methods for various commonly used generic functions (coef(), vcov(), plot(), terms(), predict()) are provided that behave like those operations on a standard modelling object.

Installation

This package is not yet on CRAN. To install the latest development version you can pull a potentially unstable version directly from GitHub:

if (!require("remotes")) {
    install.packages("remotes")
}
remotes::install_github("leeper/reggie")

reggie's People

Contributors

Stargazers

Watchers

Forkers

b-rodrigues guhjy

reggie's Issues

Reconsider Data / Formula order in call

I get the point of the experiment with argument order, but I'm not sure %>% pipelines will be the dominant use case. And if they are not, then it's probably not a good idea to mess with a strong convention that is maintained by most model-fitting R packages.

Plus, the benefit of breaking convention doesn't seem very high, given that one can just do this:

dat %>% lm(y ~ x, .)

Model support

Priority areas:

{plm}
Stata's broader xt* family

But, need to decide whether to make these separate functions (named as they are in Stata) or pass a FUN argument specifying which estimation function to call.

Add test suite

Proof of concept: reggie and texreg

This morning, another student came to my office asking how to produce tables with clustered standard errors. So here I am again.

I just wanted to point out a neat feature that I contributed to texreg a little while ago: When the package doesn't recognize the model type, it tries to extract info using broom. Thanks to that, we can easily present side-by-side estimates from reggie by creating a couple very simple broom helpers:

library(reggie)
library(texreg)
library(broom)

# Broom Helper Functions
tidy.reg <- function(x, ...) {
    ret <- data.frame('term' = row.names(x$coefficients),
                      'estimate' = x$coefficients[, 1], 
                      'std.error' = x$coefficients[, 2], 
                      'statistic' = x$coefficients[, 3], 
                      'p.value' = x$coefficients[, 4], 
                      row.names = NULL,
  stringsAsFactors = FALSE)
    return(ret)
}
glance.reg <- function(x, ...) {
    ret <- glance(x$model)
    return(ret)
}

# Simulated Data
dat <- data.frame('x' = rnorm(1000),
                  'z' = sample(letters, 1000, replace = TRUE))
dat$y <- ifelse(dat$x + rnorm(1000) > 0, 1, 0)

# Fitting models
models <- list()
models[['GLM']] <- glm(y ~ x, data = dat, family = binomial())
models[['Cluster']] <- reg(y ~ x, data = dat, family = binomial(), vcov_cluster = ~ z)
models[['HC1']] <- reg(y ~ x, data = dat, family = binomial(), vcov_type = 'HC1')

# Summary
screenreg(models, digits = 4)

============================================================
                 GLM            Cluster        HC1 
------------------------------------------------------------
(Intercept)        -0.1355        -0.1355        -0.1355
                   (0.0774)       (0.0762)       (0.0774)
x                   1.6487 ***     1.6487 ***     1.6487 *** 
                   (0.1099)       (0.1060)       (0.1138)
------------------------------------------------------------
AIC              1013.3481      1013.3481      1013.3481
BIC              1023.1637      1023.1637      1023.1637
Log Likelihood   -504.6741      -504.6741      -504.6741
Deviance         1009.3481      1009.3481      1009.3481
Num. obs.        1000
Deviance (Null)                 1382.6922      1382.6922
df.null                          999            999 
DF Resid.                        998            998 
============================================================
*** p < 0.001, ** p < 0.01, * p < 0.05

Not perfect, but pretty close...

Error in reg(ChickWeight, weight ~ Time + Diet, vcov_cluster = ~ Chick, vcov_type = "boot")

As you mentioned in the README, this regression does not work and this is the error message:

object of type 'closure' is not subsettable

I believe it's probably from line 52:

data = data[sample(seq_len(nrow(data)), nrow(data), TRUE), ], ...))

data is a function to load example datasets, and I think that this is where it goes wrong (you're trying to subset a function). I tried replacing data with dataset and when running the unit tests this is what I get:

==> devtools::test()

Loading reggie
Loading required package: testthat
Testing reggie
✔ | OK F W S | Context
✖ |  0 1     | Correct Data Structures Returned [3.5 s]
─────────────────────────────────────────────────────────────────────────────────────────────────────────
tests.R:19: error: (unknown)
object 'dataset' not found
1: reg(ChickWeight, weight ~ Time + Diet, vcov_cluster = ~Chick, vcov_type = "boot") at /home/bro/Documents/reggie/tests/testthat/tests.R:19
2: sandwich::vcovBS(x = mod, cluster = cluster_vec, R = boot_iterations) at /home/bro/Documents/reggie/R/reg.R:43
3: vcovBS.default(x = mod, cluster = cluster_vec, R = boot_iterations)
4: eval(up, envir = env, enclos = parent.frame())
5: eval(up, envir = env, enclos = parent.frame())
6: stats::glm(formula = formula, data = dataset, subset = .vcovBSenv$.vcovBSsubset)
7: eval(mf, parent.frame())
8: eval(mf, parent.frame())
9: stats::model.frame(formula = formula, data = dataset, subset = .vcovBSenv$.vcovBSsubset, drop.unused.levels = TRUE)
10: model.frame.default(formula = formula, data = dataset, subset = .vcovBSenv$.vcovBSsubset, drop.unused.levels = TRUE)
11: is.data.frame(data)

I don't know if this helps, but thought I'd mention.

Arguments weights, subset, na.action, etc. are not passed successfully to ...

reg(mtcars, mpg ~ cyl, weights = runif(nrow(mtcars), 0, 1))

fails with the crazy ..1 error. Will probably need to do some match.call() magic.

leeper / reggie Goto Github PK

reggie's Introduction

Stata-like Regression Functionality

Code Examples

Installation

reggie's People

Contributors

Stargazers

Watchers

Forkers

reggie's Issues

Reconsider Data / Formula order in call

Model support

Add test suite

Proof of concept: reggie and texreg

Error in reg(ChickWeight, weight ~ Time + Diet, vcov_cluster = ~ Chick, vcov_type = "boot")

Arguments weights, subset, na.action, etc. are not passed successfully to ...

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent