Skip to content

Commit

Permalink
Less insane naming scheme, fix linting
Browse files Browse the repository at this point in the history
  • Loading branch information
@kolesarm
kolesarm committed Jan 28, 2024
1 parent 3e4c674 commit bda24e5
Show file tree
Hide file tree
Showing 29 changed files with 539 additions and 548 deletions.
8 changes: 7 additions & 1 deletion .lintr
View file
Original file line number Diff line number Diff line change
@@ -1 +1,7 @@
linters: linters_with_defaults(line_length_linter(80), commented_code_linter=NULL, infix_spaces_linter=NULL, object_name_linter = NULL, spaces_left_parentheses_linter=NULL)
linters: linters_with_defaults(
line_length_linter(80),
infix_spaces_linter=NULL,
commented_code_linter=NULL,
object_name_linter(c("CamelCase", "camelCase", "dotted.case", "snake_case")),
indentation_linter(indent=4L)
)
2 changes: 1 addition & 1 deletion NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,6 @@ S3method(print,RDResults)
export(CVb)
export(RDHonest)
export(RDHonestBME)
export(RDScatter)
export(RDSmoothnessBound)
export(RDTEfficiencyBound)
export(plot_RDscatter)
31 changes: 16 additions & 15 deletions R/Cbound.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,9 @@
#' @references{
#'
#' \cite{Michal Kolesár and Christoph Rothe. Inference in regression
#' discontinuity designs with a discrete running variable. American Economic
#' Review, 108(8):2277—-2304, August 2018. \doi{10.1257/aer.20160945}}
#' discontinuity designs with a discrete running variable.
#' American Economic Review, 108(8):2277—-2304,
#' August 2018. \doi{10.1257/aer.20160945}}
#'
#' }
#' @examples
Expand All @@ -43,7 +44,7 @@
#' @export
RDSmoothnessBound <- function(object, s, separate=FALSE, multiple=TRUE,
alpha=0.05, sclass="H") {
d <- NPRPrelimVar.fit(object$data, se.initial="EHW")
d <- PrelimVar(object$data, se.initial="EHW")

## Curvature estimate based on jth set of three points closest to zero
Dk <- function(Y, X, xu, s2, j) {
Expand All @@ -52,16 +53,16 @@ RDSmoothnessBound <- function(object, s, separate=FALSE, multiple=TRUE,
I3 <- X >= xu[3*j*s- s+1] & X <= xu[3*j*s]

lam <- (mean(X[I3])-mean(X[I2])) / (mean(X[I3])-mean(X[I1]))
den <- if (sclass=="T") {
(1-lam)*mean(X[I3]^2) + lam*mean(X[I1]^2) + mean(X[I2]^2)
} else {
(1-lam)*mean(X[I3]^2) + lam*mean(X[I1]^2) - mean(X[I2]^2)
}
if (sclass=="T") {
den <- (1-lam)*mean(X[I3]^2) + lam*mean(X[I1]^2) + mean(X[I2]^2)
} else {
den <- (1-lam)*mean(X[I3]^2) + lam*mean(X[I1]^2) - mean(X[I2]^2)
}
## Delta is lower bound on M by lemma S2 in Kolesar and Rothe
Del <- 2*(lam*mean(Y[I1])+(1-lam)*mean(Y[I3])-mean(Y[I2])) / den
Del <- 2 * (lam*mean(Y[I1]) + (1-lam) * mean(Y[I3])-mean(Y[I2])) / den
## Variance of Delta
VD <- 4*(lam^2*mean(s2[I1])/sum(I1) +
(1-lam)^2*mean(s2[I3])/sum(I3) + mean(s2[I2])/sum(I2)) / den^2
VD <- 4 * (lam^2*mean(s2[I1])/sum(I1) + (1-lam)^2*mean(s2[I3]) /
sum(I3) + mean(s2[I2])/sum(I2)) / den^2
c(Del, sqrt(VD), mean(Y[I1]), mean(Y[I2]), mean(Y[I3]), range(X[I1]),
range(X[I2]), range(X[I3]))
}
Expand All @@ -72,8 +73,8 @@ RDSmoothnessBound <- function(object, s, separate=FALSE, multiple=TRUE,
Dpj <- function(j) Dk(d$Y[d$p], d$X[d$p], xp, d$sigma2[d$p], j)
Dmj <- function(j) Dk(d$Y[d$m], abs(d$X[d$m]), xm, d$sigma2[d$m], j)

Sp <- floor(length(xp)/(3*s))
Sm <- floor(length(xm)/(3*s))
Sp <- floor(length(xp) / (3*s))
Sm <- floor(length(xm) / (3*s))
if (min(Sp, Sm) == 0)
stop("Value of s is too big")

Expand Down Expand Up @@ -111,8 +112,8 @@ RDSmoothnessBound <- function(object, s, separate=FALSE, multiple=TRUE,
}
list(estimate=hatM, conf.low=lower,
diagnostics=c(Delta=maxt[1], sdDelta=maxt[2],
y1=maxt[3], y2=maxt[4], y3=maxt[5],
I1=maxt[6:7], I2=maxt[8:9], I3=maxt[10:11]))
y1=maxt[3], y2=maxt[4], y3=maxt[5],
I1=maxt[6:7], I2=maxt[8:9], I3=maxt[10:11]))
}

