Merge pull request #245 from hughjonesd/ivreg-methods

Ivreg methods

DESCRIPTION

@@ -86,7 +86,8 @@ Suggests:
     lsmeans,
     betareg,
     robust,
-    akima
+    akima,
+    AER
 URL: http://github.com/tidyverse/broom
 BugReports: http://github.com/tidyverse/broom/issues
 VignetteBuilder: knitr

NAMESPACE

@@ -8,6 +8,7 @@ S3method(augment,data.frame)
 S3method(augment,default)
 S3method(augment,felm)
 S3method(augment,glmRob)
+S3method(augment,ivreg)
 S3method(augment,kmeans)
 S3method(augment,lm)
 S3method(augment,lmRob)
@@ -48,6 +49,7 @@ S3method(glance,glmRob)
 S3method(glance,glmnet)
 S3method(glance,gmm)
 S3method(glance,htest)
+S3method(glance,ivreg)
 S3method(glance,kmeans)
 S3method(glance,list)
 S3method(glance,lm)
@@ -126,6 +128,7 @@ S3method(tidy,glmRob)
 S3method(tidy,glmnet)
 S3method(tidy,gmm)
 S3method(tidy,htest)
+S3method(tidy,ivreg)
 S3method(tidy,kappa)
 S3method(tidy,kde)
 S3method(tidy,kmeans)

R/ivreg_tidiers.R

@@ -0,0 +1,108 @@
+#' Tidiers for ivreg models
+#' 
+#' @param x an "ivreg" object 
+#' @param data original dataset
+#' @param conf.int whether to include a confidence interval
+#' @param conf.level confidence level of the interval, used only if
+#' \code{conf.int=TRUE}
+#' @param exponentiate whether to exponentiate the coefficient estimates
+#' and confidence intervals
+#' 
+#' @template boilerplate
+#' 
+#' @return \code{tidy.ivreg} returns a data frame with one row per
+#' coefficient, of the same form as \code{\link{tidy.lm}}.
+#' 
+#' @seealso \code{\link{lm_tidiers}}
+#' 
+#' @name ivreg_tidiers
+#' 
+#' @examples
+#' 
+#' if (require("AER", quietly = TRUE)) {
+#'     data("CigarettesSW", package = "AER")
+#'     CigarettesSW$rprice <- with(CigarettesSW, price/cpi)
+#'     CigarettesSW$rincome <- with(CigarettesSW, income/population/cpi)
+#'     CigarettesSW$tdiff <- with(CigarettesSW, (taxs - tax)/cpi)
+#'     ivr <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi),
+#'           data = CigarettesSW, subset = year == "1995")
+#'     
+#'     summary(ivr)
+#'     
+#'     tidy(ivr)
+#'     tidy(ivr, conf.int = TRUE)
+#'     tidy(ivr, conf.int = TRUE, exponentiate = TRUE)
+#'     
+#'     head(augment(ivr))
+#'     
+#'     glance(ivr)
+#' }
+#' 
+#' @export
+tidy.ivreg <- function(x, conf.int = FALSE, conf.level = .95,
+    exponentiate = FALSE, ...) {
+
+    co <- stats::coef(summary(x))
+    nn <- c("estimate", "std.error", "statistic", "p.value")
+    ret <- fix_data_frame(co, nn[1:ncol(co)])
+
+   process_lm(ret, x, conf.int = conf.int, conf.level = conf.level, exponentiate = exponentiate)
+}
+
+
+#' @rdname ivreg_tidiers
+#' 
+#' @return \code{augment} returns a data frame with one row for each
+#' initial observation, adding the columns:
+#'   \item{.fitted}{predicted (fitted) values}
+#' and if \code{newdata} is \code{NULL}:
+#'   \item{.resid}{residuals}
+#'  
+#' 
+#' @export
+augment.ivreg <- function(x, data = as.data.frame(stats::model.frame(x)), newdata, ...) {
+    augment_columns(x, data, newdata, ...)
+}
+
+
+#' @rdname ivreg_tidiers
+#' 
+#' @param ... extra arguments, not used
+#' @param diagnostics Logical. Return results of diagnostic tests.
+#' 
+#' @return \code{glance} returns a one-row data frame with columns
+#'   \item{r.squared}{The percent of variance explained by the model}
+#'   \item{adj.r.squared}{r.squared adjusted based on the degrees of freedom}
+#'   \item{statistic}{Wald test statistic}
+#'   \item{p.value}{p-value from the Wald test}
+#'   \item{df}{Degrees of freedom used by the coefficients}
+#'   \item{sigma}{The square root of the estimated residual variance}
+#'   \item{df.residual}{residual degrees of freedom}
+#' If \code{diagnostics} is \code{TRUE}, \code{glance} also returns:
+#'   \item{p.value.Sargan}{P value of Sargan test}
+#'   \item{p.value.Wu.Hausman}{P value of Wu-Hausman test}
+#'   \item{p.value.weakinst}{P value of weak instruments test}
+#' 
+#' @export
+glance.ivreg <- function(x, diagnostics = FALSE, ...) {
+    s <- summary(x, diagnostics = diagnostics)
+    ret <- with(s, data.frame(
+        r.squared     = r.squared,
+        adj.r.squared = adj.r.squared,
+        sigma         = sigma,
+        statistic     = waldtest[1],
+        p.value       = waldtest[2],
+        df            = df[1]
+    ))
+    if (diagnostics) {
+        ret <- cbind(ret, with(s, data.frame(
+            statistic.Sargan     = diagnostics["Sargan", "statistic"],
+            p.value.Sargan       = diagnostics["Sargan", "p-value"],
+            statistic.Wu.Hausman = diagnostics["Wu-Hausman", "statistic"],
+            p.value.Wu.Hausman   = diagnostics["Wu-Hausman", "p-value"],
+            statistic.weakinst   = diagnostics["Weak instruments", "statistic"],
+            p.value.weakinst     = diagnostics["Weak instruments", "p-value"]    
+        )))
+    }
+    finish_glance(ret, x)
+}

