En este ejercicio se utilizan datos de un estudio sobre cangrejos (J. Brockmann, Ethology 1996). Cada cangrejo hembra tenía un cangrejo macho adherido en su nido. El estudio investigó los factores que influyen en el hecho de que el cangrejo hembra tuviera otros machos (satelites) viviendo cerca de ella, los factores que se incluyeron el estudio fueron:
Y la variable respuesta es el número de satélites de la hembra (Sa)
Los datos del estudio pueden ser descargados aquí
crab <- read.table("crab.txt")
names(crab) <- c("Obs", "C", "S", "W", "Wt", "Sa")
crab = crab[, -1]
crab$C <- as.factor(crab$C)
crab$S <- as.factor(crab$S)
Utilizamos como guía para este ejercicio, lo presentado en un curso en línea de Penn State University. La idea básica del modelo es modelar datos de conteo que se distribuyen siguiendo una distribución Poisson (Schwartz, 2013):
\[ Y_i\sim Poisson(\mu_i) \] \[ \theta_i=\ln(\mu_i) \] \[ \theta_i=\beta_0+\beta_1X_{i1}+\beta_2 X_{i2}+...+\varepsilon \]
La idea del ejercicio es identificar las variables que influyen en el número de satélites. Veamos entonces como se distribuye la variable respuesta con respecto a las variables regresoras:
par(mfrow = c(2, 2))
plot(crab$C, crab$Sa, pch = 16, col = 2, xlab = "Color", ylab = "Satelites")
plot(crab$S, crab$Sa, pch = 16, col = 3, xlab = "Espina", ylab = "Satelites")
plot(crab$W, crab$Sa, pch = 16, col = 4, xlab = "Caparazón", ylab = "Satelites")
plot(crab$Wt, crab$Sa, pch = 16, col = 6, xlab = "Peso", ylab = "Satelites")
A primera vista en parece que tanto la condición de la espina como el peso de la hembra son las que influyen, sin embargo probaremos primero el modelo completo, para poder determinar estas cuestiones:
modelo1 = glm(Sa ~ C + S + W + Wt, family = poisson(), data = crab)
summary(modelo1)
##
## Call:
## glm(formula = Sa ~ C + S + W + Wt, family = poisson(), data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.029 -1.863 -0.599 0.933 4.945
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3572 0.9670 -0.37 0.7118
## C2 -0.2649 0.1681 -1.58 0.1151
## C3 -0.5137 0.1954 -2.63 0.0085 **
## C4 -0.5313 0.2269 -2.34 0.0192 *
## S2 -0.1504 0.2136 -0.70 0.4812
## S3 0.0874 0.1199 0.73 0.4660
## W 0.0165 0.0489 0.34 0.7358
## Wt 0.4971 0.1663 2.99 0.0028 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 632.79 on 172 degrees of freedom
## Residual deviance: 549.56 on 165 degrees of freedom
## AIC: 920.9
##
## Number of Fisher Scoring iterations: 6
Notemos que dado que las variables Color y Condición de la espina son variables categóricas, se definieron como factores, y automáticamente el paquete las introduce al modelo como variables dummies. Por lo cual se el modelo presenta estimaciones de los coeficientes para cada una de estas, para poder observar como una sola variable podemos echar mano de la función anova y dado que los estimadores son de máxima verosimilitud podemos utilizar como estadístico de prueba una \( X2 \):
anova(modelo1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Sa
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 172 633
## C 3 23.7 169 609 3.0e-05 ***
## S 2 6.9 167 602 0.0310 *
## W 1 43.6 166 559 4.1e-11 ***
## Wt 1 9.1 165 550 0.0026 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Notemos nuevamente que esta función evalua la aportación de la suma de los errores conforme al orden de inclusión de la variable en el modelo, una suma de cuadrados tipo I, contrario a la más comunmente usada suma de cuadrados tipo III, una alternativa a esta situación es utilizar la función drop1:
drop1(modelo1, ~., test = "Chisq")
## Single term deletions
##
## Model:
## Sa ~ C + S + W + Wt
## Df Deviance AIC LRT Pr(>Chi)
## <none> 550 921
## C 3 559 924 9.25 0.0262 *
## S 2 551 919 1.80 0.4069
## W 1 550 919 0.11 0.7362
## Wt 1 559 928 9.07 0.0026 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conbase a estos últimos resultados podemos concluir que son el color y el peso de la hembra las variables que influyen significativamente en el número de satélites. Sin embarego hay que recordar que la distribución Poisson tiene un uníco parámetro \( \lambda \) que es tanto su esperanza como su varianza. Por lo cual dentro del modelo se espera que la media sea igual a la esperanza. Podemos probar esto mediante una razón, pero hagamoslo con el nuevo modelo reducido al que hemos llegado:
modelo_p <- glm(Sa ~ Wt + C, family = poisson, data = crab)
summary(modelo_p)
##
## Call:
## glm(formula = Sa ~ Wt + C, family = poisson, data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.983 -1.927 -0.555 0.864 4.827
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0496 0.2331 -0.21 0.831
## Wt 0.5461 0.0681 8.02 1.1e-15 ***
## C2 -0.2051 0.1537 -1.33 0.182
## C3 -0.4497 0.1757 -2.56 0.011 *
## C4 -0.4523 0.2084 -2.17 0.030 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 632.79 on 172 degrees of freedom
## Residual deviance: 551.78 on 168 degrees of freedom
## AIC: 917.1
##
## Number of Fisher Scoring iterations: 6
Para establecer la razon podemos dividir las desviaciones residuales entre sus grados de libertad \( \frac{551.78}{168}=3.2844 \). Dado que el valor es mayor a 1 podemos interpretar como que el valor de la varianza es tres veces mayor al de la media, y por lo tanto no se cumple el supuesto del modelo de Poisson. Para solucionar esto podemos utilizar algunos modelos más flexibles en cuanto a esto: un modelo cuasi-Poisson o un modelo binomial negativo.
Another way of dealing with over-dispersion is to use the mean regression function and the variance function from the Poisson GLM but to leave the dispersion parameter \( \theta \) unrestricted. Thus\( \theta \), is not assumed to be fixed at 1 but is estimated from the data. This strategy leadsto the same coefficient estimates as the standard Poisson model but inference is adjusted forover-dispersion. Consequently, both models (quasi-Poisson and sandwich-adjusted Poisson)adopt the estimating function view of the Poisson model and do not correspond to models with fully specified likelihoods.
Es decir que la desventaja de estos modelos es que sus estimadores no son de máxima verosimilitud. Veamos el ajuste con las variables ya seleccionadas:
modelo_cp <- glm(Sa ~ Wt + C, family = quasipoisson, data = crab)
summary(modelo_cp)
##
## Call:
## glm(formula = Sa ~ Wt + C, family = quasipoisson, data = crab)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.983 -1.927 -0.555 0.864 4.827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0496 0.4159 -0.12 0.91
## Wt 0.5461 0.1215 4.50 1.3e-05 ***
## C2 -0.2051 0.2743 -0.75 0.46
## C3 -0.4497 0.3136 -1.43 0.15
## C4 -0.4523 0.3719 -1.22 0.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.183)
##
## Null deviance: 632.79 on 172 degrees of freedom
## Residual deviance: 551.78 on 168 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6
anova(modelo_cp, test = "Chisq")
## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: Sa
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 172 633
## Wt 1 71.9 171 561 2e-06 ***
## C 3 9.1 168 552 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Notemos que en este caso los coefficientes no cambian, sin embargo aplicando una prueba \( X^2 \) resulta que el color ya no es significativo.
…modeling over-dispersed count data is to assume a negative binomial (NB) distribution for \( y_i|x_i \) which can arise as a gamma mixture of Poisson distributions. One parameterization of its probability density function is: \[ f(y;\mu,\Theta)=\frac{\Gamma (y+\theta)}{\Gamma (\theta)\cdot y!}\cdot \frac{\mu^y\cdot \theta^\theta}{(\mu +\theta)^{y+\theta}} \] with mean \( \mu \) and shape parameter \( \theta;\Gamma(\cdot) \) is the gamma function. For every fixed \( \theta \), this is of type (1) and thus is another special case of the GLM framework. It also has \( \phi= 1 \) but with variance function \( V(\mu)=\mu +\frac{\mu^2}{\theta} \)
Ajustamos este modelo mediante el paquete MASS:
library(MASS)
modelo_bn <- glm.nb(Sa ~ Wt + C, data = crab)
summary(modelo_bn)
##
## Call:
## glm.nb(formula = Sa ~ Wt + C, data = crab, init.theta = 0.9596988982,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.868 -1.368 -0.317 0.426 2.279
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.425 0.538 -0.79 0.43
## Wt 0.712 0.161 4.41 1e-05 ***
## C2 -0.253 0.349 -0.73 0.47
## C3 -0.522 0.380 -1.37 0.17
## C4 -0.481 0.428 -1.12 0.26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9597) family taken to be 1)
##
## Null deviance: 220.02 on 172 degrees of freedom
## Residual deviance: 196.57 on 168 degrees of freedom
## AIC: 757.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.960
## Std. Err.: 0.175
##
## 2 x log-likelihood: -745.934
anova(modelo_bn, test = "Chisq")
## Warning: tests made without re-estimating 'theta'
## Analysis of Deviance Table
##
## Model: Negative Binomial(0.9597), link: log
##
## Response: Sa
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 172 220
## Wt 1 20.71 171 199 5.3e-06 ***
## C 3 2.74 168 197 0.43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Aquí los coeficientes de los parámetros si cambian, adicionalmente hay alertas sobre theta, y nuevamente el color resulta no significativo.
Sin embargo todaslas estimaciones son cercanas entre si, entonces la pregunta es ¿para qué tanto rollo?