Diseños metodológicos y su aplicación en el campo

Santiago Benitez-Vieyra

Distribución de Poisson

Modelos para conteos

Proceso Poisson

Número de parejas, hijos, semillas, plántulas, granos de polen, etc. Número de individuos observados en un lapso de tiempo. Número de individuos observados en un área determinada.

Registro de sucesos que ocurren en un lapso determinado de tiempo y que son independientes unos de otros.
Aproximación a eventos de naturaleza binomial (de probabilidad muy pequeña) o a sucesos de naturaleza multinomial (categóricos).

Distribución de Poisson

Tiene un único parámetro, \(\mu\) (suele aparecer como \(\lambda\).
La media es igual a la varianza.
Al aumentar la media… aumenta la varianza.

La variable respuesta se distribuye como Poisson \[Y_{i} \sim P(\mu_{i})\] Su valor esperado es \(\mu_{i}\), su varianza es igual a su valor esperado. \[E(Y_{i}) = \mu_{i}\] \[var(Y_{i} = \mu_{i})\] El enlace canónico es logaritmo \[\eta =log(\mu_{i})\] Por lo que el valor esperado es igual al exponencial de la parte sistemática. \[\mu_{i} = e^{\eta}\] Colocando las variables predictoras en la parte sistemática… \[\mu_{i} = e^{\beta_{0} + \beta_{1}X_{1i} + ... + \beta_{k}X_{ki}}\]

dat <- read.table("RoadKills.txt", header = TRUE)
fit <- lm(TOT.N ~ D.PARK, data=dat)
summary(fit)

## 
## Call:
## lm(formula = TOT.N ~ D.PARK, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.470  -9.316  -1.511   6.269  53.441 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  57.3049     4.5479  12.600  < 2e-16 ***
## D.PARK       -2.4763     0.3113  -7.955 1.95e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.29 on 50 degrees of freedom
## Multiple R-squared:  0.5586, Adjusted R-squared:  0.5498 
## F-statistic: 63.28 on 1 and 50 DF,  p-value: 1.952e-10

dat <- read.table("RoadKills.txt", header = TRUE)
fit <- glm(TOT.N ~ D.PARK, data = dat, family = poisson)
summary(fit)

## 
## Call:
## glm(formula = TOT.N ~ D.PARK, family = poisson, data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.1100  -1.6950  -0.4708   1.4206   7.3337  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.316485   0.043220   99.87   <2e-16 ***
## D.PARK      -0.105851   0.004387  -24.13   <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: 1071.4  on 51  degrees of freedom
## Residual deviance:  390.9  on 50  degrees of freedom
## AIC: 634.29
## 
## Number of Fisher Scoring iterations: 4

¿Es significativa la devianza residual?

Esto indicaría un mal ajuste del modelo, si bien es un test vago.

pchisq(390.9, 50, lower.tail=F)

## [1] 2.321894e-54

dat <- read.table("RoadKills.txt", header = TRUE)
fit <- glm(TOT.N ~ D.PARK, data = dat, family = poisson)
anova(fit, test = "Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      51     1071.4              
## D.PARK  1   680.55        50      390.9 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dat <- read.table("RoadKills.txt", header = TRUE)
fit0 <- glm(TOT.N ~ 1, data = dat, family = poisson)
fit1 <- glm(TOT.N ~ D.PARK, data = dat, family = poisson)
anova(fit0, fit1, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: TOT.N ~ 1
## Model 2: TOT.N ~ D.PARK
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        51     1071.4                          
## 2        50      390.9  1   680.55 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La diferencia de devianzas se distribuye como \(\chi^2\) con sus grados de libertad igual a la diferencia de grados de libertad entre los modelos.

Interpretación de parámetros

dat <- read.table("RoadKills.txt", header = TRUE)
fit <- glm(TOT.N ~ D.PARK, data = dat, family = poisson)
fit$coefficients

## (Intercept)      D.PARK 
##   4.3164850  -0.1058506

exp(fit$coefficients)

## (Intercept)      D.PARK 
##   74.924803    0.899559

1/exp(fit$coefficients)

## (Intercept)      D.PARK 
##  0.01334672  1.11165580

El número de muertos disminuye un 11.17% por cada kilómetro de distancia.

Interpretación de parámetros

Si b es cercano a 0, exp(b) es cercano a 1. El efecto es pequeño.
Si b es 0.13976, exp(b) es 1.15. Por cada aumento de una unidad en la variable X, la variable Y se incrementa un 15%.
Si b es -0.2231, exp(b) es 0.80, 1/exp(b) es 1.25. Por cada disminución de una unidad en la variable X, la variable Y disminuye un 25%.

Ver el patrón de serpenteo. Violación de supuestos?

Sobredispersión

La sobredispersión es provocada por variación al azar NO EXPLICADA en la variable respuesta.
Habitualmente se observa una importante devianza residual, que permanece incluso si se incorporan variables al modelo.
Esta variabilidad provoca que la relación media-varianza esperada para la familia no se cumpla.
Debe distinguirse de la sobredispersión aparente , debida a variables o interacciones faltantes, efectos no lineales no considerados, outliers en la variable respuesta o errores en la elección del enlace.

Uno de los supuestos de un MLG es que la variable dependiente se ajusta a una de las funciones de la familia exponencial. Cada una de estas familias está caracterizada por una relación media-varianza específica.

Distribución	Posición	Dispersión
Poisson	\(E(Y) = \mu\)	\(var(Y) = \mu\)
Binomial	\(E(Y) = N\pi\)	\(var(Y) = N\pi(1-\pi)\)

Sin embargo, a veces la variabilidad observada es mayor en una proporción \(\phi\) respecto a lo esperado. Si \(\phi > 1\) entonces existe sobredispersión

Distribución	Posición	Dispersión
Poisson	\(E(Y) = \mu\)	\(var(Y) = \phi\mu\)
Binomial	\(E(Y) = N\pi\)	\(var(Y) = \phi N\pi(1-\pi)\)

La sobredispersión provoca que los errores estándar de los parámetros sean la raíz cuadrada de \(\phi\) veces más chicos que lo que deberían (y aumenta por tanto el error de tipo I)

JAMÁS pueden presentar sobredispersión:

1- Datos binarios (Binomiales no agregados: 0 y 1).
2- Cuando el modelo en consideración es el saturado (como los modelos loglineales para tablas de contingencia).

Si nuestro modelo presenta un buen ajuste y aún así permanece mucha devianza residual no explicada debemos estimar \(\phi\) para probar si existe sobredispersión.

El parámetro \(\phi\) es igual a la suma de los residuos Pearson sobre sus grados de libertad (o a la devianza residual sobre los grados de libertad). Para evitar hacer la cuenta podemos ajustar un modelo de quasiverosimilitud.
En caso de ser significativa la sobre dispersión, nuestro modelo final será el de quasiverosimilitud.

No hay acuerdo en cuál es la magnitud de sobredispersión tolerable. Podemos usar una aproximación práctica (examinar si los valores P cambian significativamente), o usar el límite arbitrario de \(\phi\) > 2.

## 
## Call:
## glm(formula = TOT.N ~ D.PARK, family = poisson, data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.1100  -1.6950  -0.4708   1.4206   7.3337  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.316485   0.043220   99.87   <2e-16 ***
## D.PARK      -0.105851   0.004387  -24.13   <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: 1071.4  on 51  degrees of freedom
## Residual deviance:  390.9  on 50  degrees of freedom
## AIC: 634.29
## 
## Number of Fisher Scoring iterations: 4

## Warning in anova.glm(fit, test = "F"): using F test with a 'poisson' family
## is inappropriate

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
## NULL                      51     1071.4                     
## D.PARK  1   680.55        50      390.9 680.55 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelos de cuasiverosimilitud.

No existen la “familia quasipoisson” o “quasibinomial”, sino que especificamos una relación media-varianza y ajustamos el modelo mediante cuasiverosimilitud.

Ventaja: Controlan la sobredispersión, por lo que disminuyen el error del valor P, los valores estimados de los parámetros (y por lo tanto la interpretación biológica) no cambian.
Desventaja: Dificultan la selección e modelos, se utiliza QAIC en vez de AIC.

Modelos Binomiales Negativos.

Distribución	Posición	Dispersión
Poisson	\(E(Y) = \mu\)	\(var(Y) = \mu\)
Binomial Negativa	\(E(Y) = N\pi\)	\(var(Y) = \mu + \theta\mu^2\)

El parámetro \(\theta\) regula la sobredispersión.

library(MASS)
dat <- read.table("RoadKills.txt", header = TRUE)
nbfit <- glm.nb(TOT.N ~ D.PARK, data=dat)
summary(nbfit)

## 
## Call:
## glm.nb(formula = TOT.N ~ D.PARK, data = dat, init.theta = 3.681040094, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4160  -0.8289  -0.2116   0.4800   2.1346  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.41072    0.15476   28.50   <2e-16 ***
## D.PARK      -0.11612    0.01137  -10.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.681) family taken to be 1)
## 
##     Null deviance: 155.445  on 51  degrees of freedom
## Residual deviance:  54.742  on 50  degrees of freedom
## AIC: 393.09
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.681 
##           Std. Err.:  0.891 
## 
##  2 x log-likelihood:  -387.093

library(MASS)
dat <- read.table("RoadKills.txt", header = TRUE)
nbfit <- glm.nb(TOT.N ~ D.PARK, data=dat)
anova(nbfit, test = "Chisq")

## Warning in anova.negbin(nbfit, test = "Chisq"): tests made without re-
## estimating 'theta'

## Analysis of Deviance Table
## 
## Model: Negative Binomial(3.681), link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      51    155.445              
## D.PARK  1    100.7        50     54.742 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los modelos con la familia binomial negativa pueden estar sobredispersos (examinar la devianza residual). Pero no existen métodos de cuasi-verosimilitud para ellos.
La familia binomial negativa permite calcular el AIC, lo que facilita la elección de modelos.
Los modelos binomial negativos están anidados dentro de un modelo poisson, ya que solamente difieren en un parámetro (\(\theta\)). Puede testearse su adecuación con una cociente de verosimilitudes.

library(MASS)
dat <- read.table("RoadKills.txt", header = TRUE)
nbfit <- glm.nb(TOT.N ~ D.PARK, data = dat)
pofit <- glm(TOT.N ~ D.PARK, data = dat, family = poisson)
L.NB = logLik(nbfit); L.NB

## 'log Lik.' -193.5465 (df=3)

L.Po = logLik(pofit); L.Po

## 'log Lik.' -315.1434 (df=2)

d <- 2 * (L.NB - L.Po)
attributes(d) <- NULL # para borrar resultados innecesarios
d

## [1] 243.1939

pchisq(d, df = 1, lower.tail = FALSE)/2

## [1] 3.956301e-55

Diseños metodológicos y su aplicación en el campo

Santiago Benitez-Vieyra

Distribución de Poisson

Modelos para conteos

Proceso Poisson

Distribución de Poisson

¿Es significativa la devianza residual?

Interpretación de parámetros

Interpretación de parámetros

Sobredispersión

Sobredispersión

JAMÁS pueden presentar sobredispersión:

Modelos de cuasiverosimilitud.

Modelos Binomiales Negativos.

END