Ensayos biológicos, uso de Modelos Lineales Generalizados (MLG)

Ricardo Alves de Olinda

ricardo.estat@yahoo.com.br

http://lattes.cnpq.br/7767223263366578

Universidad del Estado de Paraíba

http://departamentos.uepb.edu.br/estatistica/corpo-docente/

Universidad de San Carlos de Guatemala

Simposio Nacional de Estadística Aplicada 2024

en las ciencias agronómicas, ambientales y forestales. En conmemoración de los 25 años de fundación del CETE

Con el apoyo de:

Use R!

R

Libro sobre modelos lineales generalizados en R Aqui

Nota

Qué es una regresión? Y un ANOVA? Cuál es la principal diferencia entre ambos? Qué supuestos estadísticos debemos asumir cuando llevemos a cabo este tipo de análisis? Estas y otras preguntas son críticas en la aplicación de modelos lineales a la resolución de problemas estadísticos.

Esquema conceptual de los pasos que deben seguirse a la hora de ajustar un modelo lineal univariante. (Luis Cayuela, 2015).

Análisis con datos de incidencia

Una aplicación de la distribución binomial es cuando tenemos datos de incidencia, es decir, presencia o ausencia, de alguna variable. Por ejemplo, presencia o ausencia de una especie o individuo en un lugar. En este caso la fórmula es diferente y el modelo es parecido a una regresión logística, veamos.

Aquí utilizaremos datos sobre la autotomía de la cola de lagartos de la especie Coleodactylus meridionalis observados en fragmentos de bosque de la Mata Atlántica en el estado de Pernambuco (C. N. Oliveira et al. 2020).

# Extensión en .csv  texto separado por espacios
lagartos = read.csv("lagarto.csv",sep=";", header=TRUE) 

# Leyendo las seis primeras líneas
head(lagartos)

##   Numero   Sexo   SVL Intact_tail_length Autotomized_tail_length Tail_state
## 1      2   Male 20.70                 NA                   12.88          0
## 2      3   Male 21.10                 NA                   13.07          0
## 3      6 Female 23.72                 NA                   17.56          0
## 4      9   Male 18.84              17.38                      NA          1
## 5     21   Male 22.20                 NA                   16.50          0
## 6     22   <NA> 20.59                 NA                   12.46          0

Descripción de variables:

1. Sexo
2. SVL: Tamaño del cuerpo del lagarto
3. Intact_tail_length: tamaño de la cola;
4. Autotomized_tail_length: cola autotomizada;
5. Tail_state: estado de cola

Pregunta

¿La probabilidad de que los lagartos de la especie Coleodactylus meridionalis pierdan (autotomicen) su cola aumenta con el tamaño del cuerpo y según el sexo del lagarto?

Predicciones

Cuanto más grande es el lagarto, mayor es la probabilidad de autotomía de la cola y esta respuesta también podría diferir entre sexos debido al dimorfismo sexual.

Variables

Variable respuesta: Presencia o ausencia de cola en lagartos encontrados por búsqueda activa.

Modelo Lineal Generalizado(MLG)

Qué son los MLG?

Los modelos lineales (regresión, ANOVA, ANCOVA), se basan en los siguientes supuestos:

1. Los errores se distribuyen normalmente.

2. La varianza es constante.

3. La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s).

• En muchas ocasiones, sin embargo, nos encontramos con que uno o varios de estos supuestos no se cumplen por la naturaleza de la información.

• Estos problemas se pueden llegar a solucionar mediante la transformación de la ariable respuesta (por ejemplo tomando logaritmos).

• Sin embargo estas transformaciones no siempre consiguen corregir la falta de normalidad, la heterocedasticidad (varianza no constante) o la no linealidad de nuestros datos.

• Además resulta muchas veces interpretar los resultados obtenidos, si utilizamos transformaciones de la variable.

• En el Modelo Lineal Generalizado (GLM) tenemos variables de respuesta asociadas a covariables

• Mientras en el Modelo Lineal combinamos aditividad de los afectos de las covariables con normalidad de las respuestas y homoscedasticidad, en el GML estas tres cosas no se satisfacen necesariamente.

• Los GLM permiten incluir respuestas no normales, como binomial, Poisson, Gamma y otras distibuciones, y en la teoría clásica la estimación se realiza mediante el método de máxima verosimilitud.

Ciertos tipos de variables respuesta sufren invariablemente la violación de estos dos supuestos de los modelos normales y los GLM ofrecen una buena alternativa para tratarlos. Especícamente, podemos considerar utilizar GLM cuando la variable respuesta es:

