add warnings to SD tests, including ordered factors

Author: rsbivand

NAMESPACE

@@ -12,7 +12,7 @@ importFrom(stats, influence.measures, lag, punif, lm, var, integrate,
              "model.frame", "model.matrix", "model.response", "na.fail",
              "na.omit", "naresid", "optimise", "printCoefmat", "resid",
              "terms", "uniroot", "weighted.residuals", "weights", "median",
-             "pbinom")
+             "pbinom", "C")
 
 importFrom(deldir, deldir)
 importFrom(boot, boot)

R/SD.RStests.R

@@ -5,8 +5,9 @@ is.formula <- function(x){
    inherits(x,"formula")
 }
 
-create_X0 <- function(X, listw, Durbin=TRUE, data=NULL, na.act=NULL) {
+create_X0 <- function(X, listw, Durbin=TRUE, data=NULL, na.act=NULL, have_factor_preds=FALSE) {
         if (isTRUE(Durbin)) {
+            if (have_factor_preds) warn_factor_preds(have_factor_preds)
             n <- NROW(X)
 	    m <- NCOL(X)
 	    # check if there are enough regressors
@@ -43,6 +44,9 @@ create_X0 <- function(X, listw, Durbin=TRUE, data=NULL, na.act=NULL) {
             }
             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")
             X0 <- try(model.matrix(Durbin, dmf), silent=TRUE)
@@ -106,14 +110,19 @@ SD.RStests <- function(model, listw, zero.policy=attr(listw, "zero.policy"), tes
 		warning("Spatial weights matrix not row standardized")
 
         if (is.formula(Durbin)) {
-            dt <- try(eval(model$call[["data"]]), silent=TRUE)
+            dt <- try(as.data.frame(model$model), silent=TRUE)
+#            dt <- try(eval(model$call[["data"]]), silent=TRUE)
             if (inherits(dt, "try-error") || !is.data.frame(dt))
-                stop("data object used to fit linear model not available for formula Durbin")
+                stop("model.frame object used to fit linear model not available for formula Durbin")
+
         }
 
-	y <- model.response(model.frame(model))
-	X <- model.matrix(terms(model), model.frame(model))
-        X0 <- create_X0(X=X, listw=listw, Durbin=Durbin, data=dt, na.act=na.act)
+	mf <- model.frame(model)
+        y <- model.response(mf)
+	X <- model.matrix(terms(model), mf)
+        have_factor_preds <- have_factor_preds_mf(mf)
+        X0 <- create_X0(X=X, listw=listw, Durbin=Durbin, data=dt, na.act=na.act,
+            have_factor_preds=have_factor_preds)
 	yhat <- as.vector(fitted(model))
 	p <- model$rank
 	p1 <- 1:p

R/cat_durbin.R

@@ -14,12 +14,35 @@ have_factor_preds_mf <- function(mf) {
         names(xlevels) <- factnames
         attr(have_factor_preds, "xlevels") <- xlevels 
         attr(have_factor_preds, "factnames") <- factnames
+        pred_contrasts <- character(length(factnames))
+        pred_ordered <- logical(length(factnames))
+        for (pred in seq(along=factnames)) {
+            contr <- C(mf[[factnames[pred]]])
+            pred_contrasts[pred] <- attr(contr, "contrasts")
+            pred_ordered[pred] <- names(attr(contr, "contrasts")) == "ordered"
+        }
+        names(pred_contrasts) <- names(pred_ordered) <- factnames
+        attr(have_factor_preds, "pred_contrasts") <- pred_contrasts
+        attr(have_factor_preds, "pred_ordered") <- pred_ordered
     }
     have_factor_preds
 }
 
 warn_factor_preds <- function(x) {
-    warning("use of spatially lagged factors (categorical variables)\n", 
+    plural <- length(attr(x, "factnames")) > 1L
+    warning("use of spatially lagged ", ifelse(plural, "factors", "factor"),
+        " (categorical ", ifelse(plural, "variables", "variable"), "\n", 
         paste(attr(x, "factnames"), collapse=", "),
         "\nis not well-understood")
+    pred_ordered <- attr(x, "pred_ordered")
+    if (any(pred_ordered)) {
+        pred_contrasts <- attr(x, "pred_contrasts")
+        ordered <- which(pred_ordered)
+        plural <- length(ordered) > 1L
+        warning("In addition ", ifelse(plural, "variables", "variable"), ":\n",
+            paste(names(pred_ordered)[ordered], collapse=", "), 
+            "\n", ifelse(plural, "are", "is"), 
+            " ordered (ordinal) with contrasts:\n",
+            paste(pred_contrasts[ordered], collapse=", "))
+    }
 }

inst/tinytest/test_Durbin_factor.R

@@ -0,0 +1,16 @@
+library(spdep)
+data(oldcol)
+lw <- 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)
+lm_obj <- lm(f, data=COL.OLD)
+expect_warning(COL.SD0 <- SD.RStests(lm_obj, lw, test="SDM", Durbin=TRUE))
+expect_warning(COL.SD1 <- SD.RStests(lm_obj, lw, test="SDM", Durbin=~ INC + HOVAL + fDISCBD*fEW))
+expect_warning(COL.SD2 <- SD.RStests(lm_obj, lw, test="SDM", Durbin=~ INC + HOVAL + fEW))
+expect_silent(COL.SD3 <- SD.RStests(lm_obj, lw, test="SDM", Durbin=~ INC + HOVAL))
+expect_warning(COL.SDE0 <- SD.RStests(lm_obj, lw, test="SDEM", Durbin=TRUE))
+expect_warning(COL.SDE1 <- SD.RStests(lm_obj, lw, test="SDEM", Durbin=~ INC + HOVAL + fDISCBD*fEW))
+expect_warning(COL.SDE2 <- SD.RStests(lm_obj, lw, test="SDEM", Durbin=~ INC + HOVAL + fEW))
+expect_silent(COL.SDE3 <- SD.RStests(lm_obj, lw, test="SDEM", Durbin=~ INC + HOVAL))
+