La regresión binomial negativa es similar en aplicación a la regresión de Poisson, pero permite una dispersión excesiva en la variable de conteo dependiente.
Este ejemplo utilizará la función glm.nb en el paquete MASS. La función Anova en el paquete car se usará para un análisis de desviación, y la función nagelkerke se usará para determinar un valor p y un valor pseudo R cuadrado para el modelo. El análisis post-hoc se puede realizar con el paquete emmeans.
Tenga en cuenta que los supuestos del modelo y las trampas de este enfoque no se analizan aquí. Se insta al lector a comprender los supuestos de este tipo de modelado antes de continuar.
Input = ("
Garden Monarchs
A 0
A 4
A 2
A 2
A 0
A 6
A 0
A 0
B 5
B 9
B 7
B 5
B 7
B 5
B 9
B 5
C 10
C 14
C 12
C 12
C 10
C 16
C 10
C 10
")
Data = read.table(textConnection(Input),header=TRUE)
Data## Garden Monarchs
## 1 A 0
## 2 A 4
## 3 A 2
## 4 A 2
## 5 A 0
## 6 A 6
## 7 A 0
## 8 A 0
## 9 B 5
## 10 B 9
## 11 B 7
## 12 B 5
## 13 B 7
## 14 B 5
## 15 B 9
## 16 B 5
## 17 C 10
## 18 C 14
## 19 C 12
## 20 C 12
## 21 C 10
## 22 C 16
## 23 C 10
## 24 C 10Data$Garden = factor(Data$Garden,
levels=unique(Data$Garden))
str(Data)## 'data.frame': 24 obs. of 2 variables:
## $ Garden : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 2 2 ...
## $ Monarchs: int 0 4 2 2 0 6 0 0 5 9 ...library(MASS)
model.nb = glm.nb(Monarchs ~ Garden,
data=Data,
control = glm.control(maxit=1000))## Warning in sqrt(1/i): Se han producido NaNs
## Warning in sqrt(1/i): Se han producido NaNslibrary(car)## Loading required package: carDataAnova(model.nb,
type="II",
test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## LR Chisq Df Pr(>Chisq)
## Garden 66.464 2 3.694e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(rcompanion)
nagelkerke(model.nb)## $Models
##
## Model: "glm.nb, Monarchs ~ Garden, Data, glm.control(maxit = 1000), 510466474900, log"
## Null: "glm.nb, Monarchs ~ 1, Data, glm.control(maxit = 1000), 1.749060129, log"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.255141
## Cox and Snell (ML) 0.776007
## Nagelkerke (Cragg and Uhler) 0.778217
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -2 -17.954 35.907 1.5952e-08
##
## $Number.of.observations
##
## Model: 24
## Null: 24
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"library(multcompView)
library(emmeans)
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geysermarginal = emmeans(model.nb,
~ Garden)
pairs(marginal,
adjust="tukey")## contrast estimate SE df z.ratio p.value
## A - B -1.312 0.301 Inf -4.358 <.0001
## A - C -1.904 0.286 Inf -6.647 <.0001
## B - C -0.592 0.173 Inf -3.426 0.0018
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatescld(marginal,
alpha = 0.05,
Letters = letters, ### Use lower-case letters for .group
type = "response", ### Report emmeans in orginal scale
adjust = "tukey") ### Tukey adjustment for multiple comparisons## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## Garden response SE df asymp.LCL asymp.UCL .group
## A 1.75 0.468 Inf 0.924 3.31 a
## B 6.50 0.901 Inf 4.668 9.05 b
## C 11.75 1.212 Inf 9.185 15.03 c
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05La regresión inflada con cero es similar en aplicación a la regresión de Poisson, pero permite una abundancia de ceros en la variable de conteo dependiente.
Este ejemplo utilizará la función zeroinfl en el paquete pscl. La función Anova en el paquete car se usará para un análisis de desviación, y la función nagelkerke se usará para determinar un valor p y un valor pseudo R cuadrado para el modelo. El análisis post-hoc se puede realizar con el paquete emmeans.
library(pscl)## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileismodel.zi = zeroinfl(Monarchs ~ Garden,
data = Data,
dist = "poisson")
### dist = "negbin" may be used
summary(model.zi)##
## Call:
## zeroinfl(formula = Monarchs ~ Garden, data = Data, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.8156 -0.5883 -0.2188 0.3112 1.9807
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2182 0.2847 4.278 1.89e-05 ***
## GardenB 0.6536 0.3167 2.064 0.039 *
## GardenC 1.2457 0.3029 4.113 3.90e-05 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.046e-02 7.363e-01 -0.096 0.924
## GardenB -2.057e+01 1.071e+04 -0.002 0.998
## GardenC -2.057e+01 1.071e+04 -0.002 0.998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 9
## Log-likelihood: -48.14 on 6 Dflibrary(car)
Anova(model.zi,
type="II",
test="Chisq")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## Df Chisq Pr(>Chisq)
## Garden 2 23.914 6.414e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(rcompanion)
nagelkerke(model.zi)## $Models
##
## Model: "zeroinfl, Monarchs ~ Garden, Data, poisson"
## Null: "zeroinfl, Monarchs ~ 1, Data, poisson"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.284636
## Cox and Snell (ML) 0.797356
## Nagelkerke (Cragg and Uhler) 0.800291
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -4 -19.156 38.311 9.6649e-08
##
## $Number.of.observations
##
## Model: 24
## Null: 24
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"library(multcompView)
library(emmeans)
marginal = emmeans(model.zi,
~ Garden)
pairs(marginal,
adjust="tukey")## contrast estimate SE df z.ratio p.value
## A - B -4.75 1.18 Inf -4.032 0.0002
## A - C -10.00 1.43 Inf -6.994 <.0001
## B - C -5.25 1.51 Inf -3.476 0.0015
##
## P value adjustment: tukey method for comparing a family of 3 estimatescld(marginal,
alpha=0.05,
Letters=letters, ### Use lower-case letters for .group
adjust="tukey") ### Tukey adjustment for multiple comparisons## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## Garden emmean SE df asymp.LCL asymp.UCL .group
## A 1.75 0.759 Inf -0.0614 3.56 a
## B 6.50 0.901 Inf 4.3477 8.65 b
## C 11.75 1.212 Inf 8.8563 14.64 c
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05La regresión robusta de Poisson es robusta a los valores atípicos en la variable dependiente.
Este ejemplo usa la función glmRob en el paquete robust. La función anova se puede utilizar para realizar un análisis de desviación. El valor p del modelo se puede encontrar comparando el modelo con un modelo nulo. Sin embargo, en el momento de escribir este artículo, no conozco ninguna forma de determinar el AIC o la pseudo R cuadrado para el modelo.
En el momento de redactar este artículo, la función glmRob solo puede utilizar las familias de modelos de Poisson y binomial.
library(robust)## Loading required package: fit.models## Registered S3 methods overwritten by 'robust':
## method from
## plot.covfm fit.models
## print.covfm fit.models
## summary.covfm fit.models
## print.summary.covfm fit.models## rlm is already registered in the fit.models registry## covfm is already registered in the fit.models registry##
## Attaching package: 'robust'## The following objects are masked from 'package:fit.models':
##
## ddPlot.covfm, distancePlot.covfm, ellipsesPlot.covfm,
## screePlot.covfmmodel.rob = glmRob(Monarchs ~ Garden,
data = Data,
family = "poisson")
anova(model.rob, test="Chisq")## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL NA NA 23 430.19850 NA
## Garden 2 400.8615 21 29.33699 4.906414e-63model.rob.null = glmRob(Monarchs ~ 1,
data = Data,
family = "poisson")
anova(model.rob.null, model.rob, test="Chisq")## Terms Resid. Df Resid. Dev Test Df Deviance Pr(>Chi)
## 1 1 23 95.12606 NA NA NA
## 2 Garden 21 29.33699 2 65.78907 5.94098e-11La regresión de cuasi-Poisson es útil ya que tiene un parámetro de dispersión variable, de modo que puede modelar datos muy dispersos. Puede ser mejor que la regresión binomial negativa en algunas circunstancias (Verhoef y Boveng. 2007).
En el momento de escribir este artículo, la regresión de Quasi-Poisson no tiene un conjunto completo de funciones de apoyo en R. Al usar la opción de familia de quasipoisson en la función glm, los resultados tendrán los mismos coeficientes de parámetro que con la opción de Poisson, pero las estadísticas de inferencia se ajustan en la función de resumen. La función Anova en el paquete car se puede usar para un análisis de la tabla de desviación, y el paquete emmeans se puede usar para comparaciones post-hoc. Dado que el modelo no produce un valor de probabilidad logarítmica, no conozco una forma de producir un valor p para la moda, para un valor pseudo R cuadrado para el modelo.
model.qp = glm(Monarchs ~ Garden,
data=Data,
family="quasipoisson")
library(car)
Anova(model.qp,
type="II",
test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## LR Chisq Df Pr(>Chisq)
## Garden 52.286 2 4.429e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(multcompView)
library(emmeans)
marginal = emmeans(model.qp,
~ Garden)
pairs(marginal,
adjust="tukey")## contrast estimate SE df z.ratio p.value
## A - B -1.312 0.339 Inf -3.866 0.0003
## A - C -1.904 0.323 Inf -5.896 <.0001
## B - C -0.592 0.195 Inf -3.038 0.0067
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatescld(marginal,
alpha=0.05,
Letters=letters,
adjust="tukey")## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## Garden emmean SE df asymp.LCL asymp.UCL .group
## A 0.56 0.301 Inf -0.16 1.28 a
## B 1.87 0.156 Inf 1.50 2.25 b
## C 2.46 0.116 Inf 2.19 2.74 c
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05La prueba de Vuong, implementada por el paquete pscl, puede probar dos modelos no anidados. Funciona con negbin, zeroinfl y algunos objetos de modelo glm que se ajustan a los mismos datos.
La hipótesis nula es que no hay diferencia en los modelos. La función produce tres pruebas, una prueba “sin procesar”, una corregida por AIC y una corregida por BIC, cualquiera de las cuales podría usarse.
Se ha sugerido que la prueba de Vuong no se utilice para probar la inflación cero (Wilson, 2015).
model.p = glm(Monarchs ~ Garden,
data=Data,
family="poisson")
library(MASS)
model.nb = glm.nb(Monarchs ~ Garden,
data=Data,
control = glm.control(maxit=10000))## Warning in sqrt(1/i): Se han producido NaNs
## Warning in sqrt(1/i): Se han producido NaNslibrary(pscl)
model.zi = zeroinfl(Monarchs ~ Garden,
data = Data,
dist = "poisson")library(pscl)
vuong(model.p,
model.nb,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw 0.03324988 model1 > model2 0.48674
## AIC-corrected 0.03324988 model1 > model2 0.48674
## BIC-corrected 0.03324988 model1 > model2 0.48674vuong(model.p,
model.zi,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.4424725 model2 > model1 0.074585
## AIC-corrected -0.4335210 model2 > model1 0.332318
## BIC-corrected 0.1607786 model1 > model2 0.436134vuong(model.nb,
model.zi,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.4424725 model2 > model1 0.074585
## AIC-corrected -0.4335210 model2 > model1 0.332318
## BIC-corrected 0.1607786 model1 > model2 0.436134