Skip to content

Commit

Permalink
Update print method
Browse files Browse the repository at this point in the history
  • Loading branch information
@kolesarm
kolesarm committed Mar 12, 2024
1 parent 1f15c77 commit 661fe6c
Show file tree
Hide file tree
Showing 7 changed files with 76 additions and 56 deletions.
18 changes: 8 additions & 10 deletions R/NPRfunctions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
## Keep only positive kernel weights
X <- drop(d$X)
W <- if (h<=0) 0*X else kern(X/h)*d$w # kernel weights

## if (!is.null(d$sigma2)) d$sigma2 <- as.matrix(d$sigma2)
Z <- outer(X, 0:order, "^")
nZ <- c("(Intercept)", colnames(d$X), paste0("I(", colnames(d$X), "^",
Expand All @@ -24,18 +23,18 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
}
r0 <- stats::lm.wfit(x=Z, y=d$Y, w=W)
class(r0) <- c(if (ncol(d$Y)>1) "mlm", "lm")

if (any(is.na(r0$coefficients))) {
return(list(estimate=0, se=NA, est_w=W*0,
sigma2=NA*d$Y, eff.obs=0, fs=NA, lm=r0))
return(list(estimate=0, se=NA, est_w=W*0, sigma2=NA*d$Y, eff.obs=0,
fs=NA, lm=r0))
}
wgt <- W*0
ok <- W!=0
wgt[ok] <- solve(qr.R(r0$qr), t(sqrt(W[W>0])*qr.Q(r0$qr)))[1, ]
## To compute effective observations, rescale against uniform kernel
q_u <- qr(sqrt(d$w * (abs(X)<=h)) * Z)
wgt_u <- solve(qr.R(q_u), t(sqrt(d$w * (abs(X)<=h)) * qr.Q(q_u)))[1, ]
eff.obs <- sum((abs(X)<=h)*d$w)*sum(wgt_u^2/d$w)/sum(wgt^2/d$w)
Wu <- d$w * (abs(X)<=h)
q_u <- qr(sqrt(Wu) * Z)
wgt_u <- solve(qr.R(q_u), t(sqrt(Wu) * qr.Q(q_u)))[1, ]
eff.obs <- sum(Wu)*sum(wgt_u^2/d$w)/sum(wgt^2/d$w)

ny <- NCOL(r0$residuals)
## Squared residuals, allowing for multivariate Y
Expand All @@ -56,8 +55,8 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
res
}

hsigma2 <- switch(se.method, nn=NN(X),
EHW=HC(as.matrix(r0$residuals)), supplied.var=d$sigma2)
hsigma2 <- switch(se.method, nn=NN(X), EHW=HC(as.matrix(r0$residuals)),
supplied.var=as.matrix(d$sigma2))

## Variance
if (is.null(d$clusterid)) {
Expand All @@ -72,7 +71,6 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
}
ret <- list(estimate=r0$coefficients[1], se=sqrt(V[1]), est_w=wgt,
sigma2=hsigma2, eff.obs=eff.obs, fs=NA, lm=r0)
## gamma=beta[seq_len(NCOL(d$covs))+NROW(r0$coefficients)-NCOL(d$covs), ])

if (inherits(d, "FRD")) {
ret$fs <- r0$coefficients[1, 2]
Expand Down
38 changes: 22 additions & 16 deletions R/RDHonest.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,36 +347,42 @@ print.RDResults <- function(x, digits = getOption("digits"), ...) {
}
fmt <- function(x) format(x, digits=digits, width=digits+1)
y <- x$coefficients
cat("Estimates (using ", y$method, " class):\n", sep="")
nm <- c("Parameter", "Estimate", "Std. Error", "Maximum Bias")
names(y)[1:4] <- nm
cat("Inference for ", y$term, " (using ", y$method,
" class), confidence level ", 100*(1-y$alpha), "%:\n", sep="")
nm <- c("Estimate", "Std. Error", "Maximum Bias")
colnames(y)[2:4] <- nm
y$"Confidence Interval" <- paste0("(", fmt(y$conf.low), ", ",
fmt(y$conf.high), ")")
y$OCI <- paste0("(-Inf, ", fmt(y$conf.high.onesided), "), (",
fmt(y$conf.low.onesided), ", Inf)")

print.data.frame(y[, c(nm, "Confidence Interval"), ],
digits=digits, row.names=FALSE)
digits=digits)
cat("\nOnesided CIs: ", y$OCI)
if (!is.null(y$bandwidth))
cat("\nBandwidth: ", fmt(y$bandwidth), ", Kernel: ", y$kernel, sep="")
else
cat("\nSmoothing parameters below and above cutoff: ",
fmt(y$bandwidth.m), ", ", fmt(y$bandwidth.p), sep="")

