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