将稳健的标准误差添加到修改后的 outreg 函数时出现问题

Question

我已经被困了几个小时了。我想向 Paul Johnsson 的 R outreg 函数添加强大的标准错误（和其他一些东西）（原始代码在这里 http://pj.freefaculty.org/R/outreg-worked.R）

我在使标准错误变得稳健方面遇到了很大的问题。在代码中，我已替换

se<-sqrt(diag(vcov(model)))[regname]

为

se<-sqrt(diag(vcovHC(model,type="HC1" )))[regname]

但无论如何，当我使用 Sweave 生成​​乳胶文件时，会出现非稳健标准错误。 我已经查看了，并且我的代码块确实返回了强大的 SE（或其他值，然后是原始代码）。我要生气了！ =)

这是我的函数的完整代码。

### Paul Johnson
### Adapted from ideas in post in r-help by Dave Armstrong May 8, 2006


###tight means one column per fitted model
###not tight means 2 columns per fitted model

###incoming= either one regression model or a list of regresion models
###title = a string
###modelLabels= a VECTOR of character strings
### varLabels= a LIST of labels linked to variable names (see examples)
### tight= BOOLEAN, indicates results should be on one tight column or two for each model
### showAIC= BOOLEAN should the AIC be displayed for each model?
### lyx=create a table suitable for inclusion in a lyx float.

