##' "Statify" regression output
##'
##' Given a regression object out a stata-like table of results.
##'
##' @title statify regression object
##' @param x a suitable object
##' @param vcov. the variance to be used
##' @param full if full output is needed
##' @param mask coefficient to be exluded
##' @param df degrees of freedom for the p-values
##' @param ...
##' @return Output
##' @author Giuseppe M. Ragusa
##' @export
statify <- function(x,
vcov. = c("HC1", "const", "HC", "HC0", "HC2", "HC3", "HC4", "HC4m", "HC5", "HAC"),
full = TRUE, mask, df=NULL, ...) {
if(missing(mask))
mask <- 1L:length(coef(x))
scipen <- options()$scipen
type <- match.arg(vcov.)
coef <- coeftest(x, df = df)
coef <- cbind(coef, confint(x, vcov. = type))
coef[,4] <- zapsmall(coef[,4], digits = 4)
coef <- coef[mask,]
vn <- rownames(coef)
mm <- match('(Intercept)', vn)
if(!is.na(mm)) {
vn[mm] <- '_cons'
vn <- vn[c(2:length(vn),1)]
coef <- coef[c(2:length(vn),1),]
}
fcoef <- regformat(format.df(coef, numeric=FALSE), len = c(7,7,3,3,7,7),
strip.zero = c(TRUE, TRUE, FALSE, FALSE, TRUE, TRUE),
width = c(11, 11, 9, 8, 13, 12))
colnames(fcoef) <- colnames(coef)
rownames(fcoef) <- rownames(coef)
cf <- fcoef
y <- names(x$model)[1]
vn <- sapply(vn, FUN=function(u) paste(substr(' ', 1, 12-nchar(u)),
u, sep=''))
if(full) {
sumx <- summary(x)
fstat <- sumx$fstatistic
info <- matrix(0,5,1)
info[1,1] <- fstat[1L]
info[2,1] <- ifelse(length(fstat)==3,
pf(fstat[1L], fstat[2L], fstat[3L], lower.tail = FALSE),
pchisq(fstat[1L], fstat[2L], lower.tail = FALSE))
info[2,1] <- zapsmall(info[2,1], digits=4)
info[3,1] <- sumx$r.squared
info[4,1] <- sumx$adj.r.squared
info[5,1] <- sqrt(sum(residuals(sumx)^2)/x$df)
info <- matrix(apply(info, 1, format, scientific=0), 1,5, byrow=TRUE)
info <- regformat(info, len = c(4, 5, 5, 5, 6), strip.zero=rep(FALSE, 5),
width=rep(7,5))
ind <- which(as.numeric(info)>10^5)
if(length(ind))
info[1,ind] <- paste(format.df(as.numeric(info[1,ind])/10^6, numeric=FALSE, dec=1), 'e+06', sep='')
rank <- formatC((x$rank-1), 'd', width = 3, digits = 0, flag = "# ")
dfr <- format(ceiling(x$df+x$rank), scientific=(1-scipen))
dfr <- paste(strtrim(' ',6-nchar(dfr)), dfr, sep='')
obs <- format(ceiling(x$df), scientific=(1-scipen))
obs <- paste(strtrim(' ',6-nchar(dfr)), dfr, sep='')
# if (type!="const")
cat('Linear regression Number of obs = ', dfr,'\n')
cat(' F(',rank,',', dfr,') = ', info[1,1],'\n', sep='')
cat(' Prob > F =', info[1,2],'\n')
cat(' R-squared =', info[1,3],'\n')
cat(' Adj R-squared =', info[1,4],'\n')
cat(' Root MSE =', info[1,5],'\n')
# else
# cat(' Source | SS df MS Number of obs =', dfr,'\n')
# cat(' Model | ', ss, length(coef(x))-1, ms, 'Prob > F =', info[1,2],'\n')
# Residual | 495615923 72 6883554.48 R-squared = 0.2196
# -------------+------------------------------ Adj R-squared = 0.2087
# Total | 635065396 73 8699525.97 Root MSE = 2623.7
}
cat('------------------------------------------------------------------------------\n')
if(type!="const")
cat(' Robust\n')
cat(' ',y,' | Coef. Std. Err. t P>|t| [95% Conf. Interval]\n', sep = '')
cat('------------------------------------------------------------------------------\n')
for(j in mask)
##cat(vn[j],' |', cf[j,1],' ',cf[j,2],' ', cf[j,3],' ',cf[j,4],' ',cf[j,5],' ',cf[j,6],'\n', sep ='')
cat(paste(vn[j], ' |', cf[j,1], cf[j,2], cf[j,3], cf[j,4], cf[j,5], cf[j,6], sep=''),'\n')
cat('------------------------------------------------------------------------------\n')
}
##' @export
eviews.regression <-
function(x, type=c("HC1", "const", "HC", "HC0", "HC2",
"HC3", "HC4", "HC4m", "HC5", "HAC"), omit.stat = FALSE, ...)
{
require(lmtest, quietly = TRUE)
type <- match.arg(type)
attx <- attributes(x)
ff <- attx$factors
yv <- as.character(x$terms[[2]])
yv <- ifelse(length(yv)>1, yv[2], yv[1])
rs <- ""
rs <- paste0("\nDependent Variable: ", yv, "\nMethod: Least Squares\n")
rs <- paste0(rs, "Date:", format(Sys.time(), "%m-%d-%Y"), ' Time:', format(Sys.time(), "%H:%M:%S"),'\n')
rs <- paste0(rs, "Included Observations: ", length(x$fitted.values), ' after adjustments\n')
if(type!="const")
rs <- paste0(rs, 'White Heteroskedasticity-Consistent Standard Errors & Covariance\n')
rs <- paste0(rs, '====================================================================\n')
rs <- paste0(rs, ' Coefficients Std.Error t-Statistic Prob.\n')
rs <- paste0(rs, '====================================================================\n')
cout <- coef(x)
sx <- summary.reg(x, type=type)
resid <- x$resid
ssq <- sum(resid^2)
n <- length(resid)
k <- length(coef(x))
mX <- model.matrix(x)
OUTbelow <- matrix(0,14,1)
OUTabove <- coef(sx)
if(rownames(OUTabove)[1]=='(Intercept)')
rownames(OUTabove)[1]='C '
colnames(OUTabove) <- rep("",4)
colnames(OUTbelow) <- ""
OUTbelow[1,] <- sx$r.squared
OUTbelow[2,] <- sx$adj.r.squared
OUTbelow[3,] <- sqrt(ssq/(n-k))
OUTbelow[4,] <- ssq
OUTbelow[5,] <- logLik(x)
OUTbelow[6,] <- sx$fstatistic[1]
if (type!="const")
OUTbelow[7,] <- pchisq(sx$fstatistic[[1]], df=sx$fstatistic[[2]])
else
OUTbelow[7,] <- pf(sx$fstatistic[[1]], df1=sx$fstatistic[[2]],
df2=n-k, lower=FALSE)
OUTbelow[8,] <- mean(x$effects)
OUTbelow[9,] <- sd(x$effects)
OUTbelow[10,] <- AIC(x)
OUTbelow[11,] <- AIC(x, k=log(n)) ## SBIC
OUTbelow[12,] <- AIC(x, k=((k/2)*log(n) + (1/2)*log(det(crossprod(mX))))/k)
OUTbelow[13,] <- dwtest(x)$stat
OUTbelow[14,] <- NaN
aa <- formatC(OUTabove,digits=6, format = 'f', width = 12)
rr <- format(rownames(aa), width=16)
for (j in 1:nrow(aa)) {
rs <- paste0(rs, paste0(rr[j], paste0(paste0(paste0(aa[j,1], aa[j,2]), aa[j,3]), aa[j,4]), "\n"))
}
rs <- paste0(rs, '====================================================================\n')
if(!omit.stat) {
no <- c('R-squared ',
'Adjusted R-squared',
'S.E. of regression',
'Sum squared resid ',
'Log likelihood ',
'F-statistic ',
'Prob(F-statistic) ',
' Mean dependent var ',
' S.D. dependent var ',
' Akaike info criterion ',
' Schwarz criterion ',
' Hannan-Quinn criter. ',
' Durbin-Watson stat ',
' ')
for (j in 1:6)
rs <- paste0(rs, no[j],
formatC(OUTbelow[j,1], digits = 6, format = 'f',
width = 12),
no[j+7],
formatC(OUTbelow[j+7,1], digits = 6, format = 'f',
width = 12), '\n')
rs <- paste0(rs, no[7],
formatC(OUTbelow[7,1], digits = 6, format = 'f',
width = 12), '\n')
rs <- paste0(rs, '====================================================================\n\n')
}
rs
}
##' Output regression in latex form on a line
##'
##' Details
##' @title lmlatexline
##' @param x object
##' @param ... other arguments
##' @rdname lmlatexline
##' @export
lmlatexline <- function(x, ...)
UseMethod('lmlatexline')
##' @S3method lmlatexline reg
lmlatexline.reg <- function(x, se = TRUE, vcov., ..., r2 = FALSE,
dmath = FALSE, inline = FALSE,
purge.factor.name = TRUE, purge.I = TRUE,
digits = max(3, getOption("digits") - 3),
scipen = options('scipen')[[1]])
{
coeff <- coef(x)
if(missing(vcov.))
vcov. <- vcov.reg
serr <- sqrt(diag(vcov.(x, ...)))
sc <- ifelse(coeff>0,"+","-")
cf <- names(coeff)
ff <- attributes(x$terms)$factors
yv <- rownames(ff)[rowSums(ff) == 0]
tmp <- cbind(abs(coeff), serr)
##tmp <- round(tmp_original, 2)
tmp1 <- round(tmp[,1], digits = digits)
tmp2 <- round(tmp[,2], digits = digits)
digits0 <- digits
while(any(tmp2==0)) {
test <- tmp2==0
digits0 <- digits0+1
tmp2[which(test)] <- round(tmp[which(test),2], digits = digits0)
}
digits0 <- digits
while(any(tmp1==0)) {
test <- tmp1==0
digits0 <- digits0+1
tmp1[which(test)] <- round(tmp[which(test),1], digits = digits0)
}
trio <- format(cbind(tmp1,tmp2), digits=max(1, digits-1), drop0trailing=TRUE, scientific = 10000)
if(nrow(trio)==1) rownames(trio) <- cf
if(r2)
r2 <- summary(x)$r.squared
else
r2 <- FALSE
if(purge.factor.name) {
topurge <- names(which(attributes(x$terms)$dataClasses=="factor"))
for(j in topurge) {
rownames(trio) <- str_replace_all(rownames(trio), j, "")
}
}
if(purge.I) {
tmp <- str_locate(rownames(trio), "I\\(")
whereI <- which(!is.na(tmp[,1]))
for(j in whereI) {
rownames(trio)[j] <- str_sub(rownames(trio)[j], start = 3, end = str_length(rownames(trio)[j])-1)
}
}
lmlatexline.print(trio, sc, yv, cf, dmath, r2, se, inline)
}
##' @S3method lmlatexline lm
lmlatexline.lm <- lmlatexline.reg
##' @S3method lmlatexline ivreg
lmlatexline.ivreg <- function(x, se = TRUE, vcov., ...,
dmath = FALSE, inline = FALSE,
digits = max(3, getOption("digits") - 3),
scipen = options('scipen')[[1]]) {
object <- x
coeff <- coef(object)
if(missing(vcov.))
vcov. <- vcov
if (is.function(vcov.))
serr <- sqrt(diag(vcov.(object, ...)))
else if (is.matrix(vcov.) && nrow(vcov.) == ncol(vcov.) && nrow(vcov.) == length(coeff))
serr <- sqrt(diag(vcov.))
else
error("Something wrong with the definition of the variance")
sc <- ifelse(coeff>0,"+","-")
cf <- names(coeff)
ff <- attributes(object$terms$full)$factor
yv <- rownames(ff)[rowSums(ff) == 0]
trio <- format(round(cbind(abs(coeff), serr), digits), digits = digits, trim = TRUE,
scientific = scipen)
r2 <- FALSE
lmlatexline.print(trio, sc, yv, cf, dmath, r2, se, inline)
}
lmlatexline.print <- function(trio, sc, yv, cf, dmath, r2, se, inline) {
pstd <- function(z)
paste('\\underset{(',z[2],')}{',z[1],'} \\,',rownames(z), sep = '')
pstd_nose <- function(z) {
paste(z[1], '\\,', rownames(z), sep = '')
}
if(cf[[1]] == '(Intercept)')
if(se) {
rs <- paste(yv, "=", "\\underset{(", trio[1,2], ")}{", trio[1,1], "}", sep = '')
} else {
rs <- paste(yv, "=", trio[1,1], sep = '')
}
else
if(se) {
rs <- paste(yv,'=', sc[1], pstd(trio[1,,drop=FALSE]))
} else {
rs <- paste(yv,'=', sc[1], trio[1,1])
}
if(NROW(trio)>1) {
for (j in 2:NROW(trio)) {
if(se) {
rs <- paste(rs, sc[j], pstd(trio[j,,drop = FALSE]), sep = '')
} else {
rs <- paste(rs, sc[j], pstd_nose(trio[j,,drop=FALSE]), sep="")
}
}
}
if(r2) {
if(dmath)
rs <- paste(rs, '\\condition[,]{$R^2 =', format(r2, digits = 3),'$}', sep = '')
else
rs <- paste(rs, ',\\quad R^2 =', format(r2, digits = 3), sep = '')
}
if(inline)
rs <- paste0("\n$",rs,"$\n")
else {
if(dmath)
rs <- paste0("\n\\begin{dmath*}\n", rs, "\n")
else
rs <- paste0("\n\\[\n", rs, "\n")
#cat("\n",rs,"\n")
if(dmath)
rs <- paste0(rs, "\\end{dmath*}\n\n")
else
rs <- paste0(rs, "\\]\n\n")
}
rs
}
format.perc <- function(probs, digits)
## Not yet exported, maybe useful in other contexts:
## quantile.default() sometimes uses a version of it
## Manually exported!!
paste(format(100 * probs, trim = TRUE, scientific = FALSE, digits = digits), "%")
fswv <- function(x, digits = 3) {
formatC(x, format = 'f', digits = digits)
}
printMatCoef <-
function(x, digits = max(3, getOption("digits") - 2),
signif.stars = getOption("show.signif.stars"),
signif.legend = signif.stars,
dig.tst = max(1, min(5, digits - 1)),
cs.ind = 1:k, tst.ind = k+1, zap.ind = integer(0),
P.values = NULL,
has.Pvalue = nc >= 4 && substr(colnames(x)[nc], 1, 3) == "Pr(",
eps.Pvalue = .Machine$double.eps,
na.print = "NA", ...)
{
## For printing ``coefficient matrices'' as they are in summary.xxx(.) where
## xxx in {lm, glm, aov, ..}. (Note: summary.aov(.) gives a class "anova").
## By Default:
## Assume: x is a matrix-like numeric object.
## ------ with *last* column = P-values --iff-- P.values (== TRUE)
## columns {cs.ind}= numbers, such as coefficients & std.err [def.: 1L:k]
## columns {tst.ind}= test-statistics (as "z", "t", or "F") [def.: k+1]
if(is.null(d <- dim(x)) || length(d) != 2)
stop("'x' must be coefficient matrix/data frame")
nc <- d[2L]
if(is.null(P.values)) {
scp <- getOption("show.coef.Pvalues")
if(!is.logical(scp) || is.na(scp)) {
warning("option \"show.coef.Pvalues\" is invalid: assuming TRUE")
scp <- TRUE
}
P.values <- has.Pvalue && scp
}
else if(P.values && !has.Pvalue)
stop("'P.values' is TRUE, but 'has.Pvalue' is not")
if(has.Pvalue && !P.values) {# P values are there, but not wanted
d <- dim(xm <- data.matrix(x[,-nc , drop = FALSE]))
nc <- nc - 1
has.Pvalue <- FALSE
} else xm <- data.matrix(x)
k <- nc - has.Pvalue - (if(missing(tst.ind)) 1 else length(tst.ind))
if(!missing(cs.ind) && length(cs.ind) > k) stop("wrong k / cs.ind")
Cf <- array("", dim=d, dimnames = dimnames(xm))
ok <- !(ina <- is.na(xm))
## zap before deciding any formats
for (i in zap.ind) xm[, i] <- zapsmall(xm[, i], digits)
if(length(cs.ind)) {
acs <- abs(coef.se <- xm[, cs.ind, drop=FALSE])# = abs(coef. , stderr)
if(any(ia <- is.finite(acs))) {
## #{digits} BEFORE decimal point -- for min/max. value:
digmin <- 1 + if(length(acs <- acs[ia & acs != 0]))
floor(log10(range(acs[acs != 0], finite = TRUE))) else 0
Cf[,cs.ind] <- format(round(coef.se, max(1,digits-digmin)),
digits=digits)
}
}
if(length(tst.ind))
Cf[, tst.ind]<- format(round(xm[, tst.ind], digits = dig.tst),
digits = digits)
if(any(r.ind <- !((1L:nc) %in%
c(cs.ind, tst.ind, if(has.Pvalue) nc))))
for(i in which(r.ind)) Cf[, i] <- format(xm[, i], digits=digits)
okP <- if(has.Pvalue) ok[, -nc] else ok
## we need to find out where Cf is zero. We can't use as.numeric
## directly as OutDec could have been set.
## x0 <- (xm[okP]==0) != (as.numeric(Cf[okP])==0)
x1 <- Cf[okP]
dec <- getOption("OutDec")
if(dec != ".") x1 <- chartr(dec, ".", x1)
x0 <- (xm[okP] == 0) != (as.numeric(x1) == 0)
if(length(not.both.0 <- which(x0 & !is.na(x0)))) {
## not.both.0==TRUE: xm !=0, but Cf[] is: --> fix these:
Cf[okP][not.both.0] <- format(xm[okP][not.both.0],
digits = max(1, digits-1))
}
if(any(ina)) Cf[ina] <- na.print
if(P.values) {
if(!is.logical(signif.stars) || is.na(signif.stars)) {
warning("option \"show.signif.stars\" is invalid: assuming TRUE")
signif.stars <- TRUE
}
if(any(okP <- ok[,nc])) {
pv <- as.vector(xm[, nc]) # drop names
Cf[okP, nc] <- base::format.pval(pv[okP],
digits = dig.tst, eps = eps.Pvalue)
signif.stars <- signif.stars && any(pv[okP] < .1)
if(signif.stars) {
Signif <- symnum(pv, corr = FALSE, na = FALSE,
cutpoints = c(0, .001,.01,.05, .1, 1),
symbols = c("***","**","*","."," "))
Cf <- cbind(Cf, format(Signif)) #format.ch: right=TRUE
}
} else signif.stars <- FALSE
} else signif.stars <- FALSE
#print.default(Cf, quote = FALSE, right = TRUE, na.print=na.print, ...)
#if(signif.stars && signif.legend)
# cat("---\nSignif. codes: ",attr(Signif,"legend"),"\n")
invisible(Cf)
}
stripzero <- function(string) {
if(is.matrix(string))
string <- apply(str, 1, function(x) if(!is.na(pmatch('0.',x))) paste('.',strsplit(x, '0.')[[1]][2], sep =''))
if(is.array(string)| length(string)>1)
for(j in 1:length(string))
if(!is.na(pmatch('0.',string[j])))
string[j] <- paste('.',strsplit(string[j], '0.')[[1]][2], sep = '')
if(!is.na(pmatch('0.',string)))
string <- paste('.',strsplit(string, '0.')[[1]][2], sep = '')
return(string)
}
regformat <- function(x, len, strip.zero, width,...) {
## Stata Default for coef
## len <- c(8,7,4,4,7,7)
## strip.zero <- c(TRUE, TRUE, FALSE, FALSE, TRUE, TRUE)
nr <- nrow(x)
nc <- ncol(x)
xout <- matrix('', ncol=nc, nrow=nr)
for(h in 1:nc) {
for(j in 1:nr){
tmp <- strsplit(x[j,h], "\\.")[[1]]
intg <- gsub("^\\s+|\\s+$", "", tmp[1])
mant <- tmp[2]
nintg <- nchar(intg)-ifelse(strip.zero[h], (intg=="0"|intg=="-0"), 0)
nmant <- nchar(mant)
decimal <- len[h]-nintg+(strsplit(intg, '-')[[1]][1]=='')
if(decimal<=0){
tmp <- format.df(ceiling(as.numeric(x[j,h])), numeric=FALSE)
} else {
tmp <- format.df(as.numeric(x[j,h]), dec = decimal, numeric.dollar=FALSE)
}
if(strip.zero[h]) {
tmp0 <- strsplit(tmp, "\\.")[[1]]
if(tmp0[1]=="0")
tmp <- paste('.', tmp0[2], sep= '')
if(tmp0[1]=="-0")
tmp <- paste('-.', tmp0[2], sep= '')
}
xout[j,h] <- paste(substr(' ', 1, width[h]-nchar(tmp)), tmp, sep='')
}
}
##matrix(format.df(as.numeric(xout), numeric=F), ncol=nc, nrow=nr)
xout
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.