Coeficientes de retrotransformación de un GLMM Gamma-log

Question

Estoy analizando los datos de un experimento de exclusión, es decir, que durante varios años se mantuvo a las cabras fuera de una valla y dentro de la valla las plantas podían crecer sin ser pastoreadas. Fuera de la valla, el pastoreo continuó como de costumbre.

Ahora ajusté un modelo GLM que contiene la biomasa en toneladas/ha por sitio durante varios años. Como la biomasa no debe ser igual o menor que cero, elegí la distribución Gamma. Requería el enlace logarítmico porque los valores de biomasa diferían mucho - no utilizar el enlace logarítmico daría lugar a residuos extraños.

Elaboré un modelo, predije el resultado y lo representé gráficamente. Los BSKAN y BSKAS de la parte inferior izquierda están fuera de una valla, mientras que los demás están dentro de una zona vallada (exclusión, Z=Zaun=valla, S=Sur, etc.).

Al observar el resumen del modelo, encontramos los coeficientes en la escala "vínculo".

Call:
glm(formula = biomass ~ date.2k * PNR, family = Gamma(link = log), 
    data = data.TRL, control = glm.control(maxit = 100))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.17969  -0.06634  -0.01269   0.05498   0.21929  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.84997    0.07475  11.370 2.26e-11 ***
date.2k          -0.07400    0.02073  -3.569 0.001484 ** 
PNRBSKAS         -0.58319    0.10572  -5.516 9.86e-06 ***
PNRBSKZN         -0.33430    0.10572  -3.162 0.004076 ** 
PNRBSKZS         -1.17416    0.10572 -11.106 3.70e-11 ***
PNRBSKZW         -0.79179    0.10572  -7.490 7.65e-08 ***
date.2k:PNRBSKAS  0.03169    0.02932   1.081 0.290129    
date.2k:PNRBSKZN  0.13376    0.02932   4.562 0.000116 ***
date.2k:PNRBSKZS  0.14976    0.02932   5.107 2.82e-05 ***
date.2k:PNRBSKZW  0.10626    0.02932   3.624 0.001292 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.01203615)

Null deviance: 4.09884  on 34  degrees of freedom
Residual deviance: 0.29665  on 25  degrees of freedom
AIC: -24.6

Number of Fisher Scoring iterations: 4

Pero, como quiero presentar la figura con los coeficientes correspondientes, necesito retrotransformar los coeficientes a la escala de respuesta original. Sé que para un enlace logarítmico esto se hace normalmente utilizando exp(coeficiente) .

Por ejemplo, para el intercepto de BSKAN ("Intercept" en la lista de coeficientes), utilizo:

 exp(0.84997) = 2.327... t/ha

La segunda intercepción se puede encontrar cuando se utiliza:

 exp(0.84997 - 0.58319) = 1.3057 t/ha

Hasta aquí todo bien. Ahora, mirando las pendientes, la cosa se pone problemática. Por supuesto, BSKAN tiene una pendiente negativa... pero tomando un exponente de un número negativo lo hace positivo...

exp(-0.074) = 0.928

No incluir el menos lleva a ...

exp(0.074)= 1.0768

Tal vez (?) esto sea correcto ...

-exp(0.074) = -1.0768 

Pero cuando intento retrotransformar la segunda pendiente me vuelvo loco...

exp(-0.074 + 0.0317) = 0.9585

o...

-exp(0.074 + 0.0317) = -1.1114

Espero que te des cuenta de que sólo estoy haciendo conjeturas. Busqué en todas las fuentes que conocía, indagué en el código r, pero sencillamente no había respuesta. El visreg con el que hice la figura consiguió transformar correctamente las pendientes. Pero quizás... ¿Debería sólo retrotransformar los interceptos y la pendiente de .0.074 es correcta?

¡H-E-L-P!

Answer 1

Dado que la respuesta ajustada no es lineal, la pendiente variará constantemente. visreg habrá producido las líneas ajustadas prediciendo a partir del modelo para un rango de valores sobre date.2k para cada nivel de PNRB .

En general exp(coef) así que exp(-0.74) es derecha y -exp(0.74) es equivocado . Para los demás grupos se suman los coeficientes de ese grupo (por lo que el date.2k coef, de 0.074 Además 0.03169 ) y luego exponenciarlo: exp(-0.074 + 0.03169) .

Es la forma en que estás interpretando estos valores como "pendientes" lo que te está derrotando. No son pendientes, desde luego no en el sentido en que usted parece utilizarlas. Son factores multiplicativos. En la escala logarítmica, el modelo es (se supone que es) lineal o aditivo. En la escala original, los términos aditivos se convierten en multiplicativos.

Por lo tanto, los valores calculados para BSKAS como

> exp(-0.074 + 0.0317)
[1] 0.9585822

indica que para un aumento unitario de date.2k la respuesta es multiplicado por ~ 0,96. Si se multiplica algo por este valor se obtiene más pequeño . En qué medida depende de cuál era el valor de la respuesta al valor de $x$ antes de aumentarlo en 1 unidad.

Ejemplo práctico

No disponemos del modelo ajustado (lo que lo habría simplificado), pero tenemos información suficiente para aproximarnos a él a partir de sus resultados. Ilustremos cómo funciona esto para su modelo, prediciendo $y$ para dos puntos de datos $x$ 2 y 3 (es decir date.2k = 2 y date.2k = 3 ) para el nivel BSKAS . Aquí hay una función (ineficiente) para hacer esto, pero al menos no necesito crear una matriz modelo...

pBSKAS <- function(x, backtf = FALSE) {
  ret <- 0.84997 + (-0.074 * x) + (-0.58319 * 1) + (0.0317 * x)
  if (backtf)
    ret <- exp(ret)
  ret
}

Esto da para $x \in {2, 3}$

> pBSKAS(c(2,3), backtf = TRUE)
[1] 1.199830 1.150136

Si ahora calculamos la magnitud del cambio entre las dos observaciones

> 1.150136 / 1.199830
[1] 0.9585825

vemos el mismo valor que hemos calculado antes exponenciando los coeficientes:

> exp(-0.074 + 0.0317)
[1] 0.9585822

(es ligeramente diferente pero eso será sólo el redondeo en los otros coefs.) Tenga en cuenta que todo el intercepto y el coeficiente de PNRBSKAS al cálculo es ponerle una escala (es decir, desplazarlo hacia arriba y hacia abajo).