[Gmm-commits] r193 - in pkg/causalOTLSE: . R

[noreply at r-forge.r-project.org]

Author: chaussep
Date: 2022-06-29 20:52:57 +0200 (Wed, 29 Jun 2022)
New Revision: 193

Modified:
   pkg/causalOTLSE/DESCRIPTION
   pkg/causalOTLSE/NAMESPACE
   pkg/causalOTLSE/R/otlse.R
Log:
Improved efficiency

Modified: pkg/causalOTLSE/DESCRIPTION
===================================================================
--- pkg/causalOTLSE/DESCRIPTION	2022-06-29 14:14:07 UTC (rev 192)
+++ pkg/causalOTLSE/DESCRIPTION	2022-06-29 18:52:57 UTC (rev 193)
@@ -1,6 +1,6 @@
 Package: causalOTLSE
 Version: 0.1-0
-Date: 2022-06-07
+Date: 2022-06-29
 Title: Optimal Thresholding Least Squares Inference for Causal Effects
 Authors at R: c(person("Pierre Chausse", "Developer", role = c("aut", "cre"),
                      email = "pchausse at uwaterloo.ca"),
@@ -8,10 +8,6 @@
 	             email = "giurcanu at uchicago.edu"))
 Description: This package includes tools to measure causal effects using least squares regressions. The number of piecewise polynomials is selected by some information criteria.
 Depends: R (>= 4.0.0)
-Imports: stats, splines, car, sandwich, cvTools, graphics
+Imports: stats, splines, sandwich, cvTools, graphics
 License: GPL (>= 2)
 NeedsCompilation: no
-Packaged: 2022-06-07 19:15:24 UTC; pierrechausse
-Author: Pierre Chausse Developer [aut, cre],
-  Mihai Giurcanu Developer [aut]
-Maintainer: Pierre Chausse Developer <pchausse at uwaterloo.ca>

Modified: pkg/causalOTLSE/NAMESPACE
===================================================================
--- pkg/causalOTLSE/NAMESPACE	2022-06-29 14:14:07 UTC (rev 192)
+++ pkg/causalOTLSE/NAMESPACE	2022-06-29 18:52:57 UTC (rev 193)
@@ -1,6 +1,5 @@
 importFrom(stats, quantile, lm, predict, coef, vcov, printCoefmat, cov, pnorm,
-           terms, model.matrix, as.formula, model.frame)
-importFrom(car, linearHypothesis)
+           terms, model.matrix, as.formula, model.frame, lm.fit, na.omit, pf)
 importFrom(graphics, grid, legend, matplot)
 importFrom(sandwich, vcovHC)
 importFrom(cvTools, cvFolds)

Modified: pkg/causalOTLSE/R/otlse.R
===================================================================
--- pkg/causalOTLSE/R/otlse.R	2022-06-29 14:14:07 UTC (rev 192)
+++ pkg/causalOTLSE/R/otlse.R	2022-06-29 18:52:57 UTC (rev 193)
@@ -128,7 +128,7 @@
     data$Xf0 <- ppSplines(form=formX, data=data, pFact=ppow,
                           method=splineMet, subGroup=id0, minZeroProp=mZeroProp)
     data$Xf1 <- ppSplines(form=formX, data=data, pFact=ppow,
-                          method=splineMet, subGroup=!id0, minZeroProp=mZeroProp)    
+                          method=splineMet, subGroup=!id0, minZeroProp=mZeroProp)
     fit <- lm(formY, data)
     knots0 <- attr(data$Xf0, "knots")
     knots1 <- attr(data$Xf1, "knots")
@@ -136,21 +136,26 @@
     p1 <- attr(data$Xf1, "p")
     cn0 <- attr(data$Xf0, "colnames")
     cn1 <- attr(data$Xf1, "colnames")
