将最大似然估计的系数放入观星表中

Stargazer 为 lm (和其他)对象生成非常漂亮的乳胶表。假设我用最大似然法拟合了一个模型。我想让天文学家为我的估算，制作一个类似 lm 的表格。我怎么能这么做？

虽然这有点古怪，但是有一种方法可能是创建一个包含我的估计值的“假”lm 对象——我认为只要总结(my.fake.lm.object)有效，这种方法就可以工作。这容易做到吗？

举个例子:

library(stargazer)


N <- 200
df <- data.frame(x=runif(N, 0, 50))
df$y <- 10 + 2 * df$x + 4 * rt(N, 4)  # True params
plot(df$x, df$y)


model1 <- lm(y ~ x, data=df)
stargazer(model1, title="A Model")  # I'd like to produce a similar table for the model below


ll <- function(params) {
## Log likelihood for y ~ x + student's t errors
params <- as.list(params)
return(sum(dt((df$y - params$const - params$beta*df$x) / params$scale, df=params$degrees.freedom, log=TRUE) -
log(params$scale)))
}


model2 <- optim(par=c(const=5, beta=1, scale=3, degrees.freedom=5), lower=c(-Inf, -Inf, 0.1, 0.1),
fn=ll, method="L-BFGS-B", control=list(fnscale=-1), hessian=TRUE)
model2.coefs <- data.frame(coefficient=names(model2$par), value=as.numeric(model2$par),
se=as.numeric(sqrt(diag(solve(-model2$hessian)))))


stargazer(model2.coefs, title="Another Model", summary=FALSE)  # Works, but how can I mimic what stargazer does with lm objects?

更准确地说: 对于 lm 对象，占星家很好地打印了表格顶部的因变量，在相应的估计值下面的括号中包括 SE，并且在表格底部有 R ^ 2和观测数量。是否有一个(n)简单的方法来获得相同的行为与“自定义”模型估计的最大可能性，如上所述？

下面是我将 Optimoutput 打扮成 lm 对象的微弱尝试:

model2.lm <- list()  # Mimic an lm object
class(model2.lm) <- c(class(model2.lm), "lm")
model2.lm$rank <- model1$rank  # Problematic?
model2.lm$coefficients <- model2$par
names(model2.lm$coefficients)[1:2] <- names(model1$coefficients)
model2.lm$fitted.values <- model2$par["const"] + model2$par["beta"]*df$x
model2.lm$residuals <- df$y - model2.lm$fitted.values
model2.lm$model <- df
model2.lm$terms <- model1$terms  # Problematic?
summary(model2.lm)  # Not working
3247 次浏览
小开

I don't know how committed you are to using stargazer, but you can try using the broom and the xtable packages, the problem is that it won't give you the standard errors for the optim model 

library(broom)
library(xtable)
xtable(tidy(model1))
xtable(tidy(model2))
小开

I was just having this problem and overcame this through the use of the coef se, and omit functions within stargazer... e.g.

stargazer(regressions, ...
coef = list(... list of coefs...),
se = list(... list of standard errors...),
omit = c(sequence),
covariate.labels = c("new names"),
dep.var.labels.include = FALSE,
notes.append=FALSE), file="")
小开

You need to first instantiate a dummy lm object, then dress it up:

#...
model2.lm = lm(y ~ ., data.frame(y=runif(5), beta=runif(5), scale=runif(5), degrees.freedom=runif(5)))
model2.lm$coefficients <- model2$par
model2.lm$fitted.values <- model2$par["const"] + model2$par["beta"]*df$x
model2.lm$residuals <- df$y - model2.lm$fitted.values
stargazer(model2.lm, se = list(model2.coefs$se), summary=FALSE, type='text')


# ===============================================
#                         Dependent variable:
#                     ---------------------------
#                                  y
# -----------------------------------------------
# const                        10.127***
#                               (0.680)
#
# beta                         1.995***
#                               (0.024)
#
# scale                        3.836***
#                               (0.393)
#
# degrees.freedom              3.682***
#                               (1.187)
#
# -----------------------------------------------
# Observations                    200
# R2                             0.965
# Adjusted R2                    0.858
# Residual Std. Error       75.581 (df = 1)
# F Statistic              9.076 (df = 3; 1)
# ===============================================
# Note:               *p<0.1; **p<0.05; ***p<0.01

(and then of course make sure the remaining summary stats are correct)


提交答案