1. un conteo de casos (p.e. abundancia de una especie);

2. un conteo de casos expresados como proporciones (p.e. porcentaje de plántulas muertas en un experimento de vivero);

3. una respuesta binaria (p.e. vivo o muerto, infectado o no infectado);

Regresando al ejemplo…

Comienza cargando las librerías y preparando los datos a utilizar.

library(ecodados)
library(visdat)
library(tidyverse)
library(lattice)
library(RVAideMemoire)
library(DHARMa)
library(performance)
library(MuMIn)
#library(piecewiseSEM)
library(MASS)
library(ggExtra)
library(Rmisc)
library(emmeans) 
library(sjPlot)
library(bbmle)
library(glmmTMB)
library(ordinal)
library(car)
library(ecolottery)
library(naniar)
library(vcd)
library(generalhoslem)

Análisis de datos exploratorios

A este conjunto de datos le faltan muchas entradas (codificadas como NA). Primero visualizaremos el conjunto de datos y luego debemos eliminar las filas que contienen datos faltantes. Aquí podemos usar la función interna de ggplot2::remove_missing() para eliminar líneas cuyas variables informadas en el argumento faltan

colnames(lagartos) <- c("numero", "sexo", "SVL", "comprimento_cauda", "cauda_autotomizada", "estado_cauda")

## Dados faltantes
vis_miss(lagartos, cluster = TRUE)

Tenga en cuenta que falta el 22,9% de los datos. Eliminemos las filas con datos faltantes para la variable sexo.

## Eliminar datos faltantes
dados_semNA <- remove_missing(lagartos, vars = "sexo") 

## Warning: Removed 84 rows containing missing values or values outside the scale
## range.

vis_miss(dados_semNA)

Ahora, siguiendo lo que estamos acostumbrados a hacer, visualicemos los datos con nuestra hipótesis.

## 
ggplot(dados_semNA, aes(SVL, estado_cauda)) +
    geom_point(aes(shape = sexo, colour = sexo), size = 4, alpha = 0.4) +
    geom_smooth(method = "glm", method.args = list(family = "binomial")) +
    labs(y = "Estado de Cola", x = "Longitud Rostro-Cloacal (mm)", shape = "Sexo", colour = "Sexo") +
    tema_livro()

## `geom_smooth()` using formula = 'y ~ x'

MLG

• Supongamos que observamos las variables de respuesta \(Y_1,\ldots,Y_n\) que son v.a. independientes relacionadas con las covariables \(x_{i1},x_{i2},\ldots,x_{ip}\), \(1 \le i \le n\).

• Genericamente pensemos en una respuesta \(Y_y\) con covariables \(x_{1},x_{2},\ldots,x_{p}\)

Componentes del Modelo

• Podemos pensar que el GLM posee tres componentes:

1. Componente Aleatorio : la variable de respuesta Y tiene distribución

\[ exp\left\{ \frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right\}, \] siendo \(\theta\) el parámetro canónico, \(\phi\) es un parámetro de dispersión y las funciones \(a()\), \(b()\), y \(c()\) son conocidas. Además, se cumple que:

