From 3faf1feee78e05b832b33fae921bfd5b773264c6 Mon Sep 17 00:00:00 2001 From: Roger Bivand Date: Mon, 18 Nov 2024 16:54:21 +0100 Subject: [PATCH] add AK test for stsls --- DESCRIPTION | 4 +- NAMESPACE | 2 +- R/kpgm_new.R | 10 +- R/s2sls.R | 763 ++++++++++++++++++++++++++------------------------ man/gstsls.Rd | 3 +- man/stsls.Rd | 15 +- 6 files changed, 412 insertions(+), 385 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1e0a96e..b14db9d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: spatialreg -Version: 1.3-5 -Date: 2024-08-19 +Version: 1.3-6 +Date: 2024-11-19 Title: Spatial Regression Analysis Encoding: UTF-8 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bivand@nhh.no", comment=c(ORCID="0000-0003-2392-6140")), diff --git a/NAMESPACE b/NAMESPACE index c8c0e0f..3f98b18 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,7 +16,7 @@ importFrom("stats", "lm.fit", "as.formula", "coef", "lm", "model.extract", importFrom("spdep", "is.symmetric.glist", "lag.listw", "card", "nb2listw", "n.comp.nb", "listw2U", "listw2mat", "Szero", "listw2sn", "mat2listw", "is.symmetric.nb", "chkIDs", "get.spChkOption", - "subset.listw") + "subset.listw", "spweights.constants") import(Matrix#, except=c("expm") ) diff --git a/R/kpgm_new.R b/R/kpgm_new.R index 2413ffb..255829b 100644 --- a/R/kpgm_new.R +++ b/R/kpgm_new.R @@ -518,9 +518,9 @@ impacts.Gmsar <- function(obj, ..., n=NULL, tr=NULL, R=NULL, listw=NULL, ####SARAR model gstsls<-function (formula, data = list(), listw, listw2=NULL, - na.action = na.fail, zero.policy = attr(listw, "zero.policy"), pars=NULL, scaleU=FALSE, - control = list(), verbose = NULL, method = "nlminb", robust = FALSE, - legacy = FALSE, W2X = TRUE ) + na.action = na.fail, zero.policy = attr(listw, "zero.policy"), pars=NULL, + scaleU=FALSE, control = list(), verbose = NULL, method = "nlminb", + robust = FALSE, legacy = FALSE, W2X = TRUE, sig2n_k=FALSE) { @@ -590,7 +590,7 @@ gstsls<-function (formula, data = list(), listw, listw2=NULL, } instr <- cbind(WX, WWX) - firststep <- tsls(y = y, yend = wy, X = x, Zinst = instr, robust = robust, legacy = legacy) + firststep <- tsls(y = y, yend = wy, X = x, Zinst = instr, robust = robust, legacy = legacy, sig2n_k=sig2n_k) ukp <- residuals(firststep) @@ -641,7 +641,7 @@ gstsls<-function (formula, data = list(), listw, listw2=NULL, colnames(xt) <- xcolnames colnames(wyt) <- c("Rho_Wy") secstep <- tsls(y = yt, yend = wyt, X = xt, Zinst = instr, - robust = robust, legacy = legacy) + robust = robust, legacy = legacy, sig2n_k=sig2n_k) rho<-secstep$coefficients[1] coef.sac<-secstep$coefficients rest.se <- sqrt(diag(secstep$var)) diff --git a/R/s2sls.R b/R/s2sls.R index a11da2b..0f7ed10 100644 --- a/R/s2sls.R +++ b/R/s2sls.R @@ -1,372 +1,391 @@ -# Copyright 2006 by Luc Anselin and Roger Bivand -# modified by Gianfranco Piras on December 11, 2009 (added the argument legacy) -# and on March 12, 2010 (added the argument W2X) -stsls <- function(formula, data = list(), listw, zero.policy=NULL, - na.action=na.fail, robust=FALSE, HC=NULL, legacy=FALSE, W2X=TRUE) { - - - if (!inherits(listw, "listw")) - stop("No neighbourhood list") - - if (is.null(zero.policy)) - zero.policy <- get("zeroPolicy", envir = .spatialregOptions) - stopifnot(is.logical(zero.policy)) - if (!inherits(formula, "formula")) formula <- as.formula(formula) - mt <- terms(formula, data = data) - mf <- lm(formula, data, na.action=na.action, method="model.frame") - na.act <- attr(mf, "na.action") - if (!is.null(na.act)) { - subset <- !(1:length(listw$neighbours) %in% na.act) - listw <- subset(listw, subset, zero.policy=zero.policy) - } - - y <- model.extract(mf, "response") - if (any(is.na(y))) stop("NAs in dependent variable") - X <- model.matrix(mt, mf) - if (any(is.na(X))) stop("NAs in independent variable") - if (robust) { - if (is.null(HC)) HC <- "HC0" - if (!any(HC %in% c("HC0", "HC1"))) - stop("HC must be one of HC0, HC1") - } -# modified to pass zero.policy Juan Tomas Sayago 100913 - Wy <- lag.listw(listw, y, zero.policy=zero.policy) - dim(Wy) <- c(nrow(X),1) - colnames(Wy) <- c("Rho") - -# WX <- lag.listw(W,X[,2:ncol(X)]) - n <- NROW(X) - m <- NCOL(X) - xcolnames <- colnames(X) - K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) - if (m > 1) { - WX <- matrix(nrow=n, ncol=(m-(K-1))) - if(W2X) WWX <- matrix(nrow = n, ncol = ncol(WX) ) - for (k in K:m) { - wx <- lag.listw(listw, X[,k], zero.policy=zero.policy) - if(W2X) wwx <- lag.listw(listw, wx, zero.policy = zero.policy) - if (any(is.na(wx))) - stop("NAs in lagged independent variable") - WX[,(k-(K-1))] <- wx - if(W2X) WWX[, (k - (K - 1))] <- wwx - } - if(W2X) inst <- cbind(WX, WWX) - else inst <- WX - } - if (K == 2 && listw$style != "W") { -# modified to meet other styles, email from Rein Halbersma - wx1 <- as.double(rep(1, n)) - wx <- lag.listw(listw, wx1, zero.policy=zero.policy) - if(W2X) wwx <- lag.listw(listw, wx, zero.policy=zero.policy) - if (m > 1) { - inst <- cbind(wx, inst) - if(W2X) inst <- cbind(wwx, inst) - } else { - inst <- matrix(wx, nrow=n, ncol=1) - if(W2X) inst <- cbind(inst, wwx) - } -# colnames(inst) <- xcolnames - - } -# if (listw$style == "W") colnames(WX) <- xcolnames[-1] - result <- tsls(y=y, yend=Wy, X=X, Zinst=inst, robust=robust, HC=HC, - legacy=legacy) - result$zero.policy <- zero.policy - result$robust <- robust - if (robust) result$HC <- HC - result$legacy <- legacy - result$listw_style <- listw$style - result$call <- match.call() - class(result) <- "Stsls" - result -} -# result <- list(coefficients=biv,var=varb,s2=s2, -# residuals=e) - -print.Stsls <- function(x, ...) { - cat("\nCall:\n") - print(x$call) - cat("\nCoefficients:\n") - print(coef(x)) - cat("\n") - invisible(x) -} - -summary.Stsls <- function(object, correlation = FALSE, ...) { - rest.se <- sqrt(diag(object$var)) -# varnames <- names(object$coefficients) - object$Coef <- cbind(object$coefficients, rest.se, - object$coefficients/rest.se, - 2*(1-pnorm(abs(object$coefficients/rest.se)))) - if (object$robust) colnames(object$Coef) <- c("Estimate", - paste(object$HC, "std. Error"), "z value", "Pr(>|z|)") - else colnames(object$Coef) <- c("Estimate", "Std. Error", - "t value", "Pr(>|t|)") - - rownames(object$Coef) <- names(object$coefficients) - if (correlation) { - object$correlation <- diag((diag(object$var)) - ^(-1/2)) %*% object$var %*% - diag((diag(object$var))^(-1/2)) - dimnames(object$correlation) <- dimnames(object$var) - } - structure(object, class=c("summary.Stsls", class(object))) -} - -print.summary.Stsls <- function(x, digits = max(5, .Options$digits - 3), - signif.stars = FALSE, ...) { - cat("\nCall:", deparse(x$call), sep = "", fill=TRUE) - cat("\nResiduals:\n") - resid <- residuals(x) - nam <- c("Min", "1Q", "Median", "3Q", "Max") - rq <- if (length(dim(resid)) == 2L) - structure(apply(t(resid), 1, quantile), dimnames = list(nam, - dimnames(resid)[[2]])) - else structure(quantile(resid), names = nam) - print(rq, digits = digits, ...) - if (x$zero.policy) { - zero.regs <- attr(x, "zero.regs") - if (!is.null(zero.regs)) - cat("Regions with no neighbours included:\n", - zero.regs, "\n") - } - cat("\nCoefficients:", x$coeftitle, "\n") - coefs <- x$Coef - printCoefmat(coefs, signif.stars=signif.stars, digits=digits, - na.print="NA") - correl <- x$correlation - cat("\n") - if (x$robust && x$legacy) cat("Asymptotic robust residual variance: ") -# if (x$legacy) cat("Asymptotic robust residual variance: ") - else cat("Residual variance (sigma squared): ") - cat(format(signif(x$sse/x$df, digits)), ", (sigma: ", - format(signif(sqrt(x$sse/x$df), digits)), ")\n", sep="") - - if (!is.null(correl)) { - p <- NCOL(correl) - if (p > 1) { - cat("\nCorrelation of Coefficients:\n") - correl <- format(round(correl, 2), nsmall = 2, - digits = digits) - correl[!lower.tri(correl)] <- "" - print(correl[-1, -p, drop = FALSE], quote = FALSE) - } - } - cat("\n") - invisible(x) - -} - -residuals.Stsls <- function(object, ...) { - if (is.null(object$na.action)) - object$residuals - else napredict(object$na.action, object$residuals) -} - -coef.Stsls <- function(object, ...) object$coefficients - -coef.summary.Stsls <- function(object, ...) object$Coef - -deviance.Stsls <- function(object, ...) object$sse - - -impacts.Stsls <- function(obj, ..., tr=NULL, R=NULL, listw=NULL, evalues=NULL, - tol=1e-6, empirical=FALSE, Q=NULL) { - if (is.null(listw) && !is.null(obj$listw_style) && - obj$listw_style != "W") - stop("Only row-standardised weights supported") - rho <- obj$coefficients[1] - beta <- obj$coefficients[-1] - icept <- grep("(Intercept)", names(beta)) - iicept <- length(icept) > 0 - if (iicept) { - P <- matrix(beta[-icept], ncol=1) - bnames <- names(beta[-icept]) - } else { - P <- matrix(beta, ncol=1) - bnames <- names(beta) - } - p <- length(beta) - n <- length(obj$residuals) - mu <- c(rho, beta) - Sigma <- obj$var - irho <- 1 - drop2beta <- 1 - res <- intImpacts(rho=rho, beta=beta, P=P, n=n, mu=mu, Sigma=Sigma, - irho=irho, drop2beta=drop2beta, bnames=bnames, interval=NULL, - type="lag", tr=tr, R=R, listw=listw, evalues=evalues, tol=tol, - empirical=empirical, Q=Q, icept=icept, iicept=iicept, p=p, - zero_fill=NULL, dvars=NULL) - attr(res, "iClass") <- class(obj) - if (!is.null(obj$robust)) { - attr(res, "robust") <- obj$robust - attr(res, "HC") <- obj$HC - } - res -} - - -# Copyright 2004 by Luc Anselin -# spatial two stage least squares -# Usage: -# stsls(listw,y,X,robust) -# Arguments: -# listw: spatial weights file as listw object -# y: dependent variable as vector -# X: explanatory variables as matrix using cbind(1,var1,...) -# robust: flag for heteroskedastic robust estimator -# Details: -# calls tsls with y as dependent variable, spatial lag of y -# as endogenous, X as exogenous variables, spatial lags of -# X as instruments and robust as specified -# Value: -# a list as returned by tsls -# coefficients: coefficient estimates -# se: (asymptotic) standard error of estimates -# t: value of asymptotic t-test statistic -# p: probability of t-test (tail, two-sided) -# var: coefficient variance matrix -# s2: residual variance (using degrees of freedom N-K) -# residuals: observed y - predicted y, to be used in diagnostics - -stsls_old <- function(W,y,X,robust=FALSE) { - Wy <- lag.listw(W,y) - dim(Wy) <- c(nrow(X),1) - colnames(Wy) <- c("Rho") - WX <- lag.listw(W,X[,2:ncol(X)]) - result <- tsls(y,Wy,X,WX,robust) - result -} - - -# Copyright 2004 by Luc Anselin -# heteroskedastic two stage least squares -# helper function, called from tsls -# Usage: -# htsls(y,Z,Q,e) -# Arguments: -# y: dependent variable as vector -# Z: matrix of endogenous and exogenous variables -# Q: matrix of instruments -# e: vector of 2SLS residuals -# Details: -# uses White consistent estimator for XOmegaX -# Value: -# a list with results -# coefficients: coefficient estimates -# se: (asymptotic) standard error of coefficients -# t: value of asymptotic t-test statistic -# p: probability of t-test (tail, two-sided) -# var: coefficient variance matrix -# s2: residual variance (using N) -# residuals: observed y - predicted y - -htsls <- function(y,Z,Q,e) { - e2 <- e^2 - oQ <- e2[,1] * Q - QoQ <- crossprod(Q,oQ) - QoQi <- solve(QoQ) - QZ <- crossprod(Q,Z) - ZQoQ <- crossprod(QZ,QoQi) - v <- ZQoQ %*% QZ - vi <- solve(v) - Qy <- crossprod(Q,y) - ZQy <- ZQoQ %*% Qy - biv <- vi %*% ZQy - yp <- Z %*% biv - e <- y - yp - biv <- biv[,1,drop=TRUE] - sse <- c(crossprod(e,e)) # / nrow(Z) - df <- nrow(Z) -# sebiv <- sqrt(diag(vi)) -# tbiv <- biv / sebiv -# pbiv <- pnorm(abs(tbiv),lower.tail=FALSE) * 2 - result <- list(coefficients=biv, -# se=sebiv,t=tbiv,p=pbiv, - var=vi,sse=sse,residuals=c(e),df=df) - result -} - - -# Copyright 2004 by Luc Anselin -# two stage least squares -# Usage: -# tsls(y,yend,X,Zinst,robust=FALSE) -# Arguments: -# y: dependent variable as vector -# yend: endogenous variables as vector or matrix (using cbind) -# X: matrix of exogenous variables, including constant -# Zinst: matrix of instruments (using cbind) -# robust: flag for heteroskedastic robust estimator -# Details: -# standard two stage least squares, using explicit two stages -# uses degrees of freedom in computation of residual variance (N-K not N) -# calls htsls when robust is TRUE -# Value: -# a list with results: -# coefficients: coefficient estimates -# se: (asymptotic) standard error of estimates -# t: value of asymptotic t-test statistic -# p: probability of t-test (tail, two-sided) -# var: coefficient variance matrix -# s2: residual variance (using degrees of freedom N-K) -# residuals: observed y - predicted y, to be used in diagnostics - -tsls <- function(y,yend,X,Zinst,robust=FALSE, HC="HC0", legacy=FALSE) { -# colnames(X) <- c("CONSTANT",colnames(X)[2:ncol(X)]) - Q <- cbind(X,Zinst) - Z <- cbind(yend,X) - df <- nrow(Z) - ncol(Z) -# QQ <- crossprod(Q,Q) - Qye <- crossprod(Q,yend) - Qr <- qr(Q) - bz <- chol2inv(Qr$qr)%*% Qye -# bz <- solve(QQ,Qye) - yendp <- Q %*% bz - Zp <- cbind(yendp,X) - Qr <- qr(Zp) -# ZpZp <- crossprod(Zp,Zp) -# ZpZpi <- solve(ZpZp) - ZpZpi <- chol2inv(Qr$qr) - Zpy <- crossprod(Zp,y) - biv <- ZpZpi %*% Zpy -# biv <- crossprod(ZpZpi,Zpy) - yp <- Z %*% biv - biv <- biv[,1,drop=TRUE] - names(biv) <- colnames(Zp) - e <- y - yp - if (robust) { - if (legacy) { - result <- htsls(y,Z,Q,e) - } else { - sse <- c(crossprod(e,e)) - if (HC == "HC0") omega <- as.numeric(e^2) - else if (HC == "HC1") - omega <- (nrow(X)/df) * as.numeric(e^2) - else stop("invalid HC choice") - ZoZ<-crossprod(Zp,(Zp*omega)) - varb<-ZpZpi%*%ZoZ%*%ZpZpi - - result <- list(coefficients=biv, - var=varb, - sse=sse, - residuals=c(e), - df=df) - - } - } else { - sse <- c(crossprod(e,e)) - s2 <- sse / df - varb <- ZpZpi * s2 -# sebiv <- sqrt(diag(varb)) -# tbiv <- biv / sebiv -# pbiv <- pnorm(abs(tbiv),lower.tail=FALSE) * 2 - result <- list(coefficients=biv, -# se=sebiv,t=tbiv,p=pbiv, - var=varb, - sse=sse, - residuals=c(e), - df=df) - } - result -} +# Copyright 2006 by Luc Anselin and Roger Bivand +# modified by Gianfranco Piras on December 11, 2009 (added the argument legacy) +# and on March 12, 2010 (added the argument W2X) +# https://github.com/r-spatial/spatialreg/issues/56 adding AK test +stsls <- function(formula, data = list(), listw, zero.policy=NULL, + na.action=na.fail, robust=FALSE, HC=NULL, legacy=FALSE, W2X=TRUE, + sig2n_k=TRUE, adjust.n=FALSE) { + + + if (!inherits(listw, "listw")) + stop("No neighbourhood list") + + if (is.null(zero.policy)) + zero.policy <- get("zeroPolicy", envir = .spatialregOptions) + stopifnot(is.logical(zero.policy)) + if (!inherits(formula, "formula")) formula <- as.formula(formula) + mt <- terms(formula, data = data) + mf <- lm(formula, data, na.action=na.action, method="model.frame") + na.act <- attr(mf, "na.action") + if (!is.null(na.act)) { + subset <- !(1:length(listw$neighbours) %in% na.act) + listw <- subset(listw, subset, zero.policy=zero.policy) + } + + y <- model.extract(mf, "response") + if (any(is.na(y))) stop("NAs in dependent variable") + X <- model.matrix(mt, mf) + if (any(is.na(X))) stop("NAs in independent variable") + if (robust) { + if (is.null(HC)) HC <- "HC0" + if (!any(HC %in% c("HC0", "HC1"))) + stop("HC must be one of HC0, HC1") + } +# modified to pass zero.policy Juan Tomas Sayago 100913 + Wy <- lag.listw(listw, y, zero.policy=zero.policy) + dim(Wy) <- c(nrow(X),1) + colnames(Wy) <- c("Rho") + +# WX <- lag.listw(W,X[,2:ncol(X)]) + n <- NROW(X) + m <- NCOL(X) + xcolnames <- colnames(X) + K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1) + if (m > 1) { + WX <- matrix(nrow=n, ncol=(m-(K-1))) + if(W2X) WWX <- matrix(nrow = n, ncol = ncol(WX) ) + for (k in K:m) { + wx <- lag.listw(listw, X[,k], zero.policy=zero.policy) + if(W2X) wwx <- lag.listw(listw, wx, zero.policy = zero.policy) + if (any(is.na(wx))) + stop("NAs in lagged independent variable") + WX[,(k-(K-1))] <- wx + if(W2X) WWX[, (k - (K - 1))] <- wwx + } + if(W2X) inst <- cbind(WX, WWX) + else inst <- WX + } + if (K == 2 && listw$style != "W") { +# modified to meet other styles, email from Rein Halbersma + wx1 <- as.double(rep(1, n)) + wx <- lag.listw(listw, wx1, zero.policy=zero.policy) + if(W2X) wwx <- lag.listw(listw, wx, zero.policy=zero.policy) + if (m > 1) { + inst <- cbind(wx, inst) + if(W2X) inst <- cbind(wwx, inst) + } else { + inst <- matrix(wx, nrow=n, ncol=1) + if(W2X) inst <- cbind(inst, wwx) + } +# colnames(inst) <- xcolnames + + } +# if (listw$style == "W") colnames(WX) <- xcolnames[-1] + result <- tsls(y=y, yend=Wy, X=X, Zinst=inst, robust=robust, HC=HC, + legacy=legacy, sig2n_k=sig2n_k) +# list(coefficients=biv, var=varb, sse=sse, residuals=c(e), df=df) + sc <- spdep::spweights.constants(listw=listw, zero.policy=zero.policy, + adjust.n=adjust.n) + u <- result$residuals + I <- c((sc$n * t(u) %*% (lag.listw(listw, u, + zero.policy=zero.policy))) / (sc$S0 * result$sse)) + Z <- attr(result, "Z") + utWZ <- t(u) %*% lag.listw(listw, Z, zero.policy=zero.policy) + A <- c(utWZ %*% (result$var/(result$sse/result$df)) %*% t(utWZ)) + s0nsq <- (sc$S0/sc$n)^2 + phisq <- c((sc$S1 + (4/(result$sse/sc$n) * A)) / (s0nsq * sc$n)) + AK <- (sc$n * I^2) / phisq + result$AK <- c(AK) + result$zero.policy <- zero.policy + result$robust <- robust + if (robust) result$HC <- HC + result$legacy <- legacy + result$listw_style <- listw$style + result$call <- match.call() + class(result) <- "Stsls" + result +} +# result <- list(coefficients=biv,var=varb,s2=s2, +# residuals=e) + +print.Stsls <- function(x, ...) { + cat("\nCall:\n") + print(x$call) + cat("\nCoefficients:\n") + print(coef(x)) + cat("\n") + invisible(x) +} + +summary.Stsls <- function(object, correlation = FALSE, ...) { + rest.se <- sqrt(diag(object$var)) +# varnames <- names(object$coefficients) + object$Coef <- cbind(object$coefficients, rest.se, + object$coefficients/rest.se, + 2*(1-pnorm(abs(object$coefficients/rest.se)))) + if (object$robust) colnames(object$Coef) <- c("Estimate", + paste(object$HC, "std. Error"), "z value", "Pr(>|z|)") + else colnames(object$Coef) <- c("Estimate", "Std. Error", + "t value", "Pr(>|t|)") + + rownames(object$Coef) <- names(object$coefficients) + if (correlation) { + object$correlation <- diag((diag(object$var)) + ^(-1/2)) %*% object$var %*% + diag((diag(object$var))^(-1/2)) + dimnames(object$correlation) <- dimnames(object$var) + } + structure(object, class=c("summary.Stsls", class(object))) +} + +print.summary.Stsls <- function(x, digits = max(5, .Options$digits - 3), + signif.stars = FALSE, ...) { + cat("\nCall:", deparse(x$call), sep = "", fill=TRUE) + cat("\nResiduals:\n") + resid <- residuals(x) + nam <- c("Min", "1Q", "Median", "3Q", "Max") + rq <- if (length(dim(resid)) == 2L) + structure(apply(t(resid), 1, quantile), dimnames = list(nam, + dimnames(resid)[[2]])) + else structure(quantile(resid), names = nam) + print(rq, digits = digits, ...) + if (x$zero.policy) { + zero.regs <- attr(x, "zero.regs") + if (!is.null(zero.regs)) + cat("Regions with no neighbours included:\n", + zero.regs, "\n") + } + cat("\nCoefficients:", x$coeftitle, "\n") + coefs <- x$Coef + printCoefmat(coefs, signif.stars=signif.stars, digits=digits, + na.print="NA") + correl <- x$correlation + cat("\n") + if (x$robust && x$legacy) cat("Asymptotic robust residual variance: ") +# if (x$legacy) cat("Asymptotic robust residual variance: ") + else cat("Residual variance (sigma squared): ") + cat(format(signif(x$sse/x$df, digits)), ", (sigma: ", + format(signif(sqrt(x$sse/x$df), digits)), ")\n", sep="") + cat("Anselin-Kelejian (1997) test for residial autocorrelation: ") + cat(format(signif(x$AK, digits)), "\n p-value: ", + format.pval(pchisq(x$AK, 1, lower.tail=FALSE), digits), sep="") + + if (!is.null(correl)) { + p <- NCOL(correl) + if (p > 1) { + cat("\nCorrelation of Coefficients:\n") + correl <- format(round(correl, 2), nsmall = 2, + digits = digits) + correl[!lower.tri(correl)] <- "" + print(correl[-1, -p, drop = FALSE], quote = FALSE) + } + } + cat("\n") + invisible(x) + +} + +residuals.Stsls <- function(object, ...) { + if (is.null(object$na.action)) + object$residuals + else napredict(object$na.action, object$residuals) +} + +coef.Stsls <- function(object, ...) object$coefficients + +coef.summary.Stsls <- function(object, ...) object$Coef + +deviance.Stsls <- function(object, ...) object$sse + + +impacts.Stsls <- function(obj, ..., tr=NULL, R=NULL, listw=NULL, evalues=NULL, + tol=1e-6, empirical=FALSE, Q=NULL) { + if (is.null(listw) && !is.null(obj$listw_style) && + obj$listw_style != "W") + stop("Only row-standardised weights supported") + rho <- obj$coefficients[1] + beta <- obj$coefficients[-1] + icept <- grep("(Intercept)", names(beta)) + iicept <- length(icept) > 0 + if (iicept) { + P <- matrix(beta[-icept], ncol=1) + bnames <- names(beta[-icept]) + } else { + P <- matrix(beta, ncol=1) + bnames <- names(beta) + } + p <- length(beta) + n <- length(obj$residuals) + mu <- c(rho, beta) + Sigma <- obj$var + irho <- 1 + drop2beta <- 1 + res <- intImpacts(rho=rho, beta=beta, P=P, n=n, mu=mu, Sigma=Sigma, + irho=irho, drop2beta=drop2beta, bnames=bnames, interval=NULL, + type="lag", tr=tr, R=R, listw=listw, evalues=evalues, tol=tol, + empirical=empirical, Q=Q, icept=icept, iicept=iicept, p=p, + zero_fill=NULL, dvars=NULL) + attr(res, "iClass") <- class(obj) + if (!is.null(obj$robust)) { + attr(res, "robust") <- obj$robust + attr(res, "HC") <- obj$HC + } + res +} + + +# Copyright 2004 by Luc Anselin +# spatial two stage least squares +# Usage: +# stsls(listw,y,X,robust) +# Arguments: +# listw: spatial weights file as listw object +# y: dependent variable as vector +# X: explanatory variables as matrix using cbind(1,var1,...) +# robust: flag for heteroskedastic robust estimator +# Details: +# calls tsls with y as dependent variable, spatial lag of y +# as endogenous, X as exogenous variables, spatial lags of +# X as instruments and robust as specified +# Value: +# a list as returned by tsls +# coefficients: coefficient estimates +# se: (asymptotic) standard error of estimates +# t: value of asymptotic t-test statistic +# p: probability of t-test (tail, two-sided) +# var: coefficient variance matrix +# s2: residual variance (using degrees of freedom N-K) +# residuals: observed y - predicted y, to be used in diagnostics + +stsls_old <- function(W,y,X,robust=FALSE) { + Wy <- lag.listw(W,y) + dim(Wy) <- c(nrow(X),1) + colnames(Wy) <- c("Rho") + WX <- lag.listw(W,X[,2:ncol(X)]) + result <- tsls(y,Wy,X,WX,robust) + result +} + + +# Copyright 2004 by Luc Anselin +# heteroskedastic two stage least squares +# helper function, called from tsls +# Usage: +# htsls(y,Z,Q,e) +# Arguments: +# y: dependent variable as vector +# Z: matrix of endogenous and exogenous variables +# Q: matrix of instruments +# e: vector of 2SLS residuals +# Details: +# uses White consistent estimator for XOmegaX +# Value: +# a list with results +# coefficients: coefficient estimates +# se: (asymptotic) standard error of coefficients +# t: value of asymptotic t-test statistic +# p: probability of t-test (tail, two-sided) +# var: coefficient variance matrix +# s2: residual variance (using N) +# residuals: observed y - predicted y + +htsls <- function(y,Z,Q,e) { + e2 <- e^2 + oQ <- e2[,1] * Q + QoQ <- crossprod(Q,oQ) + QoQi <- solve(QoQ) + QZ <- crossprod(Q,Z) + ZQoQ <- crossprod(QZ,QoQi) + v <- ZQoQ %*% QZ + vi <- solve(v) + Qy <- crossprod(Q,y) + ZQy <- ZQoQ %*% Qy + biv <- vi %*% ZQy + yp <- Z %*% biv + e <- y - yp + biv <- biv[,1,drop=TRUE] + sse <- c(crossprod(e,e)) # / nrow(Z) + df <- nrow(Z) +# sebiv <- sqrt(diag(vi)) +# tbiv <- biv / sebiv +# pbiv <- pnorm(abs(tbiv),lower.tail=FALSE) * 2 + result <- list(coefficients=biv, +# se=sebiv,t=tbiv,p=pbiv, + var=vi,sse=sse,residuals=c(e),df=df) + result +} + + +# Copyright 2004 by Luc Anselin +# two stage least squares +# Usage: +# tsls(y,yend,X,Zinst,robust=FALSE) +# Arguments: +# y: dependent variable as vector +# yend: endogenous variables as vector or matrix (using cbind) +# X: matrix of exogenous variables, including constant +# Zinst: matrix of instruments (using cbind) +# robust: flag for heteroskedastic robust estimator +# Details: +# standard two stage least squares, using explicit two stages +# uses degrees of freedom in computation of residual variance (N-K not N) +# calls htsls when robust is TRUE +# Value: +# a list with results: +# coefficients: coefficient estimates +# se: (asymptotic) standard error of estimates +# t: value of asymptotic t-test statistic +# p: probability of t-test (tail, two-sided) +# var: coefficient variance matrix +# s2: residual variance (using degrees of freedom N-K) +# residuals: observed y - predicted y, to be used in diagnostics + +tsls <- function(y,yend,X,Zinst,robust=FALSE, HC="HC0", legacy=FALSE, sig2n_k=FALSE) { +# colnames(X) <- c("CONSTANT",colnames(X)[2:ncol(X)]) + Q <- cbind(X,Zinst) + Z <- cbind(yend,X) + df <- ifelse(sig2n_k, nrow(Z) - ncol(Z), nrow(Z)) +# QQ <- crossprod(Q,Q) + Qye <- crossprod(Q,yend) + Qr <- qr(Q) + bz <- chol2inv(Qr$qr)%*% Qye +# bz <- solve(QQ,Qye) + yendp <- Q %*% bz + Zp <- cbind(yendp,X) + Qr <- qr(Zp) +# ZpZp <- crossprod(Zp,Zp) +# ZpZpi <- solve(ZpZp) + ZpZpi <- chol2inv(Qr$qr) + Zpy <- crossprod(Zp,y) + biv <- ZpZpi %*% Zpy +# biv <- crossprod(ZpZpi,Zpy) + yp <- Z %*% biv + biv <- biv[,1,drop=TRUE] + names(biv) <- colnames(Zp) + e <- y - yp + if (robust) { + if (legacy) { + result <- htsls(y,Z,Q,e) + } else { + sse <- c(crossprod(e,e)) + if (HC == "HC0") omega <- as.numeric(e^2) + else if (HC == "HC1") + omega <- (nrow(X)/df) * as.numeric(e^2) + else stop("invalid HC choice") + ZoZ<-crossprod(Zp,(Zp*omega)) + varb<-ZpZpi%*%ZoZ%*%ZpZpi + + result <- list(coefficients=biv, + var=varb, + sse=sse, + residuals=c(e), + df=df) + + } + } else { + sse <- c(crossprod(e,e)) + s2 <- sse / df + varb <- ZpZpi * s2 +# sebiv <- sqrt(diag(varb)) +# tbiv <- biv / sebiv +# pbiv <- pnorm(abs(tbiv),lower.tail=FALSE) * 2 + result <- list(coefficients=biv, +# se=sebiv,t=tbiv,p=pbiv, + var=varb, + sse=sse, + residuals=c(e), + df=df) + } + attr(result, "Z") <- Z + result +} diff --git a/man/gstsls.Rd b/man/gstsls.Rd index 3fc0093..82dd5ea 100644 --- a/man/gstsls.Rd +++ b/man/gstsls.Rd @@ -9,7 +9,7 @@ \usage{ gstsls(formula, data = list(), listw, listw2 = NULL, na.action = na.fail, zero.policy = attr(listw, "zero.policy"), pars=NULL, scaleU=FALSE, control = list(), - verbose=NULL, method="nlminb", robust=FALSE, legacy=FALSE, W2X=TRUE) + verbose=NULL, method="nlminb", robust=FALSE, legacy=FALSE, W2X=TRUE, sig2n_k=FALSE) \method{impacts}{Gmsar}(obj, \dots, n = NULL, tr = NULL, R = NULL, listw = NULL, evalues=NULL, tol = 1e-06, empirical = FALSE, Q=NULL) } @@ -34,6 +34,7 @@ neighbours, if FALSE (default) assign NA - causing \code{GMerrorsar()} to termin \item{robust}{see \code{stsls}} \item{legacy}{see \code{stsls}} \item{W2X}{see \code{stsls}} + \item{sig2n_k}{see \code{stsls}} \item{obj}{A spatial regression object created by \code{lagsarlm}, \code{lagmess} or by \code{lmSLX}; in \code{HPDinterval.LagImpact}, a LagImpact object} \item{\dots}{Arguments passed through to methods in the \pkg{coda} package} \item{tr}{A vector of traces of powers of the spatial weights matrix created using \code{trW}, for approximate impact measures; if not given, \code{listw} must be given for exact measures (for small to moderate spatial weights matrices); the traces must be for the same spatial weights as were used in fitting the spatial regression, and must be row-standardised} diff --git a/man/stsls.Rd b/man/stsls.Rd index 7236a4c..29bf7d7 100644 --- a/man/stsls.Rd +++ b/man/stsls.Rd @@ -14,7 +14,8 @@ } \usage{ stsls(formula, data = list(), listw, zero.policy = NULL, - na.action = na.fail, robust = FALSE, HC=NULL, legacy=FALSE, W2X = TRUE) + na.action = na.fail, robust = FALSE, HC=NULL, legacy=FALSE, W2X = TRUE, + sig2n_k=TRUE, adjust.n=FALSE) \method{impacts}{Stsls}(obj, \dots, tr, R = NULL, listw = NULL, evalues=NULL, tol = 1e-06, empirical = FALSE, Q=NULL) } @@ -31,8 +32,10 @@ neighbours, if FALSE (default) assign NA - causing \code{lagsarlm()} to terminat \item{na.action}{a function (default \code{na.fail}), can also be \code{na.omit} or \code{na.exclude} with consequences for residuals and fitted values - in these cases the weights list will be subsetted to remove NAs in the data. It may be necessary to set zero.policy to TRUE because this subsetting may create no-neighbour observations. Note that only weights lists created without using the glist argument to \code{nb2listw} may be subsetted.} \item{robust}{default FALSE, if TRUE, apply a heteroskedasticity correction to the coefficients covariances} \item{HC}{default NULL, if \code{robust} is TRUE, assigned \dQuote{HC0}, may take values \dQuote{HC0} or \dQuote{HC1} for White estimates or MacKinnon-White estimates respectively} - \item{legacy}{the argument chooses between two implementations of the robustness correction: default FALSE - use the estimate of Omega only in the White consistent estimator of the variance-covariance matrix, if TRUE, use the original implementation which runs a GLS using the estimate of Omega, and yields different coefficient estimates as well - see example below} - \item{W2X}{default TRUE, if FALSE only WX are used as instruments in the spatial two stage least squares; until release 0.4-60, only WX were used - see example below } + \item{legacy}{the argument chooses between two implementations of the robustness correction: default FALSE - use the estimate of Omega only in the White consistent estimator of the variance-covariance matrix, if TRUE, use the original implementation which runs a GLS using the estimate of Omega, overrides sig2n_k, and yields different coefficient estimates as well - see example below} + \item{W2X}{default TRUE, if FALSE only WX are used as instruments in the spatial two stage least squares; until release 0.4-60, only WX were used - see example below; Python \code{spreg::GM_Lag} is default FALSE} + \item{sig2n_k}{default TRUE - use n-k to calculate sigma^2, if FALSE use n; Python \code{spreg::GM_Lag} is default FALSE} + \item{adjust.n}{default FALSE, used in creating spatial weights constants for the Anselin-Kelejian (1997) test} \item{obj}{A spatial regression object created by \code{lagsarlm}, \code{lagmess} or by \code{lmSLX}; in \code{HPDinterval.LagImpact}, a LagImpact object} \item{\dots}{Arguments passed through to methods in the \pkg{coda} package} \item{tr}{A vector of traces of powers of the spatial weights matrix created using \code{trW}, for approximate impact measures; if not given, \code{listw} must be given for exact measures (for small to moderate spatial weights matrices); the traces must be for the same spatial weights as were used in fitting the spatial regression, and must be row-standardised} @@ -48,6 +51,8 @@ neighbours, if FALSE (default) assign NA - causing \code{lagsarlm()} to terminat \deqn{y = \rho W y + X \beta + \varepsilon}{y = rho W y + X beta + e} by using spatially lagged X variables as instruments for the spatially lagged dependent variable. + +From version 1.3-6, the general Anselin-Kelejian (1997) test for residual spatial autocorrelation is added. } \value{ an object of class "Stsls" containing: @@ -60,7 +65,9 @@ by using spatially lagged X variables as instruments for the spatially lagged de \references{Kelejian, H.H. and I.R. Prucha (1998). A generalized spatial two stage least squares procedure for estimating a spatial autoregressive model with autoregressive disturbances. \emph{Journal of Real Estate -Finance and Economics} 17, 99-121. +Finance and Economics} 17, 99-121. \doi{10.1023/A:1007707430416}. + +Anselin, L., & Kelejian, H. H. (1997). Testing for Spatial Error Autocorrelation in the Presence of Endogenous Regressors. International Regional Science Review, 20(1-2), 153-182. \doi{10.1177/016001769702000109}. Roger Bivand, Gianfranco Piras (2015). Comparing Implementations of Estimation Methods for Spatial Econometrics. \emph{Journal of Statistical Software}, 63(18), 1-36. \doi{10.18637/jss.v063.i18}. }