Problem with adding robust standard errors to modified outreg function

Question

I have been stuck for hours. I want to add robust standard errors (and some other stuff) to the R outreg function by Paul Johnsson (orginal code here http://pj.freefaculty.org/R/outreg-worked.R)

I have a great problem thou making the standard errors robust. In the code I have replaced

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

With

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

But anyway the nonrobust standard errors appears when I generate a latex file with Sweave. I have looking at the and my code chunk do returns robust S.E. (or other values then the orginal code anyhow). I am getting mad ! =)

Here is my complete code for the function.

### 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")
   }
 }

are you absolutely sure you are loading (sourcing) your modified, and **saved** `outreg()` in your Sweave file? I think @Joris is close with his answer, but looks more like you are getting Johnson's original and not your modified code and not realising it. Try calling it `outreg2()` in the Sewave file **and** rename your function to `outreg2()` to be sure Sweave is finding the correct function. — Gavin Simpson, Dec 07 '10 at 13:28
Yes I do source it. I have now changed the name to outreg2() but I produce the same output =/ I have made other changes to the code and they appear so probably my modified functions is used. — Marcus R, Dec 07 '10 at 13:36

score 1 · Answer 1 · edited May 23 '17 at 12:11

1

~~My guess is that Sweave just goes back to the original function in the package.~~

My guess was wrong. See Gavins answer. The problem I presumed can be a source of error, so I leave this answer for reference. But your problem lies elsewhere.

Original answer:

If you Sweave, you should explicitly load your adopted function. Even more, I'd give your adopted function a different name (like outreg2 or so) and use that. Chance is big that if you don't load it explicitly, you'll get an error saying that outreg2 could not be found.

On a side note: if you want to edit a function temporarily in your R-session, you can use one of the options discussed in this question.

edited May 23 '17 at 12:11

Community

1
1

answered Dec 07 '10 at 13:12

Joris Meys

106,551
31
221
263

There isn't an `outreg` function in sandwich, and as far as I can see all he is doing is using a function from that package, and if `outreg` weren't being found the code would fail, so some version of this function is being found and used. – Gavin Simpson Dec 07 '10 at 13:25
It doesn't seem so. I have the function pasted in to my file with source(). And the function works but the se doesn't become robust as I change the code as mentioned. – Marcus R Dec 07 '10 at 13:27
The function isn't from any package. It's used from a textfile. – Marcus R Dec 07 '10 at 13:30
@ Joris Meys. How do I make a function explicit? – Marcus R Dec 07 '10 at 13:31
@Marcus: are you **sure** your are loading the correct function? Rename it `outreg2()` and call that function in your Sweave. R should through an error if it can't find your new function. At the moment, I suspect you are loading the original `outreg()`... – Gavin Simpson Dec 07 '10 at 13:38
Thanks. I did rename it to outreg2 but still the same results. – Marcus R Dec 07 '10 at 13:44
@Marcus : with explicitly I mean save it in a .r file and do source("function.r") at the beginning of your code. You did that, so my guess was wrong. But I guess Gavin found your real problem, I missed that one. – Joris Meys Dec 07 '10 at 13:48
Yes. I did that to. I did really learn a lesson =) Thanks anyway. I love this site. – Marcus R Dec 07 '10 at 14:03

score 1 · Accepted Answer · answered Dec 07 '10 at 13:46

1

You have only edited one instance of the call to vcov() in that code. You missed this section:

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 ")
}

So, are you calling the function with the argument tight = TRUE, which is the default? I suspect this is the problem.

answered Dec 07 '10 at 13:46

Gavin Simpson

170,508
25
396
453

@Gavin: nice catch, missed that one. – Joris Meys Dec 07 '10 at 13:48
@Gavin Simpson - I love you! =) I am going to work some more on this function and eventually make it public in some way later. I have used R for just a few months so I have to learn som more first. – Marcus R Dec 07 '10 at 13:50
@Marcus: easy now, my wife may be watching ;-) – Gavin Simpson Dec 07 '10 at 13:53
@Marcus : you can show some love by upvoting the answer of Gavin and accept it as the correct one (see arrows and V-sign at the left of his answer). Guess his wife could live with that ;-) – Joris Meys Dec 07 '10 at 14:17
I have marked it as correct now =) Didn't know about that function. I can't vote as a unregistred user. Will do in a fex days. Very busy right now! – Marcus R Dec 08 '10 at 01:30

Problem with adding robust standard errors to modified outreg function

2 Answers2