\(\mu=E(Y)=b'(\theta)\) y \(Var(Y)=a(\phi)b''(\theta)\)

2. Componente Sistemático : el vector de covariables \(\mathbf{x'}=(x_1,x_2,\ldots,x_p)\) que da origen al predictor lineal \[ \eta=\sum_{j=1}^px_j\beta_j=\mathbf{x'\beta} \] siendo \(\mathbf{\beta}\) el vector a estimar

3. Función de enlace o link : relaciona los dos componentes \(\mu\) y \(\eta\) \[ g(\mu)=\eta \]

Nota: En algunos casos \(a(\phi)\) es de la forma \(a(\phi)=\frac{\phi}{w}\), siendo w un peso conocido.

Los modelos lineales generalizados permiten dos extensiones:

I. podemos tratar distribuiciones que pertenezcan a una familia exponencial

2. podemos elegir una función de enlace que sea una función monótona y diferenciable.

El Modelo Lineal Generalizado tuvo mucha difusión a partir del libro de McCullagh y Nelder (1989). En estos modelos la variable de respuesta \(Y_i\) sigue una distribución que pertenece a una familia exponencial con media \(\mu_i\) que es una función, por lo general no lineal, de \(\mathbf{x_i'\beta}\).

Ejemplos…

1. Normal: \(Y \sim N(\mu,\sigma)\) \[ f(y,\theta,\phi)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2} \right)\\ =exp\left( \frac{y\mu-\mu^2/2}{\sigma^2}-\frac{1}{2}\left[ \frac{y^2}{\sigma^2}+log(2\pi\sigma^2)\right] \right), \] por lo tanto \(\theta=\mu\), \(b(\theta)=\frac{\mu^2}{2}\), \(\phi=\sigma^2\), \(a(\phi)=\phi\) y \(c(y,\phi)=-\frac{1}{2}\left[\frac{y^2}{\phi}+log(2\pi\phi) \right]\).

\[ E(Y)=\mu \]

En el caso heteroscedástico \(Y\sim N(\mu,\frac{\sigma^2}{w})\), donde w es un peso conocido, tenemos \(\phi=\sigma^2\) y \(a(\phi)=\frac{\phi}{w}\).

2. Caso Binomial: \(Y\sim Bi(n,p)\)

Considerando \(\frac{Y}{n}=\) , la proporción de éxitos.

\[ P\left(\frac{Y}{n}=y \right)=P(Y=ny)=\binom{n}{ny}p^{ny}(1-p)^{n-ny} \] \[ P\left(\frac{Y}{n}=y \right)=exp\left(\frac{ylog\left(\frac{p}{1-p} \right)+log(1-p)}{1/n} +log\binom{n}{ny}\right), \] por lo tanto \(\theta=log\frac{p}{1-p}\) , \(b(\theta)=log(1+e^{\theta})\) , \(\phi=n\) , \(a(\phi)=\frac{1}{\phi}\) y \(c(y,\phi)=\binom{\phi}{\phi y}\) .

\[ E\left(\frac{Y}{n} \right)=p=\frac{e^\theta}{1+e^\theta}=\frac{1}{1+e^{-\theta}} \]

3. Caso Poisson: \(Y\sim P(\lambda)\) \[ P(Y=y)=e^{-\lambda}\frac{\lambda^y}{y!}\\ P(Y=y)=exp(ylog\lambda-\lambda-log(y!)) \] por lo tanto \(\theta=log\lambda\), \(b(\theta)=e^{\theta}\), \(\phi=1\) y \(c(y,\phi)=-log(y!)\) \[ E(Y)=\lambda=e^{\theta} \]

Función de enlace o link:

Esta función relaciona el predictor lineal \(\eta\) con la esperanza \(\mu\) de la respuesta Y. A diferencia del modelo lineal clásico, aquí introducimos una función uno-a-uno continua y diferenciable, \(g(\mu)\), tal que \[ \eta=g(\mu) \] En el caso Binomial, por ejemplo, tenemos que \(\mu \in (0,1)\) y la función de enlace tiene que pertenecer a la recta de los números reales. Suelen usarse 3 funciones de enlaces:

1. Logit: \(\eta=log\frac{\mu}{1-\mu}\)
2. Probit: \(\eta=\Phi^{-1}(\mu)\)
3. Complemento log-log: \(\eta=log(-log(1-\mu))\)

Links Canónicos:

En el caso normal mostramos que si \(Y\sim(\mu,\sigma^2)\) el parámetro canónico es \(\theta=\mu\).

En el caso binomial \(Y\sim Bi(n,p)\) en el que consideramos \(\frac{Y}{n}\) vimos que el vector canónico es \(\theta=logit(\pi)\). Estos son los enlaces más usados en cada caso.

Cuando usamos \(\eta=\theta\) el modelo tiene el enlace canónico convencional. Es conveniente usar el enlace convencional, ya que algunas cosas se simplifican, pero la posibilidad de usarlo dependerá de los datos con los que estemos trabajando.

• Normal: \(\eta\mu\)
• Poisson: \(\eta=log\mu\)
• Binomial: \(\eta=log\frac{\mu}{1-\mu}\)
• Gamma: \(\eta=\mu^{-1}\)

Función de Verosimilitud para el GLM:

Sea Y una v.a. con función de densidad de probabilidad perteneciente a una familia exponencial dada por: \[ f_Y(y,\theta,\phi)=exp\left\{ \frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right\}, \] para algunas funciones conocidas \(a(\phi)\), \(b(\theta)\) y \(c(y , \phi)\). \(S_i\) es un parámetro conocido, esta es una familia exponencial con parámetro canónico natural \(\theta\).

Si \(\phi\) no es conocido, ésta puede ser una familia exponencial en \((\theta, \phi)\) o no. \(\phi\) es un parámetro de dispersión o de forma.

La media \(E(Y)\) es sólo función de \(\theta\) y es por lo tanto el parámetro de interés; \(\phi\) en general es tratado como un parámetro nuisance o de ruido. En la mayoría de los casos \(\phi\) no será tratado tal como es tratado \(\theta\). Estimaremos y haremos inferencia bajo un valor asumido de \(\phi\) y si \(\phi\) necesita ser estimado, lo estimaremos y luego será tomado como un valor fijo y conocido.

Esta familia incluye distribuciones simétricas, asimétricas, discretas y continuas, tales como la distribución Normal, Binomial, Poisson o Gamma.

Momentos de una familia exponencial:

Deduciremos el primer y segundo momento de una familia exponencial a partir del logaritmo de su verosimilitud. \[ L(\theta,y)=\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi). \] Su primera derivada es: \[ L'(\theta,y)=\frac{\partial L(\theta,y)}{\partial\theta}=\frac{y-b'(\theta)}{a(\phi)}, \] mientras que su derivada segunda es: \[ L''(\theta,y)=\frac{\partial L^2(\theta,y)}{\partial^2\theta}=\frac{-b''(\theta)}{a(\phi)}. \] Como \(E\left( \frac{\partial L(\theta,y)}{\partial\theta}\right)=0\), entonces \[ 0=E(L'(\theta,y))=E\left[ \frac{y-b'(\theta)}{a(\phi)}\right] \] y por lo tanto \[ \mu=E(Y)=b'(\theta). \] Además, sabemos que \[ E(L''(\theta,y))=-E[(L'(\theta,y))^2], \] entonces \[ Var(L'(\theta,y))=E[(L'(\theta,y))^2]=-E(L''(\theta,y))=\frac{b''(\theta)}{a(\phi)}. \] Por otro lado, \[ Var(L'(\theta,y))=Var\left(\frac{y-b'(\theta)}{a(\phi)}\right)=\frac{1}{a^2(\phi)}Var(Y) \] y en consecuencia \[ Var(Y)=a(\phi)b''(\theta). \] La varianza es el producto de dos funciones: una que depende del parámetro natural, \(\theta\), y otra que depende sólo del parámetro nuisance \(\phi\). \(V(\theta) = b''(\theta)\) es llamada función de varianza del modelo.

