0

Estoy intentando obtener el valor AIC de mis modelos lineales generalizados mixtos (GLMM). Para mis modelos el valor de maceración (Grinding) y el tipo de kit (Kit) son factores fijos, el tipo de especie (Especie)es un factor aleatorio y la Conc y A1 y A2 son las variables respuesta. He realizado un modelo para cada variable respuesta. El modelo para la variable respuesta concentración sigue una distribución gamma logarítmica, los demás una distribución gausiana.

Una vez que realizo los modelos usando este código en R:

         library(lme4)

        DATA$Especie<-as.factor(DATA$Especie)
        DATA$Grinding<-as.factor(DATA$Grinding)
        DATA$Kit<-as.factor(DATA$Kit)

        glmer.conc.Grinding_Kit <- glmer(Conc ~ Grinding * Kit + (1|Especie),
                      family = Gamma(link = "log"),data = DATA)

        glmer.A1.Grinding_Kit <- lmer(A1 ~ Grinding * Kit + (1|Especie),
                            data = DATA)


        glmer.A2.Grinding_Kit <- lmer(A2 ~ Grinding * Kit + (1|Especie),
                            data = DATA)

y veo el summary:

          summary(glmer.conc.Grinding_Kit)

no obtengo el valor del delta AIC para el primer modelo

         Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
           Family: Gamma  ( log )
         Formula: Conc ~ Grinding * Kit + (1 | Especie)
            Data: DATA

              AIC      BIC   logLik deviance df.resid 
           3144.2   3183.8  -1561.1   3122.2      259 

         Scaled residuals: 
             Min      1Q  Median      3Q     Max 
         -0.8001 -0.5579 -0.2418  0.2086 11.4961 

         Random effects:
          Groups   Name        Variance Std.Dev.
          Especie  (Intercept) 0.8366   0.9147  
          Residual             1.4596   1.2081  
         Number of obs: 270, groups:  Especie, 6

         Fixed effects:
                          Estimate Std. Error t value Pr(>|z|)    
         (Intercept)        5.1652     0.3493  14.787  < 2e-16 ***
         GrindingII        -0.2720     0.2312  -1.176  0.23945    
         GrindingIII       -0.5097     0.2312  -2.204  0.02749 *  
         KitB              -0.0786     0.2306  -0.341  0.73324    
         KitC               0.3682     0.2317   1.589  0.11205    
         GrindingII:KitB   -0.8480     0.3263  -2.599  0.00936 ** 
         GrindingIII:KitB  -0.9068     0.3270  -2.773  0.00555 ** 
         GrindingII:KitC   -0.2596     0.3365  -0.771  0.44042    
         GrindingIII:KitC  -0.3282     0.3258  -1.007  0.31379    
         ---
         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         Correlation of Fixed Effects:
                     (Intr) GrndII GrnIII KitB   KitC   GII:KB GIII:KB GII:KC
         GrindingII  -0.329                                                  
         GrindingIII -0.328  0.504                                           
         KitB        -0.325  0.493  0.492                                    
         KitC        -0.330  0.492  0.487  0.490                             
         GrndngII:KB  0.228 -0.699 -0.353 -0.705 -0.342                      
         GrndnIII:KB  0.226 -0.344 -0.698 -0.705 -0.341  0.504               
         GrndngII:KC  0.222 -0.693 -0.350 -0.354 -0.672  0.490  0.247        
         GrndnIII:KC  0.231 -0.356 -0.707 -0.353 -0.695  0.253  0.499   0.497

y para los demás modelos ni tampoco el valor de AIC:

              summary(glmer.A1.Grinding_Kit)
              Linear mixed model fit by REML ['lmerMod']
              Formula: A1 ~ Grinding * Kit + (1 | Especie)
                 Data: DATA

              REML criterion at convergence: 337.4

              Scaled residuals: 
                  Min      1Q  Median      3Q     Max 
              -2.0435 -0.3311 -0.0049  0.2549 14.2163 

              Random effects:
               Groups   Name        Variance Std.Dev.
               Especie  (Intercept) 0.05648  0.2377  
               Residual             0.18008  0.4244  
              Number of obs: 270, groups:  Especie, 6

              Fixed effects:
                               Estimate Std. Error t value
              (Intercept)       1.79000    0.12416  14.416
              GrindingII       -0.03333    0.10957  -0.304
              GrindingIII      -0.13333    0.10957  -1.217
              KitB              0.20667    0.10957   1.886
              KitC              0.24000    0.10957   2.190
              GrindingII:KitB  -0.03333    0.15495  -0.215
              GrindingIII:KitB -0.03667    0.15495  -0.237
              GrindingII:KitC   0.22667    0.15495   1.463
              GrindingIII:KitC  0.12333    0.15495   0.796

              Correlation of Fixed Effects:
                          (Intr) GrndII GrnIII KitB   KitC   GII:KB GIII:KB GII:KC
              GrindingII  -0.441                                                  
              GrindingIII -0.441  0.500                                           
              KitB        -0.441  0.500  0.500                                    
              KitC        -0.441  0.500  0.500  0.500                             
              GrndngII:KB  0.312 -0.707 -0.354 -0.707 -0.354                      
              GrndnIII:KB  0.312 -0.354 -0.707 -0.707 -0.354  0.500               
              GrndngII:KC  0.312 -0.707 -0.354 -0.354 -0.707  0.500  0.250        
              GrndnIII:KC  0.312 -0.354 -0.707 -0.354 -0.707  0.250  0.500   0.500

