begin outfactoring of create_Durbin

Author: rsbivand

R/SLX_WX.R

@@ -46,7 +46,11 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
         if (is.null(weights)) weights <- rep(as.numeric(1), n)
         if (any(is.na(weights))) stop("NAs in weights")
         if (any(weights < 0)) stop("negative weights")
-        dvars <- c(NCOL(x), 0L)
+#        dvars <- c(NCOL(x), 0L)
+        if (!(isTRUE(Durbin) || is.formula(Durbin))) {
+            stop("Durbin argument neither TRUE nor formula")
+        } else {
+        if (Sys.getenv("SPATIALREG_CREATE_DURBIN") == "") {
         prefix <- "lag"
         if (isTRUE(Durbin)) {
             if (have_factor_preds) warn_factor_preds(have_factor_preds)
@@ -94,7 +98,7 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
             wxn <- substring(colnames(WX), nchar(prefix)+2,
                 nchar(colnames(WX)))
             zero_fill <- length(xn) + (which(!(xn %in% wxn)))
-        } else stop("Durbin argument neither TRUE nor formula")
+        } 
         dvars <- c(NCOL(x), NCOL(WX))
         if (is.formula(Durbin)) {
             attr(dvars, "f") <- Durbin
@@ -104,6 +108,21 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
         }
 	x <- cbind(x, WX)
 	rm(WX)
+        } else { # SPATIALREG_CREATE_DURBIN
+            res <- create_Durbin(Durbin=Durbin, 
+                have_factor_preds=have_factor_preds, x=x, listw=listw,
+                zero.policy=zero.policy, data=data, na.act=na.act)
+            x <- res$x
+            dvars <- res$dvars
+            inds <-attr(dvars, "inds") 
+            xn <- attr(dvars, "xn")
+            wxn <- attr(dvars, "wxn")
+            zero_fill <- attr(dvars, "zero_fill")
+            formula_durbin_factors <- attr(dvars, "formula_durbin_factors")
+            attr(dvars, "xn") <- NULL
+            attr(dvars, "wxn") <- NULL
+        }
+        } 
 #        WX <- create_WX(x, listw, zero.policy=zero.policy, prefix="lag")
 #        x <- cbind(x, WX)
 # 180128 Mark L. Burkey summary.lm error for SlX object
@@ -437,4 +456,67 @@ create_WX <- function(x, listw, zero.policy=NULL, prefix="") {
         WX
 }
 
+create_Durbin <- function(Durbin, have_factor_preds, x, listw, zero.policy, 
+        data, na.act) {
+        prefix <- "lag"
+        if (isTRUE(Durbin)) {
+            if (have_factor_preds) warn_factor_preds(have_factor_preds)
+            WX <- create_WX(x, listw, zero.policy=zero.policy,
+               prefix=prefix)
+        } else if (is.formula(Durbin)) {
+            data1 <- data
+            if (!is.null(na.act) && (inherits(na.act, "omit") ||
+                inherits(na.act, "exclude"))) {
+                data1 <- data1[-c(na.act),]
+            }
+            dmf <- lm(Durbin, data1, na.action=na.fail, 
+	        method="model.frame")
+	    formula_durbin_factors <- have_factor_preds_mf(dmf)
+            if (formula_durbin_factors)
+                warn_factor_preds(formula_durbin_factors)
+#	    dmf <- lm(Durbin, data, na.action=na.action, 
+#	         method="model.frame")
+            fx <- try(model.matrix(Durbin, dmf), silent=TRUE)
+            if (inherits(fx, "try-error")) 
+                 stop("Durbin variable mis-match")
+            WX <- create_WX(fx, listw, zero.policy=zero.policy,
+                prefix=prefix)
+            inds <- match(substring(colnames(WX), 5,
+	        nchar(colnames(WX))), colnames(x))
+            if (anyNA(inds)) {
+              wna <- which(is.na(inds)) #TR: continue if Durbin has intercept, but formula has not
+              if (length(wna) == 1 && grepl("Intercept", colnames(WX)[wna])
+                 && attr(terms(formula), "intercept") == 0
+                 && attr(terms(Durbin), "intercept") == 1) {
+                inds <- inds[-wna]
+              } else{
+                stop("WX variables not in X: ",
+                     paste(substring(colnames(WX), 5,
+                     nchar(colnames(WX)))[is.na(inds)], collapse=" "))
+              }
+            } 
+            icept <- grep("(Intercept)", colnames(x))
+            iicept <- length(icept) > 0L
+            if (iicept) {
+                xn <- colnames(x)[-1]
+            } else {
+                xn <- colnames(x)
+            }
+            wxn <- substring(colnames(WX), nchar(prefix)+2,
+                nchar(colnames(WX)))
+            zero_fill <- length(xn) + (which(!(xn %in% wxn)))
+        } 
+        dvars <- c(NCOL(x), NCOL(WX))
+        if (is.formula(Durbin)) {
+            attr(dvars, "f") <- Durbin
+            attr(dvars, "inds") <- inds
+            attr(dvars, "xn") <- xn
+            attr(dvars, "wxn") <- wxn
+            attr(dvars, "zero_fill") <- zero_fill
+            attr(dvars, "formula_durbin_factors") <- formula_durbin_factors
+        }
+	x <- cbind(x, WX)
+	rm(WX)
+        list(x=x, dvars=dvars)
+}
 

inst/tinytest/test_Durbin_factor.R

@@ -4,10 +4,21 @@ lw <- spdep::nb2listw(COL.nb)
 COL.OLD$fEW <- factor(COL.OLD$EW)
 COL.OLD$fDISCBD <- ordered(cut(COL.OLD$DISCBD, c(0, 1.5, 3, 4.5, 6)))
 f <- formula(CRIME ~ INC + HOVAL + fDISCBD*fEW)
+Sys.setenv("SPATIALREG_CREATE_DURBIN"="")
 expect_warning(COL.SLX0 <- lmSLX(f, data=COL.OLD, lw, Durbin=TRUE))
 expect_warning(COL.SLX1 <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fDISCBD*fEW))
 expect_warning(COL.SLX2 <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fEW))
 expect_silent(COL.SLX3 <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL))
+Sys.setenv("SPATIALREG_CREATE_DURBIN"="0")
+expect_warning(COL.SLX0a <- lmSLX(f, data=COL.OLD, lw, Durbin=TRUE))
+expect_warning(COL.SLX1a <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fDISCBD*fEW))
+expect_warning(COL.SLX2a <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fEW))
+expect_silent(COL.SLX3a <- lmSLX(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL))
+Sys.setenv("SPATIALREG_CREATE_DURBIN"="")
+expect_true(isTRUE(all.equal(COL.SLX0, COL.SLX0a)))
+expect_true(isTRUE(all.equal(COL.SLX1, COL.SLX1a)))
+expect_true(isTRUE(all.equal(COL.SLX2, COL.SLX2a)))
+expect_true(isTRUE(all.equal(COL.SLX3, COL.SLX3a)))
 expect_warning(COL.err0 <- errorsarlm(f, data=COL.OLD, lw, Durbin=TRUE))
 expect_warning(COL.err1 <- errorsarlm(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fDISCBD*fEW))
 expect_warning(COL.err2 <- errorsarlm(f, data=COL.OLD, lw, Durbin=~ INC + HOVAL + fDISCBD))