Add p-values

DESCRIPTION

@@ -28,7 +28,7 @@ Suggests:
     knitr,
     rmarkdown,
     formatR
-RoxygenNote: 7.2.0
+RoxygenNote: 7.3.1
 URL: https://github.com/kolesarm/RDHonest
 VignetteBuilder: knitr
 Language: en-US

R/NPR_lp.R

@@ -5,9 +5,6 @@
 ## @param T0bias When evaluating the maximum bias of the estimate, use the
 ##     estimate itself (if \code{T0bias==FALSE}), or use the preliminary
 ##     estimate \code{T0} (if \code{T0bias==TRUE}). Only relevant for Fuzzy RD.
-## @return Returns an object of class \code{"NPRResults"}, see descriptions in
-##     \code{\link{RDHonest}}, \code{\link{LPPHonest}}, and
-##     \code{\link{FRDHonest}}.
 NPRHonest.fit <- function(d, M, kern="triangular", h, opt.criterion, alpha=0.05,
                           beta=0.8, se.method="nn", J=3, sclass="H", T0=0,
                           T0bias=FALSE) {
@@ -54,8 +51,8 @@ NPRHonest.fit <- function(d, M, kern="triangular", h, opt.criterion, alpha=0.05,
     }
     lower <- r1$estimate - bias - stats::qnorm(1-alpha)*r1$se
     upper <- r1$estimate + bias + stats::qnorm(1-alpha)*r1$se
-    cv <- CVb(bias/r1$se, alpha)
-
+    B <- bias/r1$se
+    cv <- CVb(B, alpha)
     term <- switch(d$class, IP="Value of conditional mean",
                    SRD="Sharp RD Parameter", "Fuzzy RD Parameter")
     method <- switch(sclass, H="Holder", "Tayor")
@@ -69,7 +66,9 @@ NPRHonest.fit <- function(d, M, kern="triangular", h, opt.criterion, alpha=0.05,
                        conf.high.onesided=upper, bandwidth=h,
                        eff.obs=r1$eff.obs, leverage=max(r1$w^2)/sum(r1$w^2),
                        cv=cv, alpha=alpha, method=method, M=M[1], M.rf=M[2],
-                       M.fs=M[3], first.stage=r1$fs, kernel=kernel)
+                       M.fs=M[3], first.stage=r1$fs, kernel=kernel,
+                       p.value=stats::pnorm(B-abs(r1$estimate/r1$se))+
+                           stats::pnorm(-B-abs(r1$estimate/r1$se)))
     structure(list(coefficients=coef, data=d), class="RDResults")
 }
 