Resumiendo: \[ E(Y)=b'(\theta)\\ Var(Y)=a(\phi)b''(\theta) \] Nota: Estimación de los parámetros: Método de Newton-Raphson y Fisher-scoring

Modelado

Aquí construiremos dos modelos con la misma distribución binomial, pero con dos funciones de enlace: logit y probit. La función logit tiene colas ligeramente más planas, es decir, la curva probit se acerca a los ejes más rápidamente que la función logit. Generalmente no hay mucha diferencia entre ellos. Como no tenemos expectativas sobre qué función de enlace es mejor, podemos hacer una selección de modelo.

## Modelos
mod_log <- glm(estado_cauda ~ SVL * sexo, data = dados_semNA, family = binomial(link = "logit"))
mod_pro <- glm(estado_cauda ~ SVL * sexo, data = dados_semNA, family = binomial(link = "probit"))

# Seleção de modelos
AICctab(mod_log, mod_pro, nobs = 139)

##         dAICc df
## mod_pro 0.0   4 
## mod_log 0.1   4

Hay poca diferencia entre el modelo probit y logit. Como el modelo logit es más sencillo, sólo lo interpretaremos.

Diagnóstico de residuos del modelo.

Diagnosticemos los residuos del modelo.

## Diagnóse avançada
simulationBion <- simulateResiduals(fittedModel = mod_log, plot = T)

Realicemos una visualización de ajuste del modelo

## Ok: Aproximadamente el 100% de los residuos están dentro de los límites de error.
binned_residuals(mod_log)

## Ok: About 100% of the residuals are inside the error bounds.

Inferencia

Para modelos con parámetros de dispersión conocidos (por ejemplo, binomial y Poisson), la chi-cuadrado es la estadística más apropiada.

## Coeficientes estimados pelo modelo
summary(mod_log)

## 
## Call:
## glm(formula = estado_cauda ~ SVL * sexo, family = binomial(link = "logit"), 
##     data = dados_semNA)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    5.5834     2.5909   2.155   0.0312 *
## SVL           -0.2678     0.1178  -2.274   0.0230 *
## sexoMale       0.6977     4.4055   0.158   0.8742  
## SVL:sexoMale  -0.0442     0.2085  -0.212   0.8321  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 191.07  on 138  degrees of freedom
## Residual deviance: 181.38  on 135  degrees of freedom
## AIC: 189.38
## 
## Number of Fisher Scoring iterations: 4

anova(mod_log, test = "Chisq" )

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: estado_cauda
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                       138     191.07            
## SVL       1   9.2563       137     181.82 0.002347 **
## sexo      1   0.3920       136     181.43 0.531262   
## SVL:sexo  1   0.0454       135     181.38 0.831292   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación de resultados

La interpretación de los resultados es que el tamaño corporal (SVL) afecta negativamente la probabilidad de que la cola esté intacta, es decir, al aumentar el tamaño, la probabilidad de que la cola permanezca intacta disminuye. La interacción no fue significativa, por lo que el efecto es independiente del sexo de los lagartos.

Modelos inflados por ceros

En ecología es comum encontrar en datos de conteos valores 0. Sin embargo, sí que es mucho más usual encontrar datos de conteos en los que hay un numero de ceros mayor que el que cabría esperar de acuerdo a una distribución de Poisson o una binomial negativa.

1. Esto puede causar problemas en nuestros modelos ya que, de no tener en cuenta el exceso de ceros:

• Las estimas de los coeficientes pueden ser poco confiables.
• Puede haber sobredispersión.

Hay varios motivos por los quales se da la presecia ceros. Por ejemplo, en el contexto de espécies de pragas en un área forestal:

1. Hay errores estructurales. Es decir, que una plaga no está presente en un parche porque el hábitat no es adecuado.

2. Hay errores de diseño, debidos a un diseño experimental o muestral incorrecto. Por ejemplo, si buscamos a una especie en una época en la que los individuos se encuentran en otro sitio (e.g. espécies de plagas que atacam viveros florestais), es muy probable que nuestros conteos contengan una gran proporcion de ceros.

3. Hay errores de observador. Esto ocurre cuando dos especies son similares y el observador no sabe distinguirlas, o cuando son dificiles de detectar.

Nota

Asumir que los ceros proceden de dos procesos distintos: el proceso binomial y el proceso de Poisson. Igual que en los ZAP, se hace un glm binomial para modelar la probabilidad de medir un 0 (los falsos ceros; por ejemplo cuando el observador no ha visto ningún individuo de la especie en cuestión a pesar de que está presente). Posteriormente se modela la probabilidad de obtener el resto de valores, incluyendo ceros. Estos modelos se llaman Modelos de mescla de distribuciones (mixture models), o modelos ZIP (Zero Inflated) o ZINB.

Ejemplo de aplicación

Los biólogos estatales de vida silvestre quieren modelar cuántos peces capturan los pescadores en un parque estatal. Se pregunta a los visitantes cuánto tiempo estuvieron allí, cuántas personas había en el grupo, si había niños en el grupo y cuántos peces capturaron. Algunos visitantes no pescan, pero no hay datos sobre si una persona pescó o no. Algunos visitantes que pescaron no capturaron ningún pez, por lo que hay un exceso de ceros en los datos debido a las personas que no pescaron.

Como decimos, la regresión de Poisson inflada en cero se utiliza para modelar datos de recuento que tienen un exceso de recuentos de cero. Además, la teoría sugiere que los ceros en exceso se generan mediante un proceso separado de los valores de conteo y que los ceros en exceso se pueden modelar de forma independiente. Así, el modelo ZIP tiene dos partes, un modelo de conteo de Poisson y el modelo logit para predecir el exceso de ceros.

# Extensión en .csv  texto separado por espacios
pescado = read.csv("Pescado.csv",sep=";", header=TRUE) 

# Leyendo las seis primeras líneas
head(pescado)

##   trailer cuant_persona cuant_ninos cuantos
## 1       0             1           0       0
## 2       1             1           0       0
## 3       0             1           0       0
## 4       1             2           1       0
## 5       0             1           0       1
## 6       1             4           2       0

Variables

Tenemos datos de 250 grupos que fueron a un parque. A cada grupo se le preguntó cuántos peces pescaron (cuantos), cuántos niños había en el grupo (cuant_ninos), cuántas personas había en el grupo (cuant_persona) y si trajeron o no un remolque al parque (trailer).

Además de predecir el número de peces capturados, existe interés en predecir la existencia de ceros en exceso, es decir, la probabilidad de que un grupo haya capturado cero peces. Usaremos las variables niño, pearsons y camper en nuestro modelo.

par(mfrow=c(1,1), mar=c(3,2,1,0)+.5, mgp=c(1.6,.6,0), pch=19)
plot(table(pescado$cuantos), xlab = "Número de peces capturados", ylab = "")
grid()

Regresión de Poisson inflada a cero:

Aunque podemos ejecutar una regresión de Poisson en R usando la función glm en uno de los paquetes principales, necesitamos otro paquete para ejecutar el modelo de Poisson inflado en cero. Usamos el paquete pscl.

library(pscl)

## Warning: package 'pscl' was built under R version 4.3.3

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

m1 <- zeroinfl(cuantos ~ cuant_ninos + trailer|cuant_persona, data = pescado)
summary(m1)

## 
## Call:
## zeroinfl(formula = cuantos ~ cuant_ninos + trailer | cuant_persona, data = pescado)
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2369 -0.7540 -0.6080 -0.1921 24.0847 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.59789    0.08554  18.680   <2e-16 ***
## cuant_ninos -1.04284    0.09999 -10.430   <2e-16 ***
## trailer      0.83402    0.09363   8.908   <2e-16 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.2974     0.3739   3.470 0.000520 ***
## cuant_persona  -0.5643     0.1630  -3.463 0.000534 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 12 
## Log-likelihood: -1032 on 5 Df

El resultado se parece mucho al resultado de dos regresiones MCO en R. Debajo de la llamada del modelo, encontrará un bloque de salida que contiene los coeficientes de regresión de Poisson para cada una de las variables, junto con los errores estándar, las puntuaciones z y el valor-p para los coeficientes. A esto le sigue un segundo bloque que corresponde al modelo de inflación. Esto incluye coeficientes logit para predecir el exceso de ceros junto con sus errores estándar, puntuaciones z y valor-p.

Todos los predictores en las partes de conteo e inflación del modelo son estadísticamente significativos. Este modelo se ajusta a los datos significativamente mejor que el modelo nulo, es decir, el modelo de sólo intercepción. Para demostrar que este es el caso, podemos comparar el modelo actual con un modelo nulo sin predictores utilizando la prueba de chi-cuadrado de la diferencia en log-verosimilitudes.

mnull <- update(m1, . ~ 1)

pchisq(2*(logLik(m1) - logLik(mnull)), df = 3, lower.tail = FALSE)

## 'log Lik.' 4.041242e-41 (df=5)

Dado que tenemos tres variables predictoras en el modelo completo, los grados de libertad para la prueba de chi-cuadrado son 3. Esto da como resultado un valor p significativo alto; por lo tanto, nuestro modelo general es estadísticamente significativo.

Tenga en cuenta que el resultado del modelo anterior no indica de ninguna manera si nuestro modelo inflado a cero es una mejora con respecto a una regresión de Poisson estándar. Podemos determinar esto ejecutando el modelo estándar de Poisson correspondiente y luego ejecutando una prueba de Vuong de los dos modelos.

Definición…:

1. Caso Poisson inflada por ceros: \(Y\sim ZIP(\mu_i,\pi_i)\)
   
   

Supongamos que \(y_i,\ldots,y_n\) son Observaciones de la variable respuesta \(Y_i,i=1,2,\ldots,n\), el modelo de Poisson inflado ceros se expresa como sigue,

\[ P(Y_i=y_i|x_i) =\left\{ \begin{array}{cccccccccccccccc} \{\pi_i+(1-\pi_i)e^{-\mu_i}, \;\; y_i=0 \\ (1-\pi_i)\frac{e^{-\mu_i}\mu_i^{y_i}} {y_i!}, \;\; y_i=1,2,3,\ldots,\\ \end{array} \right. \] por lo tanto \(0<\pi_i<1\) e \(\mu_i>0\).

La media y la varianza de la distribución son: \[ E[Y_i] = (1-\pi_i)\mu_i\\ Var[Y_i] = \mu_i(1-\pi_i)(1+\pi_i\mu_i) \]

El modelo ZIP, implica que la media \(\mu_i\) de una variable de Poisson mediante un modelo de regresión de Poisson y la probabilidad de \(\pi_i\) a través de una función de enlace de regresión logística \(\eta_i = logit(\pi_i)\).

2. Caso Binomial Negativa inflada por ceros: \(Y\sim ZINB(\mu_i,\alpha,\pi_i)\)

El modelo Binomial Negativo inflado por ceros é expresado asi,
\[ P(Y_i=y_i|x_i) = \left\{ \begin{array}{ccccccc} \pi_i +(1-\pi_i)\Biggl(\frac{1}{1+\alpha\mu}\Biggr)^{\alpha^{-1}} \;\; y_i=0 \\ (1-\pi_i)\frac{\Gamma(y_i+\alpha^{-1})}{y_i!\Gamma(\alpha^{-1})}\Biggl(\frac{\alpha\mu_i}{1+\alpha\mu_i}\Biggr)^{y_i}\Biggl(\frac{1}{1+\alpha\mu_i}\Biggr)^{\alpha^{-1}} \;\; y_i=1,2,\ldots\\ \end{array} \right. \] por lo tanto \(0<\pi_i<1\) , \(\mu_i>0\).

La media y la varianza de la distribución son:
\[ E[Y_i] = (1-\pi_i)\mu_i\\ Var[Y_i] = \mu_i(1-\pi_i)(1+\pi_i\mu_i+\alpha\mu_i) \]

Regresando ao ejemplo:

summary(p1 <- glm(cuantos ~ cuant_ninos + trailer, family = poisson, 
                  data = pescado))

## 
## Call:
## glm(formula = cuantos ~ cuant_ninos + trailer, family = poisson, 
##     data = pescado)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.91026    0.08119   11.21   <2e-16 ***
## cuant_ninos -1.23476    0.08029  -15.38   <2e-16 ***
## trailer      1.05267    0.08871   11.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2958.4  on 249  degrees of freedom
## Residual deviance: 2380.1  on 247  degrees of freedom
## AIC: 2723.2
## 
## Number of Fisher Scoring iterations: 6

La prueba de Vuong compara el modelo inflado en cero con un modelo de regresión de Poisson ordinario. En este ejemplo, podemos ver que nuestra estadística de prueba es significativa, lo que indica que el modelo inflado a cero es superior al modelo estándar de Poisson.

Existen algunas críticas y debates en la literatura sobre el uso y mal uso de la prueba de Vuong para modelos no anidados en la prueba de inflación cero en modelos de conteo. El objetivo aquí se centra más en cómo implementar varios análisis en R. Consulte las referencias a continuación para obtener más información:

Desmarais Bruce A., Harden Jeffrey J., 2013. Testing for Zero Inflation in Count Models: Bias Correction for the Vuong Test. Stata Journal, 13, 4, 810-835.

Wilson P., 2015. The misuse of the Vuong test for non-nested models to test for zero-inflation. Economics Letters, 127, 51-53.

Podemos obtener intervalos de confianza para los parámetros y parámetros exponenciados de varias maneras, aquí usando bootstrapping. Para el modelo de Poisson, estos serían los índices de riesgo de incidentes; para el modelo de inflación cero, los odds ratios.

Usamos el paquete de arranque. Primero, obtenemos los coeficientes de nuestro modelo original para usarlos como valores iniciales en el modelo y acelerar el tiempo de estimación. A continuación, escribimos una función corta que toma datos e índices como entrada y devuelve los parámetros que nos interesan. Finalmente, pasamos esto a la función de arranque y hacemos 1200 replicaciones. Para obtener resultados finales, es posible que desee aumentar el número de replicaciones para ayudar a garantizar resultados estables.

coef(m1, "count")

## (Intercept) cuant_ninos     trailer 
##   1.5978883  -1.0428391   0.8340236

coef(m1, "zero")

##   (Intercept) cuant_persona 
##     1.2974403    -0.5643474

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:car':
## 
##     logit

## The following object is masked from 'package:lattice':
## 
##     melanoma

 f <- function(data, i) {
     m <- zeroinfl(cuantos ~ cuant_ninos +         trailer|cuant_persona, data = data[i, ],
                   start = list(count = c(1.598, -1.0428,   0.834), zero = c(1.297, -0.564)))
     as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
  }
 set.seed(10)
 res <- boot(pescado, f, R = 1200)
  res

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pescado, statistic = f, R = 1200)
## 
## 
## Bootstrap Statistics :
##         original        bias    std. error
## t1*   1.59788854 -0.0451465197  0.29947818
## t2*   0.08553816  0.0036239299  0.01659791
## t3*  -1.04283849  0.0112733726  0.41120836
## t4*   0.09998829  0.0042049207  0.01580662
## t5*   0.83402218  0.0003032759  0.40949791
## t6*   0.09362679  0.0041210986  0.01524233
## t7*   1.29743916  0.0211989505  0.47396484
## t8*   0.37385225  0.0069068079  0.03657407
## t9*  -0.56434716 -0.0265559374  0.26840571
## t10*  0.16296381  0.0039579583  0.03019734

Los resultados son estimaciones de parámetros alternos y errores estándar. Es decir, la primera línea contiene la estimación del primer parámetro de nuestro modelo. El segundo tiene el error estándar del primer parámetro. La tercera columna contiene los errores estándar bootstrap, que son considerablemente mayores que los estimados por zeroinfl.

Ahora podemos obtener los intervalos de confianza para todos los parámetros. Comenzamos en la escala original con intervalos de confianza (IC) ajustados por percentiles y sesgos. También comparamos estos resultados con intervalos de confianza regulares basados en errores estándar.

## Estimaciones de parámetros básicos con IC ajustados por percentiles

parms <- t(sapply(c(1, 3, 5, 7, 9), function(i) {
     out <- boot.ci(res, index = c(i, i + 1), type = c("perc", "bca"))
     with(out, c(Est = t0, pLL = percent[4], pUL = percent[5],
                 bcaLL = bca[4], bcaLL = bca[5]))
  }))
 ## agregando nombres de línea
 row.names(parms) <- names(coef(m1))
 ## mostrando los resultados
 parms

##                           Est          pLL        pUL        bcaLL       bcaLL
## count_(Intercept)   1.5978885  0.929548809  2.0891957  1.052591464  2.21402500
## count_cuant_ninos  -1.0428385 -1.778733547 -0.1686651 -1.727480382 -0.09541710
## count_trailer       0.8340222  0.008795793  1.5978058  0.001069699  1.58913696
## zero_(Intercept)    1.2974392  0.458982712  2.2596422  0.441060605  2.21684371
## zero_cuant_persona -0.5643472 -1.100952918 -0.1167136 -1.061884331 -0.03569355

## comparar con la aproximación normal
 confint(m1)

##                         2.5 %     97.5 %
## count_(Intercept)   1.4302366  1.7655400
## count_cuant_ninos  -1.2388125 -0.8468657
## count_trailer       0.6505185  1.0175288
## zero_(Intercept)    0.5647033  2.0301772
## zero_cuant_persona -0.8837505 -0.2449442

Los intervalos de confianza bootstrap son considerablemente mayores que la aproximación normal. Ahora podemos estimar el índice de riesgo de incidentes (TIR) para el modelo de Poisson y el odds ratio (OR) para el modelo logístico (inflación cero). Esto se hace usando un código casi idéntico al anterior, pero pasando una función de transformación al argumento h de boot.ci, en este caso exp para exponenciar.

## 
 expparms <- t(sapply(c(1, 3, 5, 7, 9), function(i) {
    out <- boot.ci(res, index = c(i, i + 1), type = c("perc", "bca"), h = exp)
    with(out, c(Est = t0, pLL = percent[4], pUL = percent[5],
                bcaLL = bca[4], bcaLL = bca[5]))
 }))
 ## 
 row.names(expparms) <- names(coef(m1))
 ## 
 expparms

##                          Est       pLL       pUL     bcaLL     bcaLL
## count_(Intercept)  4.9425853 2.5333659 8.0784160 2.8650662 9.1524811
## count_cuant_ninos  0.3524528 0.1688519 0.8447918 0.1777317 0.9089937
## count_trailer      2.3025614 1.0088346 4.9421772 1.0010703 4.8995186
## zero_(Intercept)   3.6599122 1.5824634 9.5796696 1.5543549 9.1783157
## zero_cuant_persona 0.5687313 0.3325540 0.8898401 0.3458036 0.9649360

Para comprender mejor nuestro modelo, podemos calcular la cantidad esperada de peces capturados para diferentes combinaciones de nuestros predictores. De hecho, dado que estamos trabajando con predictores esencialmente categóricos, podemos calcular los valores esperados para todas las combinaciones usando la función expand.grid para crear todas las combinaciones y luego la función de predicción para hacer esto. También eliminamos todas las filas donde el número de niños excede el número de personas, lo que no tiene sentido lógico, usando la función de subconjunto. Finalmente creamos un gráfico.

 newdata1 <- expand.grid(0:3, factor(0:1), 1:4)
 colnames(newdata1) <- c("cuant_ninos", "trailer", "cuant_persona")
 newdata1 <- subset(newdata1, subset=(cuant_ninos<=cuant_persona))
 newdata1$phat <- predict(m1, newdata1)
 library(ggplot2)
 ggplot(newdata1, aes(x = cuant_ninos, y = phat, colour = factor(cuant_persona))) +
     geom_point() +
     geom_line() +
     facet_wrap(~trailer) +
     labs(x = "Número de niños", y = "Previsión de pescado capturado")

R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Material sobre modelos lineales generalizados y estadística espacialAqui

Muchas Gracias!