if (separate) {
Expand Down
159 changes: 0 additions & 159 deletions R/NPR_lp.R
View file

This file was deleted.

52 changes: 24 additions & 28 deletions R/NPRfunctions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,8 @@ LPReg <- function(X, Y, h, K, order=1, se.method=NULL, sigma2, J=3,

## Squared residuals, allowing for multivariate Y
sig <- function(r) {
r[, rep(seq_len(ncol(r)), each=ncol(r))] *
r[, rep(seq_len(ncol(r)), ncol(r))]
r[, rep(seq_len(ncol(r)), each=ncol(r))] * r[, rep(seq_len(ncol(r)),
ncol(r))]
}
## For RD, compute variance separately on either side of cutoff
signn <- function(X) {
Expand All @@ -45,15 +45,15 @@ LPReg <- function(X, Y, h, K, order=1, se.method=NULL, sigma2, J=3,
wgt_unif <- ((weights)*R %*% solve(crossprod(R, weights*R)))[, 1]
eff.obs <- length(X)*sum(wgt_unif^2)/sum(wgt^2)
## Variance
V <- if (is.null(clusterid)) {
colSums(as.matrix(wgt^2 * hsigma2))
} else if (se.method == "supplied.var") {
colSums(as.matrix(wgt^2 * hsigma2))+
rho*(sum(tapply(wgt, clusterid, sum)^2)-sum(wgt^2))
} else {
us <- apply(wgt*res, 2, function(x) tapply(x, clusterid, sum))
as.vector(crossprod(us))
}
if (is.null(clusterid)) {
V <- colSums(as.matrix(wgt^2 * hsigma2))
} else if (se.method == "supplied.var") {
V <- colSums(as.matrix(wgt^2 * hsigma2))+
rho * (sum(tapply(wgt, clusterid, sum)^2)-sum(wgt^2))
} else {
us <- apply(wgt*res, 2, function(x) tapply(x, clusterid, sum))
V <- as.vector(crossprod(us))
}

list(theta=beta[1, ], sigma2=hsigma2, res=res,
var=V, w=wgt, eff.obs=eff.obs)
Expand All @@ -67,14 +67,14 @@ LPReg <- function(X, Y, h, K, order=1, se.method=NULL, sigma2, J=3,
## Calculate fuzzy or sharp RD estimate, or estimate of a conditional mean at a
## point (depending on the class of \code{d}), and its variance using local
## polynomial regression of order \code{order}.
NPRreg.fit <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
if (!is.function(kern))
kern <- EqKern(kern, boundary=FALSE, order=0)
## Keep only positive kernel weights
W <- if (h<=0) 0*d$X else kern(d$X/h) # kernel weights
if(!is.null(d$sigma2)) d$sigma2 <- as.matrix(d$sigma2)
if (!is.null(d$sigma2)) d$sigma2 <- as.matrix(d$sigma2)
r <- LPReg(d$X[W>0], as.matrix(d$Y)[W>0, ], h, kern, order, se.method,
d$sigma2[W>0, ], J, weights=d$w[W>0], RD=(d$class!="IP"),
d$sigma2[W>0, ], J, weights=d$w[W>0], RD = (d$class!="IP"),
d$rho, d$clusterid[W>0])
sigma2 <- matrix(NA, nrow=length(W), ncol=NCOL(d$Y)^2)
sigma2[W>0, ] <- r$sigma2
Expand All @@ -88,8 +88,8 @@ NPRreg.fit <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
if (d$class=="FRD") {
ret$fs <- r$theta[2]
ret$estimate <- r$theta[1]/r$theta[2]
ret$se <- sqrt(sum(c(1, -ret$estimate, -ret$estimate, ret$estimate^2) *
r$var) / ret$fs^2)
ret$se <- sqrt(sum(c(1, -ret$estimate, -ret$estimate,
ret$estimate^2) * r$var) / ret$fs^2)
}
ret
}
Expand All @@ -100,26 +100,22 @@ NPRreg.fit <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
## Use global quartic regression to estimate a bound on the second derivative
## for inference under under second order Hölder class. For RD, use a separate
## regression on either side of the cutoff
NPR_MROT.fit <- function(d) {
MROT <- function(d) {
if (d$class=="SRD") {
max(NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$p], X=d$X[d$p]), 0, "IP")),
NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$m], X=d$X[d$m]), 0, "IP")))
max(MROT(NPRData(data.frame(Y=d$Y[d$p], X=d$X[d$p]), 0, "IP")),
MROT(NPRData(data.frame(Y=d$Y[d$m], X=d$X[d$m]), 0, "IP")))
} else if (d$class=="FRD") {
c(M1=max(NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$p, 1], X=d$X[d$p]), 0,
"IP")),
NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$m, 1], X=d$X[d$m]), 0,
"IP"))),
M2=max(NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$p, 2], X=d$X[d$p]), 0,
"IP")),
NPR_MROT.fit(NPRData(data.frame(Y=d$Y[d$m, 2], X=d$X[d$m]), 0,
"IP"))))
c(M1=max(MROT(NPRData(data.frame(Y=d$Y[d$p, 1], X=d$X[d$p]), 0, "IP")),
MROT(NPRData(data.frame(Y=d$Y[d$m, 1], X=d$X[d$m]), 0, "IP"))),
M2=max(MROT(NPRData(data.frame(Y=d$Y[d$p, 2], X=d$X[d$p]), 0, "IP")),
MROT(NPRData(data.frame(Y=d$Y[d$m, 2], X=d$X[d$m]), 0, "IP"))))
} else if (d$class=="IP") {
## STEP 1: Estimate global polynomial regression
r1 <- unname(stats::lm(d$Y ~ 0 + outer(d$X, 0:4, "^"))$coefficients)
f2 <- function(x) abs(2*r1[3]+6*x*r1[4]+12*x^2*r1[5])
## maximum occurs either at endpoints, or else at the extremum,
## -r1[4]/(4*r1[5]), if the extremum is in the support
f2e <- if(abs(r1[5])<=1e-10) Inf else -r1[4]/(4*r1[5])
f2e <- if (abs(r1[5])<=1e-10) Inf else -r1[4] / (4*r1[5])
M <- max(f2(min(d$X)), f2(max(d$X)))
if (min(d$X) < f2e && max(d$X) > f2e) M <- max(f2(f2e), M)

Expand Down
Loading

0 comments on commit bda24e5

Please sign in to comment.