Allow for colllinear covariates

DESCRIPTION

@@ -35,5 +35,5 @@ Config/testthat/edition: 3
 Language: en-US
 URL: https://github.com/kolesarm/RDHonest
 BugReports: https://github.com/kolesarm/RDHonest/issues
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr

R/NPRfunctions.R

@@ -14,19 +14,32 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
     nZ <- c("(Intercept)", colnames(d$X), paste0("I(", colnames(d$X), "^",
                                                  seq_len(order), ")")[-1])
     colnames(Z) <- nZ[seq_len(order+1)]
+    Lz <- NCOL(Z)
     if (!inherits(d, "IP")) {
         ZZ <- (X>=0)*Z
         colnames(ZZ) <- c(paste0("I(", colnames(d$X), ">0)"),
                           paste0(paste0("I(", colnames(d$X), ">0):"),
                                  colnames(Z))[-1])
         Z <- cbind(ZZ, Z, d$covs)
+        Lz <- 2*Lz
     }
     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))) {
+    be <- as.matrix(r0$coefficients)
+    if (any(is.na(be[1:Lz, ]))) {
         return(list(estimate=0, se=NA, est_w=W*0, sigma2=NA*d$Y, eff.obs=0,
-                    fs=NA, lm=r0))
+                    fs=NA, lm=r0, Z=Z))
+    }
+    ## If the collinearity comes from covariates, drop them
+    if (any(is.na(be[-(1:Lz), ]))) {
+        Z <- Z[, !is.na(rowSums(be))]
+        message("The following covariates are collinear",
+                " and are dropped:\n",
+                paste(names(which(is.na(rowSums(be[-(1:Lz), , drop=FALSE])))),
+                            collapse=", "))
+        r0 <- stats::lm.wfit(x=Z, y=d$Y, w=W)
     }
+
+    class(r0) <- c(if (ncol(d$Y)>1) "mlm", "lm")
     wgt <- W*0
     ok <- W!=0
     wgt[ok] <- solve(qr.R(r0$qr), t(sqrt(W[W>0])*qr.Q(r0$qr)))[1, ]
@@ -70,7 +83,7 @@ NPReg <- function(d, h, kern="triangular", order=1, se.method="nn", J=3) {
         V <- as.vector(crossprod(us))
     }
     ret <- list(estimate=r0$coefficients[1], se=sqrt(V[1]), est_w=wgt,
-                sigma2=hsigma2, eff.obs=eff.obs, fs=NA, lm=r0)
+                sigma2=hsigma2, eff.obs=eff.obs, fs=NA, lm=r0, Z=Z)
 
     if (inherits(d, "FRD")) {
         ret$fs <- r0$coefficients[1, 2]

R/RDHonest.R

@@ -239,8 +239,12 @@ covariate_adjust <- function(d, kern, h) {
     d0$Y <- d0$Y_unadj
     r <- NPReg(d0, h, kern, order=1, se.method="EHW")
     be <- as.matrix(r$lm$coefficients)
-    L <- NCOL(d$covs)
-    d$Y <-d$Y_unadj-d$covs %*% be[seq_len(L)+NROW(be)-L, , drop=FALSE]
+    if (r$eff.obs==0)
+        stop("Too few effective observations to perform covariate adjustment")
+    ## Original
+    ## L <- NCOL(d$covs)
+    ## d$Y <-d$Y_unadj-d$covs %*% be[seq_len(L)+NROW(be)-L, , drop=FALSE]
+    d$Y <-d$Y_unadj-r$Z[, -(1:4)] %*% be[-(1:4), , drop=FALSE]
     d
 }
 

tests/testthat/test_covariates.R

@@ -54,3 +54,27 @@ test_that("Test covariates", {
                                    kern=function(u) pmax(1-abs(u), 0)))
     expect_equal(r00$coefficients, r00$coefficients)
 })
+
+test_that("Test collinear covariates", {
+
+    df <- headst[complete.cases(headst), ]
+    covlist <- "urban+black"
+    fh1 <- as.formula(paste0("mortHS ~ povrate|", covlist))
+    r0 <- RDHonest(fh1, data=df)
+    Z <- cbind(df$urban, df$black, df$urban)
+    expect_message(r1 <- RDHonest(df$mortHS~df$povrate|Z, M=r0$coefficients$M))
+    expect_lt(max(abs(r1$coefficients[2:13]-r0$coefficients[2:13])), 1e-15)
+    Z <- cbind(df$urban, df$urban, df$black, df$black, df$urban+df$black)
+    expect_message(r2 <- RDHonest(df$mortHS~df$povrate|Z))
+    expect_lt(max(abs(r2$coefficients[2:13]-r0$coefficients[2:13])), 1e-15)
+    expect_error(RDHonest(df$mortHS~df$povrate|Z, h=0.05, kern="uniform"))
+
+    ## Fuzzy
+    df <- rcp[1:1000, ]
+    Z <- cbind(df$food, df$food+df$family_size, df$family_size)
+    r0 <- RDHonest(log(c)|retired ~ elig_year|Z[, 1:2], data=df,
+                   weights=survey_year)
+    expect_message(r1 <- RDHonest(log(c)|retired ~ elig_year|Z, data=df,
+                                  weights=survey_year))
+    expect_lt(max(abs(r1$coefficients[2:13]-r0$coefficients[2:13])), 1e-15)
+})