+    v <- vcovHC(fit, type = HCtype)
     pval0 <- lapply(1:length(p0), function(i){
         if (p0[i] == 1)
             return(NA)
         sapply(1:(p0[i] - 1),  function(j){
-            t <- paste("Xf0", cn0[[i]][j], "=Xf0", cn0[[i]][j+1], sep="")
-            linearHypothesis(fit, t, vcov. = vcovHC(fit, type = HCtype))[2, 4]
+            t <- c(paste("Xf0", cn0[[i]][j], sep=""),
+                   paste("Xf0", cn0[[i]][j+1], sep=""))
+            .lintest(fit, t[1], t[2], v)
         })})
+    names(pval0) <- NULL
     names(pval0) <- names(p0)
     pval1 <- lapply(1:length(p1), function(i){
         if (p1[i] == 1)
             return(NA)        
         sapply(1:(p1[i] - 1),  function(j){
-            t <- paste("Xf1", cn1[[i]][j], "=Xf1", cn1[[i]][j+1], sep="")
-            linearHypothesis(fit, t, vcov. = vcovHC(fit, type = HCtype))[2, 4]
+            t <- c(paste("Xf1", cn0[[i]][j], sep=""),
+                   paste("Xf1", cn0[[i]][j+1], sep=""))
+            .lintest(fit, t[1], t[2], v)            
         })})
+    names(pval1) <- NULL
     names(pval1) <- names(p1)
     list(pval0 = pval0, pval1 = pval1, p0 = p0, p1 = p1, knots0 = knots0, 
          knots1 = knots1, Z=Z, Y=Y, formX=formX, formY=formY,
@@ -157,6 +162,19 @@
          treatment=colnames(Z))
 }
 
+.lintest <- function(obj, c1, c2, v)
+{
+    b <- coef(obj)
+    b <- na.omit(b)
+    c1 <- which(names(b) == c1)
+    c2 <- which(names(b) == c2)
+    s2 <- v[c1,c1]+v[c2,c2]-2*v[c1,c2]
+    ans <- 1-pf((b[c1]-b[c2])^2/s2, 1, obj$df)
+    names(ans) <- NULL
+    ans
+}
+
+
 selASY <- function (form, data, pFact = 0.3, splineMet = c("manual", "bs"),
                     HCtype="HC", mZeroProp=0.1, minPV=function(p) 1/(p*log(p)))
 {
@@ -204,22 +222,15 @@
     X1 <- as.matrix(X1)
     myK <- floor(log(n))
     cv.outi <- cvFolds(n, K = myK)
-    mspe_pred <- rep(NA, myK)    
+    mspe_pred <- rep(NA, myK)
+    X <- cbind(1-Z,Z,X0,X1)
     for(k in 1 : myK)
     {
         id.train <- cv.outi$subsets[cv.outi$which != k]
         id.valid <- cv.outi$subsets[cv.outi$which == k]
-        train.datak <- list(Yk = Y[id.train], Zk = Z[id.train], 
-                            Xf0ik = X0[id.train, , drop = FALSE],
-                            Xf1ik = X1[id.train, , drop = FALSE])
-        valid.datak <- list(Yk = Y[id.valid], Zk = Z[id.valid], 
-                            Xf0ik = X0[id.valid, , drop = FALSE],
-                            Xf1ik = X1[id.valid, , drop = FALSE])
-        lm.outk <- lm(Yk ~ 0 + factor(Zk) + Xf0ik + Xf1ik, data = train.datak)
-        pred.outk <- predict(lm.outk, newdata = valid.datak,
-                             type = "response")
-        mspe_pred[k] <- mean((valid.datak$Yk - pred.outk)^2,
-                             na.rm = TRUE)
+        lm.outk <- lm.fit(X[id.train,],Y[id.train])
+        pred.outk <- c(X[id.valid,]%*%lm.outk$coefficients)
+        mspe_pred[k] <- mean((Y[id.valid]-pred.outk)^2, na.rm=TRUE)
     }    
     mean(mspe_pred, na.rm = TRUE)
 }