Vorhersage-Methode für Felm aus lfe-Paket

Question

Hat jemand eine schöne saubere Möglichkeit, predict Verhalten für felm Modelle zu bekommen?

library(lfe) 
model1 <- lm(data = iris, Sepal.Length ~ Sepal.Width + Species) 
predict(model1, newdata = data.frame(Sepal.Width = 3, Species = "virginica")) 
# Works 

model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species) 
predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica")) 
# Does not work 

Answer 1

Ich denke, was Sie suchen, könnte das lme4 Paket sein. Ich war in der Lage ein, um Arbeit zu bekommen vorhersagen, mit diesem:

library(lme4) 
data(iris) 

model2 <- lmer(data = iris, Sepal.Length ~ (Sepal.Width | Species)) 
predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica")) 
     1 
6.610102 

Sie können sich ein wenig spielen, müssen die besonderen Effekte für Sie suchen, aber das Paket ist gut dokumentiert, so sollte es nicht ein Problem sein.

Answer 2

Dies ist möglicherweise nicht die Antwort, die Sie suchen, aber es scheint, dass der Autor keine Funktionalität hinzugefügt lfe Paket, um Vorhersagen über externe Daten mit dem angepassten felm Modell. Der primäre Fokus scheint auf der Analyse der Gruppen-Fixed-Effects zu liegen. Aber es ist interessant, dass in der Dokumentation des Pakets beachten Sie die folgenden erwähnt wird:

The object has some resemblance to an 'lm' object, and some postprocessing methods designed for lm may happen to work. It may however be necessary to coerce the object to succeed with this.

Daher könnte es möglich sein, die felm Objekt zu einem lm Objekt, um zu zwingen, einige zusätzliche lm Funktionalität zu erhalten (wenn alle erforderlichen Informationen sind in dem Objekt vorhanden, um die notwendigen Berechnungen durchzuführen).

Das lfe Paket soll sehr große Datenmengen lief werden und Anstrengungen unternommen wurden, um Speicher zu sparen: Als direkte Folge davon, das felm Objekt nicht/nicht verwenden, um eine QR-Zerlegung enthält, wie zum lm Objekt gegenüber. Leider beruht die lmpredict Prozedur auf dieser Information, um die Vorhersagen zu berechnen. Daher Nötigung die felm Objekt und die Methode fehl vorhersagen Ausführung:

> model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species) 
> class(model2) <- c("lm","felm") # coerce to lm object 
> predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica")) 
Error in qr.lm(object) : lm object does not have a proper 'qr' component. 
Rank zero or should not have used lm(.., qr=FALSE). 

Wenn Sie wirklich dieses Paket die Prognosen ausführen verwenden müssen, dann können Sie vielleicht Ihre eigene vereinfachte Version dieser Funktionalität schreiben, indem die Informationen, die Sie verfügbar im felm Objekt. Zum Beispiel sind die OLS-Regressionskoeffizienten über model2$coefficients verfügbar.

Answer 3

Als Abhilfe können, könnten Sie felm, getfe und demeanlist kombinieren wie folgt:

library(lfe) 

lm.model <- lm(data=demeanlist(iris[, 1:2], list(iris$Species)), Sepal.Length ~ Sepal.Width) 
fe <- getfe(felm(data = iris, Sepal.Length ~ Sepal.Width | Species)) 
predict(lm.model, newdata = data.frame(Sepal.Width = 3)) + fe$effect[fe$idx=="virginica"] 

Die Idee, dass Sie die Variablen zum Zentrum verwenden demeanlist ist, dann lm den Koeffizienten auf Sepal.Width mit der zentriert zu schätzen Variablen, geben Sie ein lm Objekt, über das Sie predict ausführen können. Führen Sie dann felm + getfe aus, um den bedingten Mittelwert für den festen Effekt zu erhalten, und fügen Sie diesen Wert zur Ausgabe von predict hinzu.

Answer 4

Dies sollte für Fälle funktionieren, in denen Sie die Gruppeneffekte in der Vorhersage ignorieren, für neue Xs voraussagen und nur Vertrauensbereiche wünschen. Es sucht zuerst nach einem clustervcv Attribut, dann robustvcv, dann vcv.

predict.felm <- function(object, newdata, se.fit = FALSE, 
         interval = "none", 
         level = 0.95){ 
    if(missing(newdata)){ 
    stop("predict.felm requires newdata and predicts for all group effects = 0.") 
    } 

    tt <- terms(object) 
    Terms <- delete.response(tt) 
    attr(Terms, "intercept") <- 0 

    m.mat <- model.matrix(Terms, data = newdata) 
    m.coef <- as.numeric(object$coef) 
    fit <- as.vector(m.mat %*% object$coef) 
    fit <- data.frame(fit = fit) 

    if(se.fit | interval != "none"){ 
    if(!is.null(object$clustervcv)){ 
     vcov_mat <- object$clustervcv 
    } else if (!is.null(object$robustvcv)) { 
     vcov_mat <- object$robustvcv 
    } else if (!is.null(object$vcv)){ 
     vcov_mat <- object$vcv 
    } else { 
     stop("No vcv attached to felm object.") 
    } 
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat))) 
    } 
    if(interval == "confidence"){ 
    t_val <- qt((1 - level)/2 + level, df = object$df.residual) 
    fit$lwr <- fit$fit - t_val * se.fit_mat 
    fit$upr <- fit$fit + t_val * se.fit_mat 
    } else if (interval == "prediction"){ 
    stop("interval = \"prediction\" not yet implemented") 
    } 
    if(se.fit){ 
    return(list(fit=fit, se.fit=se.fit_mat)) 
    } else { 
    return(fit) 
    } 
}