在空间计量经济学模型中计算 p.values：为什么 summary() 和 texreg() 之间存在不一致？

Question

我正在估计一些包含空间自回归项 rho 和空间误差项 lambda 的空间计量经济学模型。在尝试传达我的结果时，我使用了texreg 包，它接受我正在使用的sacsarlm 模型。但是，我注意到texreg 正在打印相同的 rho 和 lambda 参数的 p 值。 Texreg 似乎正在返回在模型对象的 model@LR1$p.value 槽中找到的 p 值。

参数 rho 和 lambda 的大小不同，并且具有不同的标准误差，因此它们不应具有等效的 p 值。如果我在模型对象上调用摘要，我会得到唯一的 p 值，但无法确定这些值在模型对象中的存储位置，尽管在 str(model) 调用中遍历了每个元素。

我的问题有两个：

1. 我认为这是 texreg（和 screenreg 等）函数中的错误是正确的，还是我的解释有误？
2. 如何计算正确的 p 值或在模型对象中找到它（我正在为 texreg 编写一个新的提取函数，需要找到正确的值）？

下面是一个显示问题的最小示例：

library(spdep)
library(texreg)
set.seed(42)
W.ran <- matrix(rbinom(100*100, 1, .3),nrow=100)
X <- rnorm(100)
Y <- .2 * X + rnorm(100) + .9*(W.ran %*% X)

W.test <- mat2listw(W.ran)
model <- sacsarlm(Y~X, type = "sacmixed",
                listw=W.test, zero.policy=TRUE)
summary(model)

Call:sacsarlm(formula = Y ~ X, listw = W.test, type = "sacmixed", 
    zero.policy = TRUE)

Residuals:
      Min        1Q    Median        3Q       Max 
-2.379283 -0.750922  0.036044  0.675951  2.577148 

Type: sacmixed 
Coefficients: (asymptotic standard errors) 
                   Estimate  Std. Error z value Pr(>|z|)
(Intercept)      0.91037455  0.65700059  1.3857   0.1659
X               -0.00076362  0.10330510 -0.0074   0.9941
lag.(Intercept) -0.03193863  0.02310075 -1.3826   0.1668
lag.X            0.89764491  0.02231353 40.2287   <2e-16

Rho: -0.0028541
Asymptotic standard error: 0.0059647
    z-value: -0.47849, p-value: 0.6323
Lambda: -0.020578
Asymptotic standard error: 0.020057
    z-value: -1.026, p-value: 0.3049

LR test value: 288.74, p-value: < 2.22e-16

Log likelihood: -145.4423 for sacmixed model
ML residual variance (sigma squared): 1.0851, (sigma: 1.0417)
Number of observations: 100 
Number of parameters estimated: 7 
AIC: 304.88, (AIC for lm: 585.63)

screenreg(model)

=================================
                      Model 1    
---------------------------------
(Intercept)              0.91    
                        (0.66)   
X                       -0.00    
                        (0.10)   
lag.(Intercept)         -0.03    
                        (0.02)   
lag.X                    0.90 ***
                        (0.02)   
---------------------------------
Num. obs.              100       
Parameters               7       
AIC (Linear model)     585.63    
AIC (Spatial model)    304.88    
Log Likelihood        -145.44    
Wald test: statistic     1.05    
Wald test: p-value       0.90    
Lambda: statistic       -0.02    
Lambda: p-value          0.00    
Rho: statistic          -0.00    
Rho: p-value             0.00    
=================================
*** p < 0.001, ** p < 0.01, * p < 0.05

显然，在示例中，Rho 和 Lambda 具有不同的 p 值，两者都不为零，因此 texreg 输出存在问题。任何有关为什么会发生这种情况或在哪里获得正确 p 值的帮助都非常感谢！

Answer 1

texreg作者在这里。谢谢你抓住这个。正如我在回复here 中所述，texreg 使用extract 方法从任何（目前支持的70 多种）模型对象类型中检索相关信息。 sarlm 对象的方法的 GOF 部分似乎存在故障。

这是该方法当前的样子（截至texreg 1.36.13）：

