Która funkcja/pakiet dla niezawodnej regresji liniowej działa z glmulti (tj. Zachowuje się jak glm)?

Question

Wstęp: Wielu modelu wnioskowania z glmulti

glmulti jest funkcją R/pakiet do automatycznego doboru modelu do ogólnych modeli liniowych, które konstruuje wszystkich możliwych ogólnych modeli liniowych podane zmienną zależną oraz zestaw czynników predykcyjnych, wpisuje je poprzez klasyczną funkcję glm i pozwala wtedy na wnioskowanie wielomodelowe (np. przy użyciu wag modeli pochodzących od AICc, BIC). glmulti działa w teorii również z innymi funkcjami, które zwracają współczynniki, prawdopodobieństwo logu modelu i liczbę wolnych parametrów (i może inne informacje?) W tym samym formacie, co robi.

Mój cel: Multi-Model wnioskowanie z solidnych błędów

Chciałbym użyć glmulti z wytrzymałej modelowania błędów zmiennej ilościowej zależnej w celu ochrony przed wpływem out błędne.

Na przykład mógłbym założyć, że błędy w modelu liniowym są dystrybuowane jako t distribution zamiast jako rozkład normalny. Dzięki parametrowi kurtozy rozkład t może mieć duże ogony, a zatem jest bardziej odporny na wartości odstające (w porównaniu do rozkładu normalnego).

Nie jestem jednak zdecydowany użyć metody t dystrybucji. Cieszę się z każdego podejścia, które daje prawdopodobieństwo logowania, a zatem działa z podejściem wielomodelowym w glmulti. Ale to znaczy, że niestety nie mogą korzystać z dobrze znanych solidnych modeli liniowych w R (np lmRob od robust lub lmrob z robustbase), ponieważ nie działają w ramach log-prawdopodobieństwa, a więc nie może pracować z glmulti.

Problem: Nie mogę znaleźć solidnego funkcji regresji, która współpracuje z glmulti

Jedyne solidne funkcja regresji liniowej dla RI stwierdził, że działa zgodnie z ramami log-prawdopodobieństwo jest heavyLm (od heavy pakiet); modeluje błędy z rozkładem t. Niestety, heavyLm nie działa z glmulti (przynajmniej nie po wyjęciu z pudełka), ponieważ nie ma metody S3 dla loglik (i ewentualnie innych rzeczy).

Aby zilustrować:

library(glmulti) 
library(heavy) 

Używanie zestawu danych stackloss

head(stackloss) 

Regularne modelu Gaussa liniowy:

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

Wielu modelu wnioskowania z glmulti nas ing GLM „s domyślna funkcja Gaussa Link

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

Model liniowy z t rozprowadzane błędu (domyślnie df = 4)

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

Wielu modelu wnioskowania z glmulti wywołującego heavyLm jako funkcja dopasowania

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

daje następujący błąd:

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

Gdybym zdefiniować następującą funkcję,

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

glmulti można uzyskać dziennik prawdopodobieństwo, ale potem następny błąd:

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

pytanie: Która funkcja/pakiet dla solidnej liniowej regresji działa z glmulti (tj. Zachowuje się jak glm)?

Nie ma chyba sposób zdefiniować dalsze funkcje, aby uzyskać heavyLm pracy z glmulti, ale przed rozpoczęciem tej podróży Chciałem zapytać, czy ktokolwiek

• zna solidnego funkcji regresji liniowej że (a) działa zgodnie z zasadami wiarygodności logów i (b) zachowuje się jak glm (i dzięki temu będzie działać bezczynnie z glmulti).
• dostał ciężki Lm już działa z glmulti.

Każda pomoc jest bardzo cenna!

Answer 1

Oto odpowiedź z użyciem heavyLm. Pomimo tego, że jest to stosunkowo stare pytanie, ten sam problem, o którym wspomniałeś, nadal występuje, gdy używana jest metoda HeavyLm (tzn. Komunikat o błędzie Error in .jcall(molly, "V", "supplyErrorDF"…).

Problem polega na tym, że glmulti wymaga stopni swobody modelu, które należy przekazać jako atrybut, który należy podać jako atrybut wartości zwracanej przez funkcję logLik.heavyLm; szczegółowe informacje można znaleźć w dokumentacji funkcji logLik.Co więcej, okazuje się, że musisz również zapewnić funkcję zwracania liczby punktów danych, które zostały użyte do dopasowania modelu, ponieważ kryteria informacyjne (AIC, BIC, ...) również zależą od tej wartości. Odbywa się to za pomocą funkcji nobs.heavyLm w poniższym kodzie.

Oto kod:

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 
} 

które gdy wprowadzone wraz z kodem, który podałeś, produkuje następujący wynik:

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. 

produkującą dlatego 2 modele w progu 2 jednostki BIC .

Ważna uwaga: nie jestem pewien, czy powyższe wyrażenia dotyczące stopni swobody są ściśle prawidłowe. Dla standardowego modelu liniowego, stopnie swobody będą równe p + 1, gdzie p jest liczbą parametrów w modelu, a dodatkowy parametr (+ 1) jest wariancją "błędu" (która służy do obliczenia prawdopodobieństwa) . W funkcji logLik.heavyLm powyżej, nie jest dla mnie jasne, czy należy również liczyć "parametr skali", który jest oszacowany przez heavyLm jako dodatkowy stopień swobody, a zatem p + 1 + 1, co miałoby miejsce, gdyby prawdopodobieństwo było również funkcją tego parametru. Niestety, nie mogę tego potwierdzić, ponieważ nie mam dostępu do odniesienia, które cytuje cytowany dokument (praca Dempstera i wsp., 1980). Z tego powodu zliczam parametr skali, tym samym zapewniając (nieco więcej) ostrożne oszacowanie złożoności modelu, penalizując "złożone" modele. Różnica ta powinna być znikoma, z wyjątkiem przypadku małej próbki.