tests/testthat/test-ivreg.R

@@ -0,0 +1,34 @@
+# test tidy, augment, glance methods from lme4-tidiers.R
+
+if (require(AER, quietly = TRUE)) {
+    context("AER::ivreg models")
+
+    data("CigarettesSW", package = "AER")
+    CigarettesSW$rprice <- with(CigarettesSW, price/cpi)
+    CigarettesSW$rincome <- with(CigarettesSW, income/population/cpi)
+    CigarettesSW$tdiff <- with(CigarettesSW, (taxs - tax)/cpi)
+    ivr <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi),
+          data = CigarettesSW, subset = year == "1995")
+
+    test_that("tidy works on AER::ivreg fits", {
+        td <- tidy(ivr)
+        td2 <- tidy(ivr, conf.int = TRUE)
+        expect_warning(tidy(ivr, exponentiate = TRUE)) # warning as we didn't use a link function, maybe this is bad?
+    })
+
+    test_that("augment works on ivreg fits", {
+        au <- augment(ivr)
+        expect_true(all(c('.resid', '.fitted') %in% names(au)))
+        expect_equivalent(au$.resid, residuals(ivr))
+        expect_equivalent(au$.fitted, fitted(ivr))
+        old_cigs <- CigarettesSW[CigarettesSW$year == "1985" & CigarettesSW$tax < 40, ]
+        au2 <- augment(ivr, newdata = old_cigs)
+        expect_true('.fitted' %in% names(au2))
+        expect_equivalent(au2$.fitted, predict(ivr, newdata = old_cigs))
+    })
+
+    test_that("glance works on ivreg fits", {
+        g <- glance(ivr)
+        g <- glance(ivr, diagnostics = TRUE)
+    })
+}