# Load the packages for reporting HC-robust SEs library(sandwich) library(lmtest) # Print an estimated coefficient # Given the model, the order of the coefficient, and no. of decimal places printcoef=function(model,num,dig){ sprintf(paste("%.",dig,"f",sep=""),round(model$coef[num],dig)) } # Print the absolute value of an estimated coefficient # Given the model, the order of the coefficient, and no. of decimal places printabscoef=function(model,num,dig){ sprintf(paste("%.",dig,"f",sep=""),abs(round(model$coef[num],dig))) } # Print the SE / HC-robust SE of an estimated coefficient # Given the model, the order of the coefficient, and no. of decimal places printse=function(model,num,dig,HC=F){ if(HC) sprintf(paste("%.",dig,"f",sep=""),round(coeftest(model,vcovHC(model,type="HC0"))[num,"Std. Error"],dig)) else sprintf(paste("%.",dig,"f",sep=""),round(summary(model)$coef[num,"Std. Error"],dig)) } # Print the R^2 of a model # Given the model and no. of decimal places printr=function(model,dig){ sprintf(paste("%.",dig,"f",sep=""),round(summary(model)$r.squared,dig)) } # Print the adjusted R^2 of a model # Given the model and no. of decimal places printrbar=function(model,dig){ sprintf(paste("%.",dig,"f",sep=""),round(summary(model)$adj.r.squared,dig)) } # Print the t stat / HC-robust t stat of an estimated coefficient # Given the model, and the order and no. of decimal places of the coefficient printt=function(model,num,dig,HC=F,accrte=F,tdig=2){ # The default decimal place is 2; it can be altered by changing tdig. if(accrte) round(summary(model)$coef[num,"t value"],tdig) # For now, "reporting accurate t-stat" only applies to usual t-stats, not robust ones. else if(HC) round(as.numeric(printcoef(model,num,dig))/as.numeric(printse(model,num,dig,T)),tdig) else round(as.numeric(printcoef(model,num,dig))/as.numeric(printse(model,num,dig)),tdig) } # Calculate the no. of integer digits # It's optional whether to count the minus sign (if any) as a digit ninteger=function(x,countminus=F){ isPos=x>=0 x=abs(x) if(x<=1) temp=1 else temp=trunc(log10(x))+1 if(!isPos&countminus) temp=temp+1 return(temp) } # Report the result of a regression estimation # Given the regression model and a vector containing the decimal places of relevant figures reportreg=function(model,dig,suffix="hat",adj=F,HC=F,intercept=T){ if(intercept) p=0 else p=-3 num=nrow(summary(model)$coef) # Number of coefficients numlist=c(1:num) # Initialize a matrix to keep the results if(HC) result=matrix("",3,3*num+2) else result=matrix("",2,3*num+2) if(intercept) variables=c("",all.vars(formula(model))[2:num]) else variables=all.vars(formula(model))[2:(num+1)] result[1,1]=paste(all.vars(formula(model))[1],suffix,sep="") result[1,2]="=" for (n in numlist) { if(model$coefficients[n]>0) { if(n>1) result[1,3*n]="+" } else result[1,3*n]="-" result[1,3*n+2]=variables[n] nspace=max(ninteger(model$coefficients[n]),ninteger(as.numeric(printse(model,n,dig[n])))+1,ninteger(as.numeric(printse(model,n,dig[n],T)))+1) result[1,3*n+1]=paste(paste(rep(" ",nspace-ninteger(model$coefficients[n])),collapse=""),printabscoef(model,n,dig[n])," ",sep="") result[2,3*n+1]=paste(paste(rep(" ",nspace-ninteger(as.numeric(printse(model,n,dig[n])))-1),collapse=""),"(",printse(model,n,dig[n]),")",sep="") if(HC) result[3,3*n+1]=paste(paste(rep(" ",nspace-ninteger(as.numeric(printse(model,n,dig[n],T)))-1),collapse=""),"[",printse(model,n,dig[n],T),"]",sep="") } colnames(result)=rep("",3*num+2) if(HC) rownames(result)=rep("",3) else rownames(result)=rep("",2) # Remove the col/row names of the matrix, for display purpose if(intercept) result=result[,-5] if(result[1,3]=="") result=result[,-3] print(result,quote=F) cat("\tn = ",nrow(model$model),", R^2 = ",printr(model,dig[num+1]),sep="") if (adj) cat(", Rbar^2 = ",printrbar(model,dig[num+1]),sep="") }