Welche Funktion/Paket für robuste lineare Regression funktioniert mit Glmulti (d. H. Verhält sich wie Glm)?

Question

Hintergrund: Multimodell-Inferenz mit glmulti

glmulti ist ein R-Funktion/Paket für automatisierte Modellauswahl für allgemeine lineare Modelle, die eine abhängige Variable und einen Satz von Prädiktoren gegeben möglich Lineares Modell konstruiert, paßt sie über die klassische glm-Funktion und ermöglicht dann eine Multi-Modell-Inferenz (z. B. unter Verwendung von Modellgewichten, die von AICc, BIC abgeleitet sind). glmulti funktioniert in der Theorie auch mit jeder anderen Funktion, die Koeffizienten, die Log-Likelihood des Modells und die Anzahl der freien Parameter (und vielleicht andere Informationen?) Im gleichen Format, das Glm tut.

Mein Ziel: Mehr Modell Inferenz mit robusten Fehlern

Ich mag glmulti mit robuster Modellierung des Fehlers einer quantitativen abhängigen Variablen verwenden, gegen die Wirkung aus Ausreißern zu schützen.

Zum Beispiel könnte ich davon ausgehen, dass die Fehler im linearen Modell als t distribution anstatt als eine normale Verteilung verteilt sind. Mit ihrem Kurtosis-Parameter kann die t-Verteilung schwere Schwänze aufweisen und ist daher robuster gegenüber Ausreißern (im Vergleich zur Normalverteilung).

Allerdings bin ich nicht verpflichtet, den T-Verteilungsansatz zu verwenden. Ich bin mit jedem Ansatz zufrieden, der eine Log-Likelihood zurückgibt und somit mit dem Multimodelansatz in glmulti arbeitet. Aber das bedeutet, dass ich kann leider nicht die bekannten robusten linearen Modelle in R verwenden (zB lmRob von robust oder lmrob von robustbase), weil sie unter dem Log-Likelihood-Rahmen arbeiten nicht und somit kann nicht funktionieren mit glmulti.

Das Problem: Ich kann nicht eine robuste Regressionsfunktion, die mit arbeitet finden glmulti

Die einzige robuste Funktion lineare Regression für RI festgestellt, dass arbeitet unter dem Log-Likelihood Rahmen ist heavyLm (von der heavy Paket); Es modelliert die Fehler mit einer t-Verteilung. Leider funktioniert healthLm nicht mit glmulti (zumindest nicht out of the box), weil es keine S3-Methode für loglik (und möglicherweise andere Dinge) hat.

Zur Veranschaulichung:

summary(glm(stack.loss ~ ., data = stackloss)) 

Multi-Modell-Inferenz mit glmulti uns:

library(glmulti) 
library(heavy) 

den Datensatz

head(stackloss) 

Regular Gaußsche lineare Modell stackloss Verwendung ing glm ‚s Standard Gaussian Verknüpfungsfunktion

stackloss.glmulti <- glmulti(stack.loss ~ ., data = stackloss, level=1, crit=bic) 
print(stackloss.glmulti) 
plot(stackloss.glmulti) 

Linear-Modell mit Fehlern t verteilt (default df = 4)

summary(heavyLm(stack.loss ~ ., data = stackloss)) 

Multimodell-Inferenz mit glmulti Aufruf heavyLm als Die Anpassungsfunktion

stackloss.heavyLm.glmulti <- glmulti(stack.loss ~ ., 
data = stackloss, level=1, crit=bic, fitfunction=heavyLm) 

gibt der folgende Fehler:

Initialization... 
    Error in UseMethod("logLik") : 
    no applicable method for 'logLik' applied to an object of class "heavyLm". 

Wenn ich die folgende Funktion definieren,

logLik.heavyLm <- function(x){x$logLik} 

glmulti kann das Log-Likelihood, erhalten, aber dann ist der nächste Fehler auftritt:

