|
| 1 | +# Copyright 2023-4 by Roger Bivand |
| 2 | +# |
| 3 | + |
| 4 | +is.formula <- function(x){ |
| 5 | + inherits(x,"formula") |
| 6 | +} |
| 7 | + |
| 8 | +create_X0 <- function(X, listw, Durbin=TRUE, data=NULL, na.act=NULL) { |
| 9 | + if (isTRUE(Durbin)) { |
| 10 | + n <- NROW(X) |
| 11 | + m <- NCOL(X) |
| 12 | + # check if there are enough regressors |
| 13 | + xcolnames <- colnames(X) |
| 14 | + stopifnot(!is.null(xcolnames)) |
| 15 | + K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) |
| 16 | + vars <- NULL |
| 17 | + xI <- NULL |
| 18 | + X0 <- NULL |
| 19 | + if (K == 2) { |
| 20 | + # unnormalized weight matrices |
| 21 | + if (!(listw$style == "W")) { |
| 22 | + xI <- as.double(rep(1, n)) |
| 23 | + vars <-"X0.(Intercept)" |
| 24 | + } |
| 25 | + } |
| 26 | + if (m > 1 || (m == 1 && K == 1)) { |
| 27 | + X0 <- matrix(as.numeric(NA), nrow=n, |
| 28 | + ncol=ifelse(m==1, 1, (m-(K-1)))) |
| 29 | + for (k in K:m) { |
| 30 | + j <- ifelse(k==1, 1, k-(K-1)) |
| 31 | + X0[,j] <- X[,xcolnames[k]] |
| 32 | + vars <- c(vars, xcolnames[k]) |
| 33 | + } |
| 34 | + } |
| 35 | + if (!is.null(xI)) X0 <- cbind(xI, X0) |
| 36 | + colnames(X0) <- vars |
| 37 | + rownames(X0) <- rownames(X) |
| 38 | + } else if (is.formula(Durbin)) { |
| 39 | + data1 <- data |
| 40 | + if (!is.null(na.act) && (inherits(na.act, "omit") || |
| 41 | + inherits(na.act, "exclude"))) { |
| 42 | + data1 <- data1[-c(na.act),] |
| 43 | + } |
| 44 | + dmf <- lm(Durbin, data1, na.action=na.fail, |
| 45 | + method="model.frame") |
| 46 | +# dmf <- lm(Durbin, data, na.action=na.action, |
| 47 | +# method="model.frame") |
| 48 | + X0 <- try(model.matrix(Durbin, dmf), silent=TRUE) |
| 49 | + if (inherits(X0, "try-error")) |
| 50 | + stop("Durbin variable mis-match") |
| 51 | + |
| 52 | + inds <- match(colnames(X0), colnames(X)) |
| 53 | + if (anyNA(inds)) { |
| 54 | + wna <- which(is.na(inds)) #TR: continue if Durbin has intercept, but formula has not |
| 55 | + if (length(wna) == 1 && grepl("Intercept", colnames(X0)[wna]) |
| 56 | + && attr(terms(Durbin), "intercept") == 1) { |
| 57 | + inds <- inds[-wna] |
| 58 | + } else { |
| 59 | + stop("X0 variables not in X: ", |
| 60 | + paste(colnames(X0)[is.na(inds)], collapse=" ")) |
| 61 | + } |
| 62 | + } |
| 63 | + icept <- grep("(Intercept)", colnames(X0)) |
| 64 | + if (length(icept) == 1L && listw$style == "W") |
| 65 | + X0 <- X0[, -icept, drop=FALSE] |
| 66 | + } else stop("Durbin argument neither TRUE nor formula") |
| 67 | + X0 |
| 68 | +} |
| 69 | + |
| 70 | +SD.RStests <- function(model, listw, zero.policy=attr(listw, "zero.policy"), test="SDM", Durbin=TRUE) { |
| 71 | + |
| 72 | + if (inherits(model, "lm")) na.act <- model$na.action |
| 73 | + else na.act <- attr(model, "na.action") |
| 74 | + |
| 75 | + listw_name <- deparse(substitute(listw)) |
| 76 | + |
| 77 | + SDM.tests <- c("SDM_RSlag", "SDM_adjRSlag", "SDM_RSWX", "SDM_adjRSWX", "SDM_Joint") |
| 78 | + SDEM.tests <- c("SDEM_RSerr", "SDEM_RSWX", "SDEM_Joint") |
| 79 | + all.tests <- c(SDM.tests, SDEM.tests) |
| 80 | + if (test[1] == "SDM") test <- SDM.tests |
| 81 | + if (test[1] == "SDEM") test <- SDEM.tests |
| 82 | + if (test[1] == "all") test <- all.tests |
| 83 | + if (!all(test %in% all.tests)) |
| 84 | + stop("Invalid test selected - must be either \"all\", \"SDM\", \"SDEM\" or a vector of tests") |
| 85 | + nt <- length(test) |
| 86 | + if (nt < 1) stop("non-positive number of tests") |
| 87 | + |
| 88 | + if (!inherits(listw, "listw")) stop(paste(listw_name, |
| 89 | + "is not a listw object")) |
| 90 | + if (is.null(zero.policy)) |
| 91 | + zero.policy <- get("zeroPolicy", envir = .spdepOptions) |
| 92 | + stopifnot(is.logical(zero.policy)) |
| 93 | + if (!is.null(na.act)) { |
| 94 | + subset <- !(1:length(listw$neighbours) %in% na.act) |
| 95 | + listw <- subset(listw, subset, zero.policy=zero.policy) |
| 96 | + } |
| 97 | + |
| 98 | + if(!inherits(model, "lm")) stop(paste(deparse(substitute(model)), |
| 99 | + "not an lm object")) |
| 100 | + N <- length(listw$neighbours) |
| 101 | + u <- resid(model) |
| 102 | + if (N != length(u)) stop("objects of different length") |
| 103 | + u <- as.vector(u) |
| 104 | + |
| 105 | + if (is.null(attr(listw$weights, "W")) || !attr(listw$weights, "W")) |
| 106 | + warning("Spatial weights matrix not row standardized") |
| 107 | + |
| 108 | + if (is.formula(Durbin)) { |
| 109 | + dt <- try(eval(model$call[["data"]]), silent=TRUE) |
| 110 | + if (inherits(dt, "try-error") || !is.data.frame(dt)) |
| 111 | + stop("data object used to fit linear model not available for formula Durbin") |
| 112 | + } |
| 113 | + |
| 114 | + y <- model.response(model.frame(model)) |
| 115 | + X <- model.matrix(terms(model), model.frame(model)) |
| 116 | + X0 <- create_X0(X=X, listw=listw, Durbin=Durbin, data=dt, na.act=na.act) |
| 117 | + yhat <- as.vector(fitted(model)) |
| 118 | + p <- model$rank |
| 119 | + p1 <- 1:p |
| 120 | + nacoefs <- which(is.na(coefficients(model))) |
| 121 | +# fixed after looking at TOWN dummy in Boston data |
| 122 | + if (length(nacoefs) > 0L) X <- X[,-nacoefs] |
| 123 | + XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE]) |
| 124 | + sigma2 <- c(t(u) %*% u) / N |
| 125 | + TrW <- tracew(listw) |
| 126 | + Wu <- lag.listw(listw, u, zero.policy) |
| 127 | + Wy <- lag.listw(listw, y, zero.policy) |
| 128 | + dr <- (t(Wy) %*% u)/sigma2 # lagged y |
| 129 | + dl <- (t(Wu) %*% u)/sigma2 # lagged residuals |
| 130 | + Wyhat <- lag.listw(listw, yhat, zero.policy) |
| 131 | + WX0 <- lag.listw(listw, X0, zero.policy) |
| 132 | + dg <- c(t(WX0) %*% u)/sigma2 |
| 133 | + k <- ncol(X) |
| 134 | + k0 <- ncol(X0) |
| 135 | + J_11 <- rbind(cbind((crossprod(X)/(N*sigma2)), rep(0, k)), |
| 136 | + cbind(t(rep(0, k)), (1/(2*(sigma2^2))))) |
| 137 | + invJ_11 <- solve(J_11) |
| 138 | + Jrp <- rbind((t(X) %*% Wyhat)/(N*sigma2), t(rep(0, 1))) |
| 139 | + Jgb <- (t(X) %*% WX0)/(N*sigma2) |
| 140 | + Jgp <- rbind(Jgb, t(rep(0, k0))) |
| 141 | + J_12 <- cbind(Jrp, Jgp) |
| 142 | + Jrr <- (c(crossprod(Wyhat)) + TrW*sigma2)/(N*sigma2) |
| 143 | + Jgg <- crossprod(WX0)/(N*sigma2) |
| 144 | + Jrg <- (t(WX0) %*% Wyhat)/(N*sigma2) |
| 145 | + J_22 <- rbind(cbind(Jrr, t(Jrg)), cbind(Jrg, Jgg)) |
| 146 | + Jrg.p <- t(Jrg) - c(t(Jrp) %*% invJ_11 %*% Jgp) |
| 147 | + Jr.p <- Jrr - c(t(Jrp) %*% invJ_11 %*% Jrp) |
| 148 | + Jg.p <- Jgg - (t(Jgp) %*% invJ_11 %*% Jgp) |
| 149 | + invJg.p <- solve(Jg.p) |
| 150 | + dr_adj <- dr - (Jrg.p %*% invJg.p %*% dg) |
| 151 | + Jr.p_adj <- Jr.p - (Jrg.p %*% invJg.p %*% t(Jrg.p)) |
| 152 | + dg_adj <- dg - c(dr * (1/Jr.p)) * Jrg.p |
| 153 | + Jg.p_adj <- Jg.p - ((1/Jr.p) * crossprod(Jrg.p)) |
| 154 | + J.22 <- solve(J_22 - t(J_12) %*% invJ_11 %*% J_12) |
| 155 | + invJg.b <- solve(Jgg - t(Jgb) %*% solve(crossprod(X)/(N*sigma2)) %*% |
| 156 | + Jgb) |
| 157 | + tres <- vector(mode="list", length=nt) |
| 158 | + names(tres) <- test |
| 159 | + for (i in 1:nt) { |
| 160 | + testi <- test[i] |
| 161 | + zz <- switch(testi, |
| 162 | + SDM_RSlag = vec <- c((1/N) * ((dr^2) * 1/Jr.p), 1), |
| 163 | + SDM_adjRSlag = vec <- c((1/N)*((dr_adj^2)*(1/Jr.p_adj)), 1), |
| 164 | + SDM_RSWX = vec <- c((1/N) * (t(dg) %*% invJg.p %*% dg), |
| 165 | + ncol(X0)), |
| 166 | + SDM_adjRSWX = vec <- c((1/N) * (dg_adj %*% solve(Jg.p_adj) %*% |
| 167 | + t(dg_adj)), ncol(X0)), |
| 168 | + SDM_Joint = vec <- c(((1/N) * (t(c(dr, dg)) %*% |
| 169 | + J.22 %*% c(dr, dg))), ncol(X0)+1), |
| 170 | + SDEM_RSerr = vec <- c((dl^2) / TrW, 1), |
| 171 | + SDEM_RSWX = vec <- c(((t(dg) %*% invJg.b %*% dg) / N), |
| 172 | + ncol(X0)), |
| 173 | + SDEM_Joint = vec <- c(((t(dg) %*% invJg.b %*% dg) / N) + |
| 174 | + ((dl^2) / TrW), ncol(X0)+1) |
| 175 | + ) |
| 176 | + if (is.null(zz)) stop(paste(testi, ": no such test", sep="")) |
| 177 | + statistic <- vec[1] |
| 178 | + names(statistic) <- testi |
| 179 | + parameter <- vec[2] |
| 180 | + names(parameter) <- "df" |
| 181 | + p.value <- 1 - pchisq(statistic, parameter) |
| 182 | + if (!is.finite(p.value) || p.value < 0 || p.value > 1) |
| 183 | + warning("Out-of-range p-value: reconsider test arguments") |
| 184 | + names(p.value) <- "" |
| 185 | + method <- "Rao's score test spatial Durbin diagnostics" |
| 186 | + Durf <- "" |
| 187 | + if (is.formula(Durbin)) |
| 188 | + Durf <- paste0("Durbin: ", paste(as.character(Durbin), |
| 189 | + collapse=" "), "\n") |
| 190 | + data.name <- paste("\n", paste(strwrap(paste("model: ", |
| 191 | + gsub("[ ]+", " ", paste(deparse(model$call), |
| 192 | + sep="", collapse="")))), collapse="\n"), |
| 193 | + "\nweights: ", listw_name, "\n", Durf, sep="") |
| 194 | + tres[[i]] <- list(statistic=statistic, parameter=parameter, |
| 195 | + p.value=p.value, method=method, data.name=data.name) |
| 196 | + class(tres[[i]]) <- "htest" |
| 197 | + } |
| 198 | + class(tres) <- "RStestlist" |
| 199 | + tres |
| 200 | +} |
| 201 | + |
| 202 | + |
0 commit comments