cat("\nNumber of effective observations:", fmt(y$eff.obs))
par <- paste0(tolower(substr(y$Parameter, 1, 1)), substring(y$Parameter, 2))
par <- paste0(tolower(substr(y$term, 1, 1)), substring(y$term, 2))
cat("\nMaximal leverage for ", par, ": ", fmt(y$leverage),
sep="")
if (!is.null(y$first.stage) && !is.na(y$first.stage))
cat("\nFirst stage estimate:", fmt(y$first.stage),
"\nFirst stage smoothness constant M:", fmt(y$M.fs),
"\nReduced form smoothness constant M:", fmt(y$M.rf),
"\n")
else
cat("\nSmoothness constant M:", fmt(y$M),
"\n")
cat("P-value:", fmt(y$p.value), "\n")
"\nReduced form smoothness constant M:", fmt(y$M.rf))
else if (!is.null(y[["M"]])) {
cat("\nSmoothness constant M:", fmt(y$M))
}

cat("\nP-value:", fmt(y$p.value), "\n")

if (!is.null(y$bandwidth)) {
cat("\nBased on local regression with bandwidth: ", fmt(y$bandwidth),
", kernel: ", y$kernel, "\nRegression coefficients:\n", sep="")
print.default(format(coef(x$lm), digits = max(3L, digits - 3L)),
print.gap = 2L, quote = FALSE)
} else {
cat("\nSmoothing parameters below and above cutoff: ",
fmt(y$bandwidth.m), ", ", fmt(y$bandwidth.p), sep="")
}
if (inherits(x$na.action, "omit"))
cat(length(x$na.action), "observations with missing values dropped\n")