outreg <- function(incoming, title="My Regression", label="", modelLabels=NULL, varLabels=NULL, tight=TRUE, showAIC=TRUE, lyx=TRUE){

  modelList <- NULL

  ## was input just one model, or a list of models?  ###
  if ( "lm" %in% class(incoming)) { ##just one model input
    nmodels <- 1
    modelList <- list(modl1=incoming)
  } else {
    nmodels <- length(incoming)
    modelList <- incoming
  } 

  ##TODO modelLabels MUST have same number of items as "incoming"


  ## Get a regression summary object for each fitted model
  summaryList <- list()
  fixnames <- vector()
  myModelClass <- vector()

  i <-  1
  for (model in modelList){
    summaryList[[i]] <- summary(model)

    fixnames <- unique( c( fixnames, names(coef(model))))
    myModelClass[i] <- class(model)[1]
    i <- i+1
  }



  ###If you are just using LaTeX, you need these
if (lyx == FALSE){
 cat("\\begin{table}\n ")
  cat("\\caption{",title,"}\\label{",label,"}\n ")
 }
  cat("\\begin{center}\n ")
  nColumns <- ifelse(tight, 1+nmodels, 1 + 2*nmodels) 
  cat(paste("\\begin{tabular}{*{",nColumns,"}{l}}\n ", sep=""))
  cat("\\hline\n ")

  ### Put model labels on top of each model column, if modelLabels were given
  if (!is.null(modelLabels)){
    cat("     ")
    for (modelLabel in modelLabels){
      if (tight == T) {
        cat(paste("&", modelLabel))
      }else{
        cat(paste("&\\multicolumn{2}{c}{",modelLabel,"}",sep=""))
      }
    }
    cat (" \\\\\n ")
  }

  ### Print the headers "Estimate" and "(S.E.)", output depends on tight or other format 
  if (tight == T){
    cat("             ")
    for (i in 1:nmodels) { cat (" & Estimate ") }
    cat(" \\\\\n")

    cat("             ")
    for (i in 1:nmodels) {  cat (" & (S.E.) ") }
    cat(" \\\\\n")
  }else{

    cat("             ")
    for (i in 1:nmodels) { cat (" & Estimate & S.E.") }
    cat(" \\\\\n")
  }


  cat("\\hline \n \\hline\n ")


  ### Here come the regression coefficients
  for (regname in fixnames){
    if ( !is.null(varLabels[[regname]]) ) { cat(paste("",varLabels[[regname]]), sep="")}
    else {cat(paste("", regname), sep="")}

    for (model in modelList) {
      est <- coef(model)[regname]
      se <- sqrt(diag(vcov(model)))[regname]
      if ( !is.na(est) ) {
        cat (paste("   &   ", round(est,3)))
        pval <- pt(abs(est/se), lower.tail=F, df = model$df.residual)
        if (pval < 0.025) cat("*")

        if (tight == F) {
          cat (paste("   &   (", round(se,3),")",sep=""))
        }
      } else {
        cat ("   & . ")
        if (tight == F) cat (" & " )
      }
    }
    cat (" \\\\\n ")

    if (tight == T){
      for (model in modelList) {
        est <- coef(model)[regname]
        if (!is.na(est)) cat (paste("   &   (",round(sqrt(diag(vcov(model)))[regname],3)),")",sep="")
        else cat("   &  ")
      }
      cat (" \\\\\n ")
    }
   }
  cat("\\hline \n")


  ### Print a row for the number of cases
  cat(paste("N"), sep="")
  for (model in summaryList) {
    myDF <- sum( model$df[-3] ) #omit third value from df vector
    cat (paste("   &   ", myDF))
    if (tight == F) cat("    &")
  }
  cat (" \\\\\n ")


  ### Print a row for the root mean square error
  if ("lm" %in% myModelClass) {
       cat(paste("$RMSE$"),sep="")
       for (model in summaryList) {
         cat( paste("       &", if(is.numeric(model$sigma)) round(model$sigma,3)))
         if (tight == F) cat("    &")
       }
       cat ("  \\\\\n ")
     }


  ### Print a row for the R-square
  if ("lm" %in% myModelClass) {
     cat(paste("$R^2$"),sep="")
     for (model in summaryList) {
       cat( paste("       &", if(is.numeric(model$r.square))round(model$r.square,3)))
       if (tight == F) cat("    &")
     }
     cat ("  \\\\\n ")
   }


  ## Print a row for the model residual deviance
   if ("glm" %in% myModelClass) {
    cat(paste("$Deviance$"),sep="")
    for (model in summaryList) {
      cat (paste("      &", if(is.numeric(model$deviance))round(model$deviance,3)))
      if (tight == F) cat("      &")
    }
    cat ("  \\\\\n ")
  }

  ### Print a row for the model's fit, as -2LLR
  if ("glm" %in% myModelClass) {    
    cat (paste("$-2LLR (Model \\chi^2)$"),sep="")
    for (model in modelList) {
      if (is.numeric(model$deviance)){
        n2llr <- model$null.deviance - model$deviance 
        cat (paste("      &", round(n2llr,3)))
        gmdf <- model$df.null - model$df.residual + 1

        if (pchisq(n2llr, df= gmdf, lower.tail=F) < 0.05) {cat ("*")}
      }

      else {
        cat ("    &")
      }
      if (tight == F) cat("      &")
    }
    cat ("  \\\\\n ")
  }



  ## Print a row for the model's fit, as -2 LLR
  ### Can't remember why I was multiplying by -2

  if (showAIC == T) {
    cat(paste("$AIC$"),sep="")
    for (model in modelList) {
      cat (paste("      &", if(is.numeric(AIC(model)))round(AIC(model),3)))
      if (tight == F) cat("      &")
    }
    cat ("  \\\\\n ")
  }



   cat("\\hline\\hline\n")
   cat ("* $p \\le 0.05$")
   cat("\\end{tabular}\n")
   cat("\\end{center}\n")
   if (lyx == FALSE){ 
      cat("\\end{table}\n")
   }
 }

Answer 1

您仅编辑了该代码中 vcov() 调用的一个实例。您错过了本节：

if (tight == T) {
    for (model in modelList) {
        est <- coef(model)[regname]
        if (!is.na(est)) cat (paste("   &   (",round(sqrt(diag(vcov(model))[regname],3)),")",sep="")
        else cat("   &  ")
    }
    cat (" \\\\\n ")
}

那么，您是否使用默认参数 tight = TRUE 来调用该函数？我怀疑这就是问题所在。

Answer 2

我的猜测是，Sweave 只是返回到包中的原始函数。

我的猜测是错误的。参见加文斯的回答。我认为这个问题可能是错误的根源，所以我留下这个答案供参考。但你的问题出在别处。

---

原始答案：

如果您使用 Sweave，则应该显式加载您采用的函数。更重要的是，我会给你采用的函数一个不同的名称（如 outreg2 等）并使用它。如果您不显式加载它，您很可能会收到一条错误消息，指出无法找到 outreg2。

附带说明：如果您想在 R 会话中临时编辑函数，可以使用 这个问题。