# extension for sarlm objects (spdep package)
extract.sarlm <- function(model, include.nobs = TRUE, include.aic = TRUE, 
    include.loglik = TRUE, include.wald = TRUE, include.lambda = TRUE, 
    include.rho = TRUE, ...) {
  s <- summary(model, ...)

  names <- rownames(s$Coef)
  cf <- s$Coef[, 1]
  se <- s$Coef[, 2]
  p <- s$Coef[, ncol(s$Coef)]

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()

  if (include.nobs == TRUE) {
    n <- length(s$fitted.values)
    param <- s$parameters
    gof <- c(gof, n, param)
    gof.names <- c(gof.names, "Num.\ obs.", "Parameters")
    gof.decimal <- c(gof.decimal, FALSE, FALSE)
  }
  if (include.aic == TRUE) {
    aic <- AIC(model)
    aiclm <- s$AIC_lm.model
    gof <- c(gof, aiclm, aic)
    gof.names <- c(gof.names, "AIC (Linear model)", "AIC (Spatial model)")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.loglik == TRUE) {
    ll <- s$LL
    gof <- c(gof, ll)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.wald == TRUE) {
    waldstat <- s$Wald1$statistic
    waldp <- s$Wald1$p.value
    gof <- c(gof, waldstat, waldp)
    gof.names <- c(gof.names, "Wald test: statistic", "Wald test: p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.lambda == TRUE && !is.null(s$lambda)) {
    lambda <- s$lambda
    LRpval <- s$LR1$p.value[1]
    gof <- c(gof, lambda, LRpval)
    gof.names <- c(gof.names, "Lambda: statistic", "Lambda: p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.rho == TRUE && !is.null(s$rho)) {
    rho <- s$rho
    LRpval <- s$LR1$p.value[1]
    gof <- c(gof, rho, LRpval)
    gof.names <- c(gof.names, "Rho: statistic", "Rho: p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  tr <- createTexreg(
      coef.names = names, 
      coef = cf, 
      se = se, 
      pvalues = p, 
      gof.names = gof.names, 
      gof = gof, 
      gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("sarlm", "spdep"), 
    definition = extract.sarlm)

我认为 lambda 和 rho 部分需要一些更新是对的。 sacsarlm 函数不会将其 summary 方法打印的结果存储在任何对象中，因此您正确地指出，任何使用 str 的尝试似乎都不会显示真实的 p 值等。

因此，看看spdep 包中的print.summary.sarlm 函数在打印摘要时的实际作用是有意义的。我在source code of the package on CRAN 的文件R/summary.spsarlm.R 中找到了这个函数的代码。它看起来像这样：

print.summary.sarlm <- function(x, digits = max(5, .Options$digits - 3),
    signif.stars = FALSE, ...)
{
    cat("\nCall:", deparse(x$call), sep = "", fill=TRUE)
       if (x$type == "error") if (isTRUE(all.equal(x$lambda, x$interval[1])) ||
            isTRUE(all.equal(x$lambda, x$interval[2]))) 
            warning("lambda on interval bound - results should not be used")
       if (x$type == "lag" || x$type == "mixed")
            if (isTRUE(all.equal(x$rho, x$interval[1])) ||
            isTRUE(all.equal(x$rho, x$interval[2]))) 
            warning("rho on interval bound - results should not be used")
    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, ...)
    cat("\nType:", x$type, "\n")
    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")
    }
        if (!is.null(x$coeftitle)) {
        cat("Coefficients:", x$coeftitle, "\n")
        coefs <- x$Coef
        if (!is.null(aliased <- x$aliased) && any(x$aliased)){
        cat("    (", table(aliased)["TRUE"], 
            " not defined because of singularities)\n", sep = "")
        cn <- names(aliased)
        coefs <- matrix(NA, length(aliased), 4, dimnames = list(cn, 
                    colnames(x$Coef)))
                coefs[!aliased, ] <- x$Coef
        }
        printCoefmat(coefs, signif.stars=signif.stars, digits=digits,
        na.print="NA")
    }
#   res <- LR.sarlm(x, x$lm.model)
    res <- x$LR1
        pref <- ifelse(x$ase, "Asymptotic", "Approximate (numerical Hessian)")
    if (x$type == "error") {
        cat("\nLambda: ", format(signif(x$lambda, digits)),
            ", LR test value: ", format(signif(res$statistic,
                        digits)), ", p-value: ", format.pval(res$p.value,
                        digits), "\n", sep="")
        if (!is.null(x$lambda.se)) {
                    if (!is.null(x$adj.se)) {
                        x$lambda.se <- sqrt((x$lambda.se^2)*x$adj.se)   
                    }
            cat(pref, " standard error: ", 
                format(signif(x$lambda.se, digits)),
            ifelse(is.null(x$adj.se), "\n    z-value: ",
                               "\n    t-value: "), format(signif((x$lambda/
                x$lambda.se), digits)),
            ", p-value: ", format.pval(2*(1-pnorm(abs(x$lambda/
                x$lambda.se))), digits), "\n", sep="")
            cat("Wald statistic: ", format(signif(x$Wald1$statistic, 
            digits)), ", p-value: ", format.pval(x$Wald1$p.value, 
            digits), "\n", sep="")
        }
    } else if (x$type == "sac" || x$type == "sacmixed") {
        cat("\nRho: ", format(signif(x$rho, digits)), "\n",
                    sep="")
                if (!is.null(x$rho.se)) {
                    if (!is.null(x$adj.se)) {
                        x$rho.se <- sqrt((x$rho.se^2)*x$adj.se)   
                    }
          cat(pref, " standard error: ", 
            format(signif(x$rho.se, digits)), 
                        ifelse(is.null(x$adj.se), "\n    z-value: ",
                               "\n    t-value: "), 
            format(signif((x$rho/x$rho.se), digits)),
            ", p-value: ", format.pval(2 * (1 - pnorm(abs(x$rho/
                x$rho.se))), digits), "\n", sep="")
                }
        cat("Lambda: ", format(signif(x$lambda, digits)), "\n", sep="")
        if (!is.null(x$lambda.se)) {
                    pref <- ifelse(x$ase, "Asymptotic",
                        "Approximate (numerical Hessian)")
                    if (!is.null(x$adj.se)) {
                        x$lambda.se <- sqrt((x$lambda.se^2)*x$adj.se)   
                    }
            cat(pref, " standard error: ", 
                format(signif(x$lambda.se, digits)),
            ifelse(is.null(x$adj.se), "\n    z-value: ",
                               "\n    t-value: "), format(signif((x$lambda/
                x$lambda.se), digits)),
            ", p-value: ", format.pval(2*(1-pnorm(abs(x$lambda/
                x$lambda.se))), digits), "\n", sep="")
                }
                cat("\nLR test value: ", format(signif(res$statistic, digits)),
            ", p-value: ", format.pval(res$p.value, digits), "\n",
                    sep="")
        } else {
        cat("\nRho: ", format(signif(x$rho, digits)), 
                    ", LR test value: ", format(signif(res$statistic, digits)),
            ", p-value: ", format.pval(res$p.value, digits), "\n",
                    sep="")
                if (!is.null(x$rho.se)) {
                    if (!is.null(x$adj.se)) {
                        x$rho.se <- sqrt((x$rho.se^2)*x$adj.se)   
                    }
          cat(pref, " standard error: ", 
            format(signif(x$rho.se, digits)), 
                        ifelse(is.null(x$adj.se), "\n    z-value: ",
                               "\n    t-value: "), 
            format(signif((x$rho/x$rho.se), digits)),
            ", p-value: ", format.pval(2 * (1 - pnorm(abs(x$rho/
                x$rho.se))), digits), "\n", sep="")
                }
        if (!is.null(x$Wald1)) {
            cat("Wald statistic: ", format(signif(x$Wald1$statistic, 
            digits)), ", p-value: ", format.pval(x$Wald1$p.value, 
            digits), "\n", sep="")
        }

    }
    cat("\nLog likelihood:", logLik(x), "for", x$type, "model\n")
    cat("ML residual variance (sigma squared): ", 
        format(signif(x$s2, digits)), ", (sigma: ", 
        format(signif(sqrt(x$s2), digits)), ")\n", sep="")
        if (!is.null(x$NK)) cat("Nagelkerke pseudo-R-squared:",
            format(signif(x$NK, digits)), "\n")
    cat("Number of observations:", length(x$residuals), "\n")
    cat("Number of parameters estimated:", x$parameters, "\n")
    cat("AIC: ", format(signif(AIC(x), digits)), ", (AIC for lm: ",
        format(signif(x$AIC_lm.model, digits)), ")\n", sep="")
    if (x$type == "error") {
        if (!is.null(x$Haus)) {
            cat("Hausman test: ", format(signif(x$Haus$statistic, 
            digits)), ", df: ", format(x$Haus$parameter),
                        ", p-value: ", format.pval(x$Haus$p.value, digits),
                        "\n", sep="")
        }
        }
    if ((x$type == "lag" || x$type ==  "mixed") && x$ase) {
        cat("LM test for residual autocorrelation\n")
        cat("test value: ", format(signif(x$LMtest, digits)),
            ", p-value: ", format.pval((1 - pchisq(x$LMtest, 1)), 
            digits), "\n", sep="")
    }
        if (x$type != "error" && !is.null(x$LLCoef)) {
        cat("\nCoefficients: (log likelihood/likelihood ratio)\n")
        printCoefmat(x$LLCoef, signif.stars=signif.stars,
            digits=digits, na.print="NA")
        }
        correl <- x$correlation
        if (!is.null(correl)) {
            p <- NCOL(correl)
            if (p > 1) {
                    cat("\n", x$correltext, "\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)
}

您可以在那里看到，该函数首先区分不同的子模型（error vs. sac/sacmixed vs. else），然后决定使用哪些标准错误，然后计算 p 值即时进行，无需将它们保存在任何地方。

所以这也是我们需要在 extract 方法中做的事情，以便获得与 spdep 包中的 summary 方法相同的结果。我们还需要将它从 GOF 块移到表的系数块（参见下面的 cmets 部分进行讨论）。这是我在extract 方法中采用他们的方法的尝试：

# extension for sarlm objects (spdep package)
extract.sarlm <- function(model, include.nobs = TRUE, include.loglik = TRUE,
    include.aic = TRUE, include.lr = TRUE, include.wald = TRUE, ...) {
  s <- summary(model, ...)

  names <- rownames(s$Coef)
  cf <- s$Coef[, 1]
  se <- s$Coef[, 2]
  p <- s$Coef[, ncol(s$Coef)]

  if (model$type != "error") {  # include coefficient for autocorrelation term
    rho <- model$rho
    cf <- c(cf, rho)
    names <- c(names, "$\\rho$")
    if (!is.null(model$rho.se)) {
      if (!is.null(model$adj.se)) {
        rho.se <- sqrt((model$rho.se^2) * model$adj.se)   
      } else {
        rho.se <- model$rho.se
      }
      rho.pval <- 2 * (1 - pnorm(abs(rho / rho.se)))
      se <- c(se, rho.se)
      p <- c(p, rho.pval)
    } else {
      se <- c(se, NA)
      p <- c(p, NA)
    }
  }

  if (!is.null(model$lambda)) {
    cf <-c(cf, model$lambda)
    names <- c(names, "$\\lambda$")  
    if (!is.null(model$lambda.se)) {
      if (!is.null(model$adj.se)) {
        lambda.se <- sqrt((model$lambda.se^2) * model$adj.se)   
      } else {
        lambda.se <- model$lambda.se
      }
      lambda.pval <- 2 * (1 - pnorm(abs(model$lambda / lambda.se)))
      se <- c(se, lambda.se)
      p <- c(p, lambda.pval)
    } else {
      se <- c(se, NA)
      p <- c(p, NA)
    }
  }

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()

  if (include.nobs == TRUE) {
    n <- length(s$fitted.values)
    param <- s$parameters
    gof <- c(gof, n, param)
    gof.names <- c(gof.names, "Num.\ obs.", "Parameters")
    gof.decimal <- c(gof.decimal, FALSE, FALSE)
  }
  if (include.loglik == TRUE) {
    ll <- s$LL
    gof <- c(gof, ll)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.aic == TRUE) {
    aic <- AIC(model)
    aiclm <- s$AIC_lm.model
    gof <- c(gof, aiclm, aic)
    gof.names <- c(gof.names, "AIC (Linear model)", "AIC (Spatial model)")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.lr == TRUE && !is.null(s$LR1)) {
    gof <- c(gof, s$LR1$statistic[[1]], s$LR1$p.value[[1]])
    gof.names <- c(gof.names, "LR test: statistic", "LR test: p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.wald == TRUE && !is.null(model$Wald1)) {
    waldstat <- model$Wald1$statistic
    waldp <- model$Wald1$p.value
    gof <- c(gof, waldstat, waldp)
    gof.names <- c(gof.names, "Wald test: statistic", "Wald test: p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  tr <- createTexreg(
      coef.names = names, 
      coef = cf, 
      se = se, 
      pvalues = p, 
      gof.names = gof.names, 
      gof = gof, 
      gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("sarlm", "spdep"), 
    definition = extract.sarlm)

您可以在运行时执行此代码来更新texreg 处理这些对象的方式。如果您仍然认为我没有发现任何故障，请告诉我。如果 cmets 中没有报告，我将在下一个 texreg 版本中包含此更新的 extract 方法。

通过这些更改，调用 screenreg(model, single.row = TRUE) 会产生以下输出：

=======================================
                     Model 1           
---------------------------------------
(Intercept)             0.91 (0.66)    
X                      -0.00 (0.10)    
lag.(Intercept)        -0.03 (0.02)    
lag.X                   0.90 (0.02) ***
rho                    -0.00 (0.01)    
lambda                 -0.02 (0.02)    
---------------------------------------
Num. obs.             100              
Parameters              7              
Log Likelihood       -145.44           
AIC (Linear model)    585.63           
AIC (Spatial model)   304.88           
LR test: statistic    288.74           
LR test: p-value        0.00           
=======================================
*** p < 0.001, ** p < 0.01, * p < 0.05