Quale funzione/pacchetto per la regressione lineare robusta funziona con glmulti (cioè si comporta come glm)?

Question

Priorità: Multi-modello inferenza con glmulti

glmulti è una funzione R/pacchetto per selezione del modello automatizzato per modelli lineari generali che costruisce tutti i possibili modelli lineari generali dati una variabile dipendente e un insieme di predittori, li inserisce tramite la classica funzione glm e consente quindi l'inferenza su più modelli (ad es. utilizzando pesi modello derivati ​​da AICc, BIC). glmulti funziona in teoria anche con qualsiasi altra funzione che restituisca i coefficienti, la log-verosimiglianza del modello e il numero di parametri liberi (e forse altre informazioni?) Nello stesso formato che fa glm.

Il mio obiettivo: Multi-modello di inferenza con errori robusti

vorrei utilizzare glmulti con la modellazione solida degli errori di una variabile dipendente quantitativa in guardia contro l'effetto fuori valori anomali.

Ad esempio, posso supporre che gli errori nel modello lineare siano distribuiti come t distribution anziché come una distribuzione normale. Con il parametro kurtosis, la distribuzione t può avere code pesanti ed è quindi più robusta rispetto ai valori anomali (rispetto alla distribuzione normale).

Tuttavia, non sono impegnato a utilizzare l'approccio di distribuzione t. Sono contento di qualsiasi approccio che restituisca una probabilità di log e quindi funziona con l'approccio multimodel in glmulti. Ma questo significa, che purtroppo non posso usare i ben noti modelli lineari robusti in R (ad esempio, lmRob da robust o lmrob da robustbase), perché non operano nel quadro log-verosimiglianza e quindi possono non funzionare con glmulti.

Il problema: non riesco a trovare una funzione di regressione robusto che funziona con glmulti

L'unica robusta funzione di regressione lineare per RI trovato che opera nel quadro di log-verosimiglianza è heavyLm (dal Pacchetto heavy); modella gli errori con una distribuzione t. Sfortunatamente, heavyLm non funziona con glmulti (almeno non fuori dalla scatola) perché non ha il metodo S3 per loglik (e probabilmente altre cose).

Facciamo un esempio:

library(glmulti) 
library(heavy) 

Utilizzando il set di dati stackloss

head(stackloss) 

regolare modello lineare gaussiano:

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

Multi-modello di inferenza con glmulti noi ing glm s' default funzione di collegamento gaussiana

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

modello lineare con t distribuito errore (di default è df = 4)

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

Multi-modello di inferenza con glmulti chiamando heavyLm come la funzione di raccordo

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

fornisce il seguente errore:

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

Se io definisco la seguente funzione,

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

glmulti può ottenere il log-verosimiglianza, ma poi l'errore successivo si verifica:

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

La domanda: Quale funzione/pacchetto per una regressione lineare robusta funziona con glmulti (cioè si comporta come glm)?

C'è probabilmente un modo per definire ulteriori funzioni per ottenere heavyLm lavorare con glmulti, ma prima di intraprendere questo viaggio volevo chiedere se qualcuno

• sa di una robusta funzione di regressione lineare che (a) opera sotto il quadro log-verosimiglianza e (b) si comporta come glm (e funzionerà quindi con glmulti out-of-the-box).
• ottenuto heavyLm già funzionante con glmulti.

Qualsiasi aiuto è molto apprezzato!

Answer 1

Questa è una risposta utilizzando heavyLm. Anche se si tratta di una domanda relativamente vecchia, lo stesso problema che hai menzionato si verifica ancora quando si utilizza heavyLm (ad esempio, il messaggio di errore Error in .jcall(molly, "V", "supplyErrorDF"…).

Il problema è che glmulti richiede i gradi di libertà del modello, da passare come attributo di cui è necessario fornire come attributo del valore restituito dalla funzione logLik.heavyLm; consultare la documentazione per la funzione logLik per i dettagli.Inoltre, risulta necessario anche fornire una funzione per restituire il numero di punti dati utilizzati per il montaggio del modello, poiché i criteri di informazione (AIC, BIC, ...) dipendono anche da questo valore. Questo è fatto dalla funzione nobs.heavyLm nel codice qui sotto.

Ecco il codice:

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 
} 

che, quando messo insieme con il codice che hai fornito, produce il seguente risultato:

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. 

producendo quindi 2 modelli entro la soglia 2 unità BIC .

Un'osservazione importante però: non sono sicuro che l'espressione sopra per i gradi di libertà sia strettamente corretta. Per un modello lineare standard, i gradi di libertà sarebbero uguali a p + 1, dove p è il numero di parametri nel modello e il parametro extra (il + 1) è la varianza "errore" (che viene utilizzata per calcolare la probabilità) . Nella funzione logLik.heavyLm di cui sopra, non è chiaro se si debba anche contare il "parametro di scala" stimato da heavyLm come ulteriore grado di libertà, e quindi lo p + 1 + 1, che sarebbe il caso se la verosimiglianza fosse anche una funzione di questo parametro. Sfortunatamente, non posso confermarlo, poiché non ho accesso al riferimento che cita heavyLm (il documento di Dempster et al., 1980). Per questo motivo, sto contando il parametro di scala, fornendo quindi una stima (leggermente più) conservativa della complessità del modello, penalizzando i modelli "complessi". Questa differenza dovrebbe essere trascurabile, tranne nel caso di un piccolo campione.