2012-08-13 8 views
5

Antecedentes: Multi-modelo inferencia con glmulti

glmulti es una función R/paquete para la selección del modelo automatizado para los modelos lineales generales que construye todos los posibles modelos lineales generales dadas una variable dependiente y un conjunto de predictores, ellos encaja a través de la función clásica glm y permite entonces la inferencia de modelos múltiples (por ejemplo, utilizando pesos modelo derivados de AICc, BIC). glmulti funciona en teoría también con cualquier otra función que arroje coeficientes, la logaritmo-verosimilitud del modelo y el número de parámetros libres (y tal vez otra información?) En el mismo formato que glm.¿Qué función/paquete para la regresión lineal robusta funciona con glmulti (es decir, se comporta como glm)?

Mi objetivo: la inferencia multi-modelo con errores robustos

me gustaría utilizar glmulti con el modelado robusta de los errores de una variable dependiente cuantitativa para protegerse contra el efecto cabo valores atípicos.

Por ejemplo, podría suponer que los errores en el modelo lineal se distribuyen como t distribution en lugar de como una distribución normal. Con su parámetro de curtosis, la distribución t puede tener colas pesadas y, por lo tanto, es más robusta a los valores atípicos (en comparación con la distribución normal).

Sin embargo, no me comprometo a usar el enfoque de distribución t. Estoy contento con cualquier enfoque que devuelva una probabilidad logarítmica y por lo tanto funciona con el enfoque multimodal en glmulti. Pero eso significa, que, por desgracia no puedo usar los modelos lineales sólidas bien conocidas en R (por ejemplo, lmRob de robust o lmrob de robustbase) debido a que no operan bajo el marco de probabilidad logarítmica y por lo tanto puede no funcionar con glmulti.

El problema: No puedo encontrar una función de regresión robusta que trabaja con glmulti

La única robusta función de regresión lineal para RI encontró que opera bajo el marco de probabilidad logarítmica es heavyLm (desde el heavy paquete); modela los errores con una distribución t. Por desgracia, heavyLm no funciona con glmulti (al menos no fuera de la caja), ya que no tiene un método para S3 loglik (y posiblemente otras cosas).

Para ilustrar:

library(glmulti) 
library(heavy) 

Utilizando el conjunto de datos stackloss

head(stackloss) 

Regular modelo gaussiano lineal:

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

multi-modelo de inferencia con glmulti nos ing GLM 's predeterminado función de enlace Gaussian

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

modelo lineal con t distribuido error (por defecto se df = 4)

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

Multi-modelo inferencia con glmulti llamando heavyLm como la función de ajuste

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

da el siguiente error:

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

Si defino la siguiente función,

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

glmulti puede obtener el logaritmo de la verosimilitud, pero entonces se produce el siguiente error:

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 pregunta: ¿Qué función/paquete para la regresión lineal robusta funciona con glmulti (es decir, se comporta como glm)?

Probablemente hay una manera de definir otras funciones para obtener heavyLm trabajar con glmulti, pero antes de embarcarse en este viaje que quería preguntar si alguien

  • sabe de una función de regresión lineal robusta que (a) opera bajo el marco log-verosimilitud y (b) se comporta como glm (y por lo tanto funcionará con glmulti out-of-the-box).
  • tengo heavyLm que ya está trabajando con glmulti.

¡Toda ayuda es muy apreciada!

Respuesta

1

Aquí hay una respuesta usando heavyLm. Aunque esta es una pregunta relativamente antigua, el mismo problema que mencionó aún se produce cuando se usa heavyLm (es decir, el mensaje de error Error in .jcall(molly, "V", "supplyErrorDF"…).

El problema es que glmulti requiere los grados de libertad del modelo, para pasar como un atributo de lo que necesita proporcionar como un atributo del valor devuelto por la función logLik.heavyLm; ver la documentación para la función logLik para más detalles.Además, resulta que también debe proporcionar una función para devolver la cantidad de puntos de datos que se usaron para ajustar el modelo, ya que los criterios de información (AIC, BIC, ...) también dependen de este valor. Esto se hace con la función nobs.heavyLm en el siguiente código.

Aquí está el código:

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 
} 

que, cuando se pone junto con el código que ha proporcionado, produce el siguiente resultado:

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. 

produciendo por lo tanto 2 modelos dentro del umbral de 2 unidades BIC .

Sin embargo, una observación importante: no estoy seguro de que la expresión anterior para los grados de libertad sea estrictamente correcta. Para un modelo lineal estándar, los grados de libertad serían iguales a p + 1, donde p es el número de parámetros en el modelo, y el parámetro adicional (+ 1) es la varianza de "error" (que se usa para calcular la probabilidad) . En la función logLik.heavyLm anterior, no me queda claro si también se debe contar el "parámetro de escala" estimado por heavyLm como un grado extra de libertad, y por lo tanto el p + 1 + 1, que sería el caso si la probabilidad también es una función de este parámetro Desafortunadamente, no puedo confirmar esto, ya que no tengo acceso a la referencia que heavyLm cita (el documento de Dempster et al., 1980). Debido a esto, estoy contando el parámetro de escala, proporcionando así una estimación (un poco más) conservadora de la complejidad del modelo, penalizando los modelos "complejos". Esta diferencia debe ser insignificante, excepto en el caso de muestra pequeña.

+0

¡Muchas gracias! – jonlemon

Cuestiones relacionadas