Uzyskaj współczynniki oszacowane według maksymalnego prawdopodobieństwa do tabeli obserwatorów gwiazd

Stargazer tworzy bardzo ładne lateksowe tabele dla lm (i innych) obiektów. Załóżmy, że dopasowałem model według największego prawdopodobieństwa. Chciałbym, aby Stargazer utworzył tabelę podobną do lm do moich szacunków. Jak mogę to zrobić?

Chociaż jest to trochę skomplikowane, jednym ze sposobów może być utworzenie „fałszywego” obiektu lm zawierającego moje szacunki - myślę, że będzie to działać tak długo, jak długo będzie działać podsumowanie (my.fake.lm.object). Czy jest to łatwe do wykonania?

Przykład:

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?

Dokładniej: z obiektami lm, Stargazer ładnie drukuje zmienną zależną na górze tabeli, zawiera SE w nawiasach poniżej odpowiednich oszacowań i ma R ^ 2 i liczbę obserwacji na dole tabeli. Czy istnieje (n łatwy) sposób na uzyskanie takiego samego zachowania z „niestandardowym” modelem szacowanym według największego prawdopodobieństwa, jak powyżej?

Oto moje słabe próby ubierania mojego wyjścia optymalnego jako obiektu 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

Właśnie miałem ten problem i przezwyciężyłem go dzięki użyciu funkcji coef sei omitw programie Stargazer ... np

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

Musisz najpierw utworzyć instancję lmobiektu zastępczego , a następnie go ubierać:

#...
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

(a potem oczywiście upewnij się, że pozostałe statystyki podsumowujące są poprawne)

Nie wiem, jak bardzo jesteś zaangażowany w używanie Stargazera, ale możesz spróbować użyć miotły i pakietów xtable, problem w tym, że nie da ci to standardowych błędów dla modelu optim

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