Holen Sie sich die durch die maximale Wahrscheinlichkeit geschätzten Koeffizienten in eine Sternguckertabelle

83

Stargazer produziert sehr schöne Latex-Tische für lm (und andere) Objekte. Angenommen, ich habe ein Modell mit maximaler Wahrscheinlichkeit angepasst. Ich möchte, dass Stargazer eine lm-ähnliche Tabelle für meine Schätzungen erstellt. Wie kann ich das machen?

Obwohl es ein bisschen hackig ist, könnte eine Möglichkeit darin bestehen, ein "falsches" lm-Objekt zu erstellen, das meine Schätzungen enthält. Ich denke, dies würde funktionieren, solange die Zusammenfassung (my.fake.lm.object) funktioniert. Ist das leicht machbar?

Ein Beispiel:

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?

Genauer gesagt: Bei lm-Objekten druckt stargazer die abhängige Variable oben in der Tabelle gut aus, enthält SEs in Klammern unter den entsprechenden Schätzungen und hat R ^ 2 und die Anzahl der Beobachtungen unten in der Tabelle. Gibt es eine (n einfache) Möglichkeit, dasselbe Verhalten mit einem "benutzerdefinierten" Modell zu erzielen, das wie oben anhand der maximalen Wahrscheinlichkeit geschätzt wird?

Hier sind meine schwachen Versuche, meine optimale Ausgabe als lm-Objekt zu verkleiden:

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
r optimization lm stargazer Adrian
quelle
6
Ich habe versucht, etwas Ähnliches mit dem texregPaket. Aufgrund der Faulheit habe ich die Koeffizienten und Standardfehler eines anderen Modells überschrieben, wodurch ich die gewünschte Ausgabe erhalten habe. In Ihrem Fall könnten Sie zB die Koeffizienten und Standardfehler von überschreiben model1. Dies ist zwar keine ausgefeilte Lösung, sollte aber funktionieren. Unnötig zu sagen, ich bin gespannt, ob es bessere Lösungen gibt ...
coffeinjunky
1
Sie können sich die Stargazer-Funktion ansehen, mit der das schwere Heben ausgeführt wird stargazer:::.stargazer.wrap. Es sieht aus wie ein Container mit einer Reihe weiterer Funktionen zusätzlich zu dem Code, der die Tabellen formatiert. Und es scheint, als würde es einige Komponenten für lm(und glm) bewerten , die es sehr schwierig machen würden, Ihre optim()Ergebnisse zu verbessern .
Andybega
3
In texregwürde es ausreichen, ein texregObjekt mit der createTexregFunktion zu erstellen . Sie übergeben im Grunde nur die Koeffizienten, SEs usw. Siehe ?createTexreg. Das texregObjekt kann dann in die zugeführt werden texreg, htmlreg, screenregund plotregFunktionen. Alternativ wird in Abschnitt 6 des JSS-Artikels beschrieben, wie Methoden für neue Modelltypen geschrieben und registriert werden, falls Sie dieselbe Vorlage später wiederverwenden möchten.
Philip Leifeld

Antworten:

2

Ich hatte gerade dieses Problem und habe es durch die Verwendung von coef seund omitFunktionen innerhalb von Stargazer überwunden ... z

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="")
Biege
quelle
1

Sie müssen zuerst ein Dummy- lmObjekt instanziieren und es dann verkleiden:

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

(und dann natürlich sicherstellen, dass die verbleibenden zusammenfassenden Statistiken korrekt sind)

luffe
quelle
0

Ich weiß nicht, wie sehr Sie sich für die Verwendung von Stargazer engagieren, aber Sie können versuchen, den Besen und die xtable-Pakete zu verwenden. Das Problem ist, dass Sie nicht die Standardfehler für das Optim-Modell erhalten 

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