y lo mismo para el tercer modelo (tiene pinta que al ser familia gaussina ha tenido que cambiar algo que no permite que me devuelva el valor AIC.

He intentado buscar información pero no encuentro nada. Por otro lado he intentado una manera alternativa de obtener estos resultados:

           Cand.models <- list( )

           Cand.models[[1]] <- glmer(Conc ~ Grinding * Kit + (1|Especie),
                      family = Gamma(link = "log"),data = DATA)


           Cand.models[[2]] <- lmer(A1 ~ Grinding * Kit + (1|Especie),
                     data = DATA)

           Cand.models[[3]] <- lmer(A2 ~ Grinding * Kit + (1|Especie),
                     data = DATA)



           Modnames <- paste("mod", 1:length(Cand.models), sep = " ")

           aictab(cand.set=Cand.models,modnames=modnames,second.ord=true,nobs=null,sort=true) 

Pero me dice esto:

            Error in formatCands(cand.set) : 
            Functions do not support mixture of model classes

Me gustaría saber como obtener estos valores pues necesito saber cual es el valor de AIC y delta AIC de todos los modelos. ¿Qué estoy haciendo mal? ¿Qué maneras alternativas hay de hacerlo?

Gracias de antemano.

1 respuesta 1

1

No pude reproducir el código porque no tengo los datos, sin embargo mirando la documentación de la función y algunos hilos sobre el tema te comento lo que encontré.

  1. El delta AIC es la diferencia entre los AIC de dos o más modelos. No se puede calcular para un solo modelo. Esto es muy razonable considerando que el AIC es, en su uso práctico, una medida relativa. Sirve para comparar el ajuste de dos o más modelos con especificaciones muy similares, penalizando al que estima más cantidad de parámetros. Generalmente modelos iguales a los que agregamos/quitamos predictores o les cambiamos el lugar en la ecuación (fijo, aleatorio, etc). El valor en sí mismo no dice nada interpretable, solo tiene sentido en comparación con otro.

  2. Por lo que leí en esta buena respuesta https://stats.stackexchange.com/a/131280/136344 lmer no reporta AIC cuando se usa Restricted Maximum Likelihood Estimation (RMEL). Por defecto ese argumento está como TRUE. Podrías probar con RMEL = FALSE y entonces el sumario sí regresaría los criterios de información Bayesiano y de Akaike.

  3. Creo que el problema que tienes con aictab es el que indica el mensaje de error: no se pueden comparar modelos de clases diferentes. De manera estrictamente técnica y de R estás poniendo objetos de clases diferentes: unos son de la clase lmerMod y otros de la clase glmerMod. En cualquier caso siempre se puede obtener el valor de AIC con la función AIC(modelo_1, modelo_2) o de manera más rústica pero efectiva compara a los modelos con ANOVA. Por sólidos motivos estadísticos no es posible "engañar" a glmer para que ajuste un generalizado con distribución gaussiana, forzando al output a pertenecer a la clase glmerMod.

Creo que con estas alternativas (RMEL = FALSE y usar glmer aún para los gaussianos) resolverías el problema técnico. Muchas suerte.

En nota aparte y quizás fuera del tema de este foro estarían los problemas propiamente estadísticos. Decía más arriba que el AIC/BIC se usa para comparar modelos con especificaciones muy similares. Al estar comparando modelos con variables dependientes que miden cosas diferentes en mi opinión la comparación por AIC no tiene mucho sentido. Admito que no sé cuál es tu objetivo de investigación, quizás en ese contexto sí tenga sentido saber cuál es la variable dependiente que mejor ajuste tiene para una especificación de predictores siempre igual. Aún en ese caso no ese caso no sería necesario aplicar la penalización del AIC, porque no estás cambiando la cantidad de parámetros que ajustas para cada modelo.

3
  • Muchas gracias @mpaladino por esta pregunta. He seguido tus consejos y al añadir REML=FALSE obtengo el valor AIC. Respecto al aictab, sigo obteniendo el mismo error "Error in formatCands(cand.set) : Functions do not support mixture of model classe" . Muchas gracias por la respuesta y por describirla tan detalladamente. No entiendo porque te ha evaluado alguien del foro ajeno a mí con -1 cuando en me has resuelto la duda. Commented el 7 mar. 2021 a las 9:31
  • 1
    Revisé con los datos del ejemplo de lme4::glmer() y en efecto el output de la clase el lmerMod y no glmerMod. De hecho emite un warning muy informativo al respecto. El error con aictab va a persistir, voy a editar la respuesta. Más allá de este problema técnico y de programación, sigue el problema estadístico de fondo: no tiene sentido matemático comparar estos dos modelos por AIC/BIC.
    – mpaladino
    Commented el 7 mar. 2021 a las 16:54
  • 1
    Sí, vi que la respuesta fue reportada como de baja calidad, aunque finalmente invalidaron el reporte. Es un problema que se presenta bastante en la etiqueta R, las preguntas muchas veces no son de programación "pura" y suelen tener partes importantes que son problemas estadísticos. Yo siempre que puedo trato de responder igual, sobre todo porque no hay un SO en español específico de estadística (como CrossValidated). Creo que sería un buen tema para plantear en meta.
    – mpaladino
    Commented el 7 mar. 2021 a las 17:00

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.