Initialization... 
    Error in .jcall(molly, "V", "supplyErrorDF", 
    as.integer(attr(logLik(fitfunc(as.formula(paste(y, : 
    method supplyErrorDF with signature ([I)V not found 

Die Frage: Welche Funktion/Paket für robuste lineare Regression arbeitet mit glmulti (dh verhält sich wie glm)?

Es gibt wohl eine Möglichkeit, weitere Funktionen zu definieren heavyLm zum Laufen zu bringen mit glmulti, aber vor dieser Reise einsteigen wollte ich, ob jemand fragen

• kennt eine robusten Funktion linearen Regression, dass (a) arbeitet unter dem Log-Likelihood-Framework und (b) verhält sich wie glm (und wird daher mit glmulti Out-of-the-Box arbeiten).
• haben healthLm bereits arbeiten mit glmulti.

Jede Hilfe wird sehr geschätzt!

Answer 1

Hier ist eine Antwort mit heavyLm. Obwohl dies eine relativ alte Frage ist, tritt das gleiche Problem, das Sie erwähnten, immer noch auf, wenn Sie heavyLm verwenden (d. H. Die Fehlermeldung Error in .jcall(molly, "V", "supplyErrorDF"…).

Das Problem ist, dass glmulti erfordert die Freiheitsgrade des Modells, um als ein Attribut von Ihnen als Attribut des Werts zurückgegeben werden zurückgegeben von der Funktion ; Einzelheiten finden Sie in der Dokumentation zur Funktion logLik.Außerdem müssen Sie eine Funktion bereitstellen, die die Anzahl der Datenpunkte zurückgibt, die für die Anpassung des Modells verwendet wurden, da die Informationskriterien (AIC, BIC, ...) auch von diesem Wert abhängen. Dies geschieht durch die Funktion nobs.heavyLm im folgenden Code. Hier

ist der Code:

nobs.heavyLm <- function(mdl) mdl$dims[1] # the sample size (number of data points) 

logLik.heavyLm <- function(mdl) { 
    res <- mdl$logLik 
    attr(res, "nobs") <- nobs.heavyLm(mdl) # this is not really needed for 'glmulti', but is included to adhere to the format of 'logLik' 
    attr(res, "df") <- length(mdl$coefficients) + 1 + 1 # I am also considering the scale parameter that is estimated; see mdl$family 
    class(res) <- "logLik" 
    res 
} 

die, wenn sie zusammen mit dem Code, den Sie zur Verfügung gestellt, das folgende Ergebnis erhalten:

Initialization... 
TASK: Exhaustive screening of candidate set. 
Fitting... 
Completed. 

> print(stackloss.glmulti) 
glmulti.analysis 
Method: h/Fitting: glm/IC used: bic 
Level: 1/Marginality: FALSE 
From 8 models: 
Best IC: 117.892471265874 
Best model: 
[1] "stack.loss ~ 1 + Air.Flow + Water.Temp" 
Evidence weight: 0.709174196998897 
Worst IC: 162.083142797858 
2 models within 2 IC units. 
1 models to reach 95% of evidence weight. 

daher zwei Modelle innerhalb der 2 BIC Einheiten Schwelle Herstellung .

Eine wichtige Bemerkung: Ich bin mir nicht sicher, ob der obige Ausdruck für die Freiheitsgrade genau richtig ist. Für ein standardmäßiges lineares Modell wären die Freiheitsgrade gleich p + 1, wobei p die Anzahl der Parameter im Modell ist und der zusätzliche Parameter (+ 1) die "Fehlervarianz" ist (die zur Berechnung der Wahrscheinlichkeit verwendet wird). . In obiger Funktion ist es für mich nicht klar, ob man auch den "Skalenparameter", der durch heavyLm geschätzt wird, als einen zusätzlichen Freiheitsgrad zählen sollte, und somit die p + 1 + 1, was der Fall wäre, wenn die Wahrscheinlichkeit auch eine Funktion ist von diesem Parameter. Leider kann ich das nicht bestätigen, da ich keinen Zugang zu der Referenz habe, die heavyLm zitiert (die Veröffentlichung von Dempster et al., 1980). Aus diesem Grund zähle ich den Skalenparameter und sorge so für eine (etwas mehr) konservative Schätzung der Modellkomplexität, die "komplexe" Modelle benachteiligt. Dieser Unterschied sollte vernachlässigbar sein, außer im kleinen Beispielfall.