Expand Down
24 changes: 14 additions & 10 deletions R/RD_bme.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,18 +125,22 @@ RDHonestBME <- function(formula, data, subset, cutoff=0, na.action,
oci_l <- stats::coef(m1)[order+2]+dev-stats::qnorm(1-alpha)*se #
oci_u <- stats::coef(m1)[order+2]+dev+stats::qnorm(1-alpha)*se

wt <- solve(qr.R(m1$qr), t(qr.Q(m1$qr)))[order+2, ]

l <- which.min(ci_l)
u <- which.max(ci_u)
coef <- data.frame(term="Sharp RD parameter",
estimate=unname(m1$coefficients[order +2]),
std.error=sqrt(vdt["e2", "e2"]),
maximum.bias=max(abs(c(dev[u], dev[l]))),
conf.low=ci_l[l], conf.high=ci_u[u],
conf.low.onesided=min(oci_l),
conf.high.onesided=max(oci_u), bandwidth=h,
eff.obs=length(x), cv=NA, alpha=alpha, method="BME",
kernel="uniform")
structure(list(coefficients=coef, call=cl,
co <- data.frame(term="Sharp RD parameter",
estimate=m1$coefficients[order +2],
std.error=sqrt(vdt["e2", "e2"]),
maximum.bias=max(abs(c(dev[u], dev[l]))), conf.low=ci_l[l],
conf.high=ci_u[u], conf.low.onesided=min(oci_l),
conf.high.onesided=max(oci_u), bandwidth=h,
eff.obs=length(x), leverage=max(wt^2)/sum(wt^2),
cv=NA, alpha=alpha,
method="BME", kernel="uniform")
co$p.value <- stats::pnorm(co$maximum.bias-abs(co$estimate/co$std.error))+
stats::pnorm(-co$maximum.bias-abs(co$estimate/co$std.error))
structure(list(coefficients=co, call=cl, lm=m1,
na.action=attr(mf, "na.action")),
class="RDResults")
}
1 change: 1 addition & 0 deletions R/RD_opt.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@ RDTEstimator <- function(d, f, alpha, se.method, J) {
r.u <- NPReg(d, max(abs(d$X[W!=0])), kern="uniform")
eff.obs <- r.u$eff.obs*sum(r.u$w^2)/sum(W^2)
d$est_w <- W
names(Lhat) <- colnames(d$X)
co <- data.frame(term="Sharp RD parameter", estimate=Lhat, std.error=sd,
maximum.bias=maxbias, conf.low=NA, conf.high=NA,
conf.low.onesided=NA, conf.high.onesided=NA,
Expand Down
10 changes: 6 additions & 4 deletions tests/testthat/test_frd.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,10 +86,12 @@ test_that("FRD interface", {
M=as.numeric(p1.1$coefficients[16:17]),
kern="triangular", opt.criterion="OCI",
T0=p1.1$coefficients$estimate)
eo <- c(paste0(" Fuzzy RD parameter -0.5536",
" 0.674 0.4361 (-2.109, 1.002)"),
eo <- c(paste0("retired -0.5536 0.674 0.4361",
" (-2.109, 1.002)"),
"Number of effective observations: 175",
"Maximal leverage for fuzzy RD parameter: 0.03235",
"First stage estimate: 0.2562 ")
expect_equal(capture.output(print(p2, digits=4))[c(9, 14:15)],
"First stage estimate: 0.2562 ", "P-value: 0.5018 ",
"elig_year -0.0125 0.0403")
expect_equal(capture.output(print(p2, digits=4))[c(9, 12:14, 17, 25)],
eo)
})
9 changes: 7 additions & 2 deletions tests/testthat/test_lpp.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,13 @@ test_that("Optimal bandwidth calculations", {

r1 <- RDHonest(voteshare~margin, data=lee08, subset=margin>0, M=2*Mh,
opt.criterion="MSE", point.inference=TRUE)
r <- capture.output(print(r1, digits=4))
expect_equal(r[11], "Bandwidth: 13.41, Kernel: triangular")
r <- capture.output(print(r1, digits=6))
expect_equal(r[c(8, 11, 12)],
c(paste0("(Intercept) 52.3928 0.880413",
" 0.492905 (50.4288, 54.3568)"),
"Number of effective observations: 662.044",
paste0("Maximal leverage for value of",
" conditional mean: 0.00936479")))

re$d$sigma2 <- NULL
r2 <- OptBW(re$d, M=2*Mh, opt.criterion="MSE")
Expand Down
32 changes: 18 additions & 14 deletions tests/testthat/test_rd.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,12 +78,12 @@ test_that("Honest inference in Lee and LM data", {
expect_equal(r1$coefficients$estimate, -1.1982581306)
expect_equal(c(r1$coefficients$std.error, r2$coefficients$std.error),
c(0.6608206598, 0.6955768728))
m1 <- paste0(" Sharp RD parameter -1.198 ",
"0.6608 4.796 (-7.081, 4.684)")
expect_equal(r1o[8], m1)
expect_equal(r1o[8], paste0("I(povrate>0) -1.198 0.6608",
" 4.796 (-7.081, 4.684)"))
expect_equal(r1o[10], "Onesided CIs: (-Inf, 4.684), (-7.081, Inf)")
expect_equal(r1o[11], "Number of effective observations: 954")
expect_equal(r2o[10], "Onesided CIs: (-Inf, 4.742), (-7.138, Inf)")
expect_equal(r1o[16], "24 observations with missing values dropped")
expect_equal(r1o[22], "24 observations with missing values dropped")

d <- RDHonest(mortHS~povrate, data=headst, h=10, M=2)$d
d <- PrelimVar(d, se.initial="Silverman")
Expand Down Expand Up @@ -138,7 +138,7 @@ test_that("Honest inference in Lee and LM data", {
r <- RDHonest(voteshare ~ margin, data=lee08, kern="optimal", M=0.04,
opt.criterion="FLCI", se.method="supplied.var")
ro <- capture.output(print(r, digits=8))
expect_equal(ro[8], paste0(" Sharp RD parameter 5.8864375 1.4472117",
expect_equal(ro[8], paste0("margin 5.8864375 1.4472117",
" 0.80247208 (2.6646087, 9.1082664)"))
## Try nn
r1 <- RDHonest(voteshare ~ margin, data=lee08, kern="optimal", M=0.04,
Expand Down Expand Up @@ -198,17 +198,20 @@ test_that("Honest inference in Lee and LM data", {
expect_message(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, 16)],
c("Bandwidth: 3.98048, Kernel: uniform",
expect_equal(r1[c(8, 11, 12, 22)],
c(paste0("I(povrate>0) -3.17121 1.44439",
" 0.759228 (-6.35207, 0.00964571)"),
"Number of effective observations: 239",
"Maximal leverage for sharp RD parameter: 0.019615",
"24 observations with missing values dropped"))
expect_message(r2 <- RDHonest(mortHS ~ povrate, data=headst,
kern="epanechnikov",
point.inference=TRUE, cutoff=4))
expect_equal(capture.output(print(r2, digits=6))[c(11, 12, 16)],
c("Bandwidth: 9.42625, Kernel: epanechnikov",
expect_equal(capture.output(print(r2, digits=6))[c(8, 11, 19)],
c(paste0("(Intercept) 2.63181 0.300132 ",
"0.152587 (1.97505, 3.28856)"),
"Number of effective observations: 373.642",
"24 observations with missing values dropped"))
" 2.631805 -0.000162 "))
})

test_that("BME CIs match paper", {
Expand Down Expand Up @@ -237,10 +240,11 @@ test_that("BME CIs match paper", {
r6 <- RDHonestBME(duration~age, data=rebp, subset = (period==1 & female==0),
order=3, h=1, cutoff=50, alpha=0.1)
r5 <- capture.output(print(r5, digits=5))
est <- paste0(" Sharp RD parameter 14.798 ",
"2.234 24.315 (-31.823, 60.724)")
expect_equal(r5[8], est)
expect_equal(r5[10], "Onesided CIs: (-Inf, 57.249), (-28.236, Inf)")
expect_equal(r5[c(8, 10, 12)],
c(paste0("I(x >= 0)TRUE 14.798 2.234 ",
"24.315 (-31.823, 60.724)"),
"Onesided CIs: (-Inf, 57.249), (-28.236, Inf)",
"Maximal leverage for sharp RD parameter: 0.00039197"))
expect_equal(as.numeric(r6$coefficients[c(2:6, 10)]),
c(12.206023620, 8.878539149, 17.030704793, -24.727726605,
46.856041887, 3030))
Expand Down

0 comments on commit 661fe6c

Please sign in to comment.