Merge pull request #49 from r-spatial/slx_ml

Slx ml

R/ML_models.R

@@ -13,7 +13,7 @@ errorsarlm <- function(formula, data = list(), listw, na.action, weights=NULL,
             super=NULL, spamPivot="MMD", in_coef=0.1, type="MC",
             correct=TRUE, trunc=TRUE, SE_method="LU", nrho=200,
             interpn=2000, small_asy=TRUE, small=1500, SElndet=NULL,
-            LU_order=FALSE, pre_eig=NULL, glht=FALSE)
+            LU_order=FALSE, pre_eig=NULL, return_impacts=TRUE)
         nmsC <- names(con)
         con[(namc <- names(control))] <- control
         if (length(noNms <- namc[!namc %in% nmsC])) 
@@ -29,7 +29,7 @@ errorsarlm <- function(formula, data = list(), listw, na.action, weights=NULL,
 #        stopifnot(is.logical(con$super))
         stopifnot(is.logical(con$compiled_sse))
         stopifnot(is.character(con$spamPivot))
-        stopifnot(is.logical(con$glht))
+        stopifnot(is.logical(con$return_impacts))
         if (!inherits(formula, "formula")) formula <- as.formula(formula)
 #	mt <- terms(formula, data = data)
 #	mf <- lm(formula, data, na.action=na.action, method="model.frame")
@@ -42,7 +42,7 @@ errorsarlm <- function(formula, data = list(), listw, na.action, weights=NULL,
         mf <- eval(mf, parent.frame())
         mt <- attr(mf, "terms")
         if (attr(mt, "intercept") == 1 && !any(attr(mt, "factors") == 1) &&
-            (is.formula(Durbin) || isTRUE(Durbin))) {
+            (!missing(Durbin)) && (is.formula(Durbin) || isTRUE(Durbin))) {
             warning("intercept-only model, Durbin invalid and set FALSE")
             Durbin <- FALSE
         }
@@ -270,7 +270,8 @@ errorsarlm <- function(formula, data = list(), listw, na.action, weights=NULL,
 	names(coef.lambda) <- xcolnames
         sum_lm_target <- summary.lm(lm.target, correlation = FALSE)
         emixedImps <- NULL
-	if (etype == "emixed") {
+        if (any(sum_lm_target$aliased)) warning("aliased variables found")
+	if (con$return_impacts && etype == "emixed") {
           if (isTRUE(Durbin)) {
             m.1 <- m > 1
             if (m.1 && K == 2) {

R/SLX_WX.R

@@ -1,6 +1,6 @@
 
 
-lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin=TRUE, zero.policy=NULL) {
+lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin=TRUE, zero.policy=NULL, return_impacts=TRUE) {
         if (is.null(zero.policy))
             zero.policy <- get("zeroPolicy", envir = .spatialregOptions)
         stopifnot(is.logical(zero.policy))
@@ -108,17 +108,19 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
             lm.model <- lm(formula(paste("y ~ 0 + ", paste(colnames(x), collapse="+"))), data=as.data.frame(x), weights=weights)
         }
         sum_lm_model <- summary.lm(lm.model, correlation = FALSE)
+        if (any(sum_lm_model$aliased)) warning("aliased variables found")
         mixedImps <- NULL
-	K <- ifelse(isTRUE(grep("Intercept",
+        if (return_impacts) {
+	  K <- ifelse(isTRUE(grep("Intercept",
             names(coefficients(lm.model))[1]) == 1L), 2, 1)
-        if (isTRUE(Durbin)) {
-          m <- length(coefficients(lm.model))
-          m.1 <- m > 1
-          if (m.1 && K == 2) { #TR: without intercept and m.1 use m/2
+          if (isTRUE(Durbin)) {
+            m <- length(coefficients(lm.model))
+            m.1 <- m > 1
+            if (m.1 && K == 2) { #TR: without intercept and m.1 use m/2
               m2 <- (m-1)/2
-          } else {
+            } else {
               m2 <- m/2
-          }
+            }
             cm <- matrix(0, ncol=m, nrow=m2)
             if (K == 2) {
                 if (m.1) {
@@ -143,7 +145,7 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
             }
             suppressWarnings(lc <- summary(multcomp::glht(lm.model, linfct=cm)))
             totImps <- cbind("Estimate"=lc$test$coefficients, "Std. Error"=lc$test$sigma)
-      } else if (is.formula(Durbin)) {
+        } else if (is.formula(Durbin)) {
 #FIXME
             LI <- ifelse(listw$style != "W" 
                          && attr(terms(Durbin), "intercept") == 1, 1, 0) #TR: lagged intercept if not W and in Durbin formula
@@ -195,10 +197,11 @@ lmSLX <- function(formula, data = list(), listw, na.action, weights=NULL, Durbin
                     }
                   }
                 rownames(totImps) <- xn
-        } else stop("undefined Durbin state")
-        mixedImps <- list(dirImps=dirImps, indirImps=indirImps,
+          } else stop("undefined Durbin state")
+          mixedImps <- list(dirImps=dirImps, indirImps=indirImps,
             totImps=totImps)
 
+        }
         attr(lm.model, "mixedImps") <- mixedImps
         attr(lm.model, "dvars") <- dvars
         if (is.formula(Durbin)) attr(lm.model, "Durbin") <- deparse(Durbin)

R/sarlm_functions.R

@@ -297,15 +297,15 @@ getVmate <- function(coefs, env, s2, trs, tol.solve=1.0e-10, optim=FALSE,
     if (optim) {
       if (optimM == "nlm") {
            options(warn=-1)
-           opt <- nlm(f=f_laglm_hess_nlm, p=coefs, env=env, hessian=TRUE)
+           opt <- nlm(f=f_errlm_hess_nlm, p=coefs, env=env, hessian=TRUE)
            options(warn=0)
            mat <- opt$hessian
-#        opt <- optimHess(par=coefs, fn=f_laglm_hess, env=env)
+#        opt <- optimHess(par=coefs, fn=f_errlm_hess, env=env)
 #        mat <- opt
        } else if (optimM == "optimHess") {
-           mat <- optimHess(par=coefs, fn=f_laglm_hess, env=env)
+           mat <- optimHess(par=coefs, fn=f_errlm_hess, env=env)
        } else {
-           opt <- optim(par=coefs, fn=f_laglm_hess, env=env, method=optimM,
+           opt <- optim(par=coefs, fn=f_errlm_hess, env=env, method=optimM,
            hessian=TRUE)
            mat <- opt$hessian
       }
@@ -354,6 +354,11 @@ f_errlm_hess <- function(coefs, env) {
     ret
 }
 
+f_errlm_hess_nlm <- function(coefs, env) {
+    ret <- f_errlm_hess(coefs, env)
+    -ret
+}
+
 insert_asye <- function(coefs, env, s2, mat, trs) {
     lambda <- coefs[1]
     p <- length(coefs)-1L

man/ML_models.Rd

@@ -147,6 +147,7 @@ Because numerical optimisation is used to find the values of lambda and rho in \
   \item{SElndet}{default NULL, may be used to pass a pre-computed SE toolbox style matrix of coefficients and their lndet values to the "SE_classic" and "SE_whichMin" methods}
   \item{LU_order}{default FALSE; used in \dQuote{LU_prepermutate}, note warnings given for \code{lu} method}
   \item{pre_eig}{default NULL; may be used to pass a pre-computed vector of eigenvalues}
+  \item{return_impacts}{default TRUE; may be set FALSE to avoid problems calculating impacts with aliased variables}
   \item{OrdVsign}{default 1; used to set the sign of the final component to negative if -1 (alpha times ((sigma squared) squared) in Ord (1975) equation B.1).}
   \item{opt_method:}{default \dQuote{nlminb}, may be set to \dQuote{L-BFGS-B} to use box-constrained optimisation in \code{optim}}
   \item{opt_control:}{default \code{list()}, a control list to pass to \code{nlminb} or \code{optim}}

man/SLX.Rd

@@ -16,7 +16,7 @@
 }
 \usage{
 lmSLX(formula, data = list(), listw, na.action, weights=NULL, Durbin=TRUE,
- zero.policy=NULL)
+ zero.policy=NULL, return_impacts=TRUE)
 \method{print}{SlX}(x, digits = max(3L, getOption("digits") - 3L), ...)
 \method{summary}{SlX}(object, correlation = FALSE, symbolic.cor = FALSE, ...)
 \method{print}{summary.SlX}(x, digits = max(3L, getOption("digits") - 3L),
@@ -39,6 +39,7 @@ is called.}
   \item{weights}{an optional vector of weights to be used in the fitting process. Non-NULL weights can be used to indicate that different observations have different variances (with the values in weights being inversely proportional to the variances); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations (including the case that there are w_i observations equal to y_i and the data have been summarized) - \code{\link{lm}}}
   \item{Durbin}{default TRUE for \code{lmSLX} (Durbin model including WX); if TRUE, full spatial Durbin model; if a formula object, the subset of explanatory variables to lag}
   \item{zero.policy}{default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA}
+  \item{return_impacts}{default TRUE; may be set FALSE to avoid problems calculating impacts with aliased variables}
   \item{digits}{the number of significant digits to use when printing}
   \item{correlation}{logical; if \code{TRUE}, the correlation matrix of the estimated parameters is returned and printed}
   \item{symbolic.cor}{logical. If \code{TRUE}, print the correlations in a symbolic form (see 'symnum') rather than as numbers}