@@ -149,6 +148,7 @@ print.RDResults <- function(x, digits = getOption("digits"), ...) {
     else
         cat("\nSmoothness constant M:", fmt(y$M),
             "\n")
+    cat("P-value:", fmt(y$p.value), "\n")
 
     if (inherits(x$na.action, "omit"))
         cat(length(x$na.action), "observations with missing values dropped\n")

R/RDHonest.R

@@ -72,9 +72,9 @@
 #'     optimal bandwidth. Only relevant for Fuzzy RD.
 #' @param sigma2 Supply variance. Ignored when kernel is optimal.
 #' @param clusterid Cluster id for cluster-robust standard errors
-#' @return Returns an object of class \code{"NPRResults"}. The function
+#' @return Returns an object of class \code{"RDResults"}. The function
 #'     \code{print} can be used to obtain and print a summary of the results. An
-#'     object of class \code{"NPRResults"} is a list containing the following
+#'     object of class \code{"RDResults"} is a list containing the following
 #'     components
 #'
 #'     \describe{
@@ -152,6 +152,8 @@ RDHonest <- function(formula, data, subset, weights, cutoff=0, M,
     stopifnot(length(formula)[1] == 1L, length(formula)[2] <= 2)
     mf$formula <- formula
 
+    ## http://madrury.github.io/jekyll/update/statistics/2016/07/20/lm-in-R.html
+
     mf[[1L]] <- quote(stats::model.frame)
     mf <- eval(mf, parent.frame())
 

R/plots.R

@@ -57,11 +57,11 @@ plot_RDscatter <- function(formula, data, subset, cutoff=0, na.action, avg=10,
         bd$count <- avg
     }
     bd$x <- bd$x+d$orig.cutoff
-
     if (propdotsize) {
-        p <- ggplot2::qplot(x=bd$x, y=bd$y, size=bd$count)
+        p <- ggplot2::ggplot()+
+            ggplot2::geom_point(ggplot2::aes(x=bd$x, y=bd$y, size=bd$count))
     } else {
-        p <- ggplot2::qplot(x=bd$x, y=bd$y)
+        p <- ggplot2::ggplot()+ggplot2::geom_point(ggplot2::aes(x=bd$x, y=bd$y))
     }
 
     p <- p + ggplot2::theme(legend.position = "none")

README.md

@@ -13,9 +13,6 @@ Kolesár (2020)](https://doi.org/10.3982/QE1199)
 [RDHonest-vStata](https://github.com/tbarmstr/RDHonest-vStata) for a Stata
 version of this package.
 
-
-
-
 See vignette [RDHonest](doc/RDHonest.pdf) for description of the package
 (available through `vignette("RDHonest")` once package is installed), and the
 package [manual](doc/manual.pdf) for documentation of the package functions.

doc/RDHonest.R

@@ -1,4 +1,4 @@
-## ---- include=FALSE, cache=FALSE----------------------------------------------
+## ----include=FALSE, cache=FALSE-----------------------------------------------
 library("knitr")
 knitr::opts_knit$set(self.contained = FALSE)
 knitr::opts_chunk$set(tidy = TRUE, collapse=TRUE, comment = "#>",
@@ -7,15 +7,15 @@ knitr::opts_chunk$set(tidy = TRUE, collapse=TRUE, comment = "#>",
 ## -----------------------------------------------------------------------------
 library("RDHonest")
 
-## ---- fig.width=4.5, fig.height=3.5, fig.cap="Lee (2008) data"----------------
+## ----fig.width=4.5, fig.height=3.5, fig.cap="Lee (2008) data"-----------------
 ## plot 25-bin averages in for observations 50 at most points
 ## away from the cutoff. See Figure 1.
 plot_RDscatter(voteshare~margin, data=lee08, subset=abs(lee08$margin)<=50,
             avg=50, propdotsize=FALSE,
             xlab="Margin of victory",
             ylab="Vote share in next election")
 
-## ---- fig.width=4.5, fig.height=3.5, fig.cap="Oreopoulos (2006) data"---------
+## ----fig.width=4.5, fig.height=3.5, fig.cap="Oreopoulos (2006) data"----------
 ## see Figure 2
 f2 <- plot_RDscatter(I(log(earnings))~yearat14, data=cghs, cutoff=1947,
     avg=Inf, xlab="Year aged 14", ylab="Log earnings",

tests/testthat/test_clustering.R

@@ -55,16 +55,17 @@ test_that("Test clustering formulas", {
     ##                       cadjust=FALSE)[1, 1] 651509.586
     expect_equal(p1c$coefficients$std.error^2, 723178.2655)
     expect_equal(p1h$coefficients$std.error, 792.2299439)
+
     p2c <- RDHonest(c~elig_year, data=rr, clusterid=clusterid,
                     point.inference=TRUE)
-    h <- p2c$coefficients$bandwidth
-    m3 <- lm(c~elig_year, data=rr, subset=abs(elig_year)<=h, weights
-             =(1 - abs(elig_year/h)))
+    h <- p2c$coefficients$bandwidth #8.27778
+    m3 <- lm(c~elig_year, data=rr, subset=abs(elig_year)<=h,
+             weights =(1 - abs(elig_year/h)))
     ## sqrt(sandwich::vcovCL(m3, type="HC0",
     ##                       cluster=clusterid[abs(rr$elig_year)<=h],
     ##                       cadjust=FALSE)[1, 1])
     expect_equal(as.numeric(p2c$coefficients[2:3]),
-                 c(unname(m3$coefficients[1]), 919.1556447))
+                 c(unname(m3$coefficients[1]), 919.162944))
 
     f1c <- RDHonest(c~retired | elig_year, data=rr, clusterid=clusterid,
                     kern="uniform")

tests/testthat/test_frd.R

@@ -35,6 +35,10 @@ test_that("FRD data example check", {
                         T0=0)$coefficients
     r2 <- NPRHonest.fit(d, M, kern="triangular", opt.criterion="MSE",
                         T0=r1$estimate)$coefficients
+    r1a <- NPRHonest.fit(d, M, kern="triangular", opt.criterion="MSE",
+                        T0=0, alpha=r1$p.value)$coefficients
+    expect_equal(r1a$conf.high, 0)
+    expect_equal(max(abs(expect_equal(r1a$conf.high, 0))), 0)
     ## With positive T0, expect greater effective M, and thus smaller bandwidth
     expect_lt(r2$bandwidth, r1$bandwidth)
     ## On travis, only the first 6 digits match, not sure why

tests/testthat/test_rd.R

@@ -91,7 +91,7 @@ test_that("Honest inference in Lee and LM data",  {
     expect_equal(r1o[8], m1)
     expect_equal(r1o[10], "Onesided CIs:  (-Inf, 4.684), (-7.081, Inf)")
     expect_equal(r2o[10], "Onesided CIs:  (-Inf, 4.742), (-7.138, Inf)")
-    expect_equal(r1o[15], "24 observations with missing values dropped")
+    expect_equal(r1o[16], "24 observations with missing values dropped")
 
     d <- NPRData(headst[!is.na(headst$mortHS), c("mortHS", "povrate")],
                 cutoff=0, "SRD")
@@ -202,13 +202,13 @@ test_that("Honest inference in Lee and LM data",  {
     r1 <- RDHonest(mortHS ~ povrate, data=headst, kern="uniform",
                    na.action="na.omit")
     r1 <- capture.output(print(r1, digits=6))
-    expect_equal(r1[c(11, 12, 15)],
+    expect_equal(r1[c(11, 12, 16)],
                  c("Bandwidth: 3.98048, Kernel: uniform",
                    "Number of effective observations:     239",
                    "24 observations with missing values dropped"))
     r2 <- RDHonest(mortHS ~ povrate, data=headst, kern="epanechnikov",
                    point.inference=TRUE, cutoff=4)
-    expect_equal(capture.output(print(r2, digits=6))[c(11, 12, 15)],
+    expect_equal(capture.output(print(r2, digits=6))[c(11, 12, 16)],
                  c("Bandwidth: 9.42625, Kernel: epanechnikov",
                    "Number of effective observations: 373.642",
                    "24 observations with missing values dropped"))