¿Cómo maneja glmnet la sobredispersión?

Tengo una pregunta sobre cómo modelar texto sobre datos de conteo, en particular cómo podría usar la lassotécnica para reducir las características.

Digamos que tengo N artículos en línea y el recuento de páginas vistas para cada artículo. Extraje 1 gramo y 2 gramos para cada artículo y quería ejecutar una regresión sobre los 1,2 gramos. Dado que las características (1,2 gramos) son mucho más que el número de observaciones, el lazo sería un buen método para reducir el número de características. Además, he encontrado que glmnetes muy útil para ejecutar el análisis de lazo.

Sin embargo, el número de recuento de páginas vistas se overdispersed (varianza> media), pero glmnetno ofrece quasipoisson(explícitamente) o negative binomialsino poissonpara datos de recuento. La solución en la que he pensado es log transformcontar los datos de conteo (un método comúnmente utilizado entre los científicos sociales) y hacer que la variable de respuesta siga aproximadamente una distribución normal. Como tal, posiblemente podría modelar los datos con la familia gaussiana usando glmnet.

Entonces mi pregunta es: ¿es apropiado hacerlo? O, ¿debo usar poisson para las manijas de las glmnetcajas ? ¿O hay otros paquetes R que manejan esta situación?glmnetquasipoisson

¡Muchas gracias!

Respuesta corta

¡La sobredispersión no importa al estimar un vector de coeficientes de regresión para la media condicional en un modelo cuasi / poisson! Estará bien si olvida la sobredispersión aquí, use glmnet con la familia de Poisson y solo concéntrese en si su error de predicción con validación cruzada es bajo.

La calificación sigue a continuación.

Poisson, cuasi-Poisson y funciones de estimación:

Digo lo anterior porque la sobredispersión (OD) en un modelo de poisson o cuasi-poisson influye en cualquier cosa que tenga que ver con la dispersión (o varianza o escala o heterogeneidad o propagación o como quiera llamarlo) y como tal tiene un efecto en el estándar errores e intervalos de confianza, pero deja intactas las estimaciones para la media condicional de (llamada μ ). Esto se aplica particularmente a las descomposiciones lineales de la media, como x ⊤$y$$\mu$$x^\top\beta$ .

Esto viene del hecho de que las ecuaciones de estimación para los coeficientes de la media condicional son prácticamente las mismas para los modelos de poisson y cuasi-poisson. Cuasi-poisson especifica la función de varianza en términos de la media y un parámetro adicional (digamos ) como V a r ( y ) = θ μ (con Poisson θ = 1), pero el θ no resulta relevante cuando optimizando la ecuación de estimación. Por lo tanto, θ no juega ningún papel en la estimación de β cuando la media condicional y la varianza son proporcionales. Por lo tanto las estimaciones puntuales beta$\theta$$Var(y)=\theta\mu$$\theta$$\theta$$\theta$$\beta$$\hat{\beta}$ son idénticos para los modelos cuasi y poisson!

Permítanme ilustrar con un ejemplo (observe que uno necesita desplazarse para ver todo el código y la salida):

> library(MASS)
> data(quine) 
> modp <- glm(Days~Age+Sex+Eth+Lrn, data=quine, family="poisson")
> modqp <- glm(Days~Age+Sex+Eth+Lrn, data=quine, family="quasipoisson")
> summary(modp)

Call:
glm(formula = Days ~ Age + Sex + Eth + Lrn, family = "poisson", 
    data = quine)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-6.808  -3.065  -1.119   1.819   9.909  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.71538    0.06468  41.980  < 2e-16 ***
AgeF1       -0.33390    0.07009  -4.764 1.90e-06 ***
AgeF2        0.25783    0.06242   4.131 3.62e-05 ***
AgeF3        0.42769    0.06769   6.319 2.64e-10 ***
SexM         0.16160    0.04253   3.799 0.000145 ***
EthN        -0.53360    0.04188 -12.740  < 2e-16 ***
LrnSL        0.34894    0.05204   6.705 2.02e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2073.5  on 145  degrees of freedom
Residual deviance: 1696.7  on 139  degrees of freedom
AIC: 2299.2

Number of Fisher Scoring iterations: 5

> summary(modqp)

Call:
glm(formula = Days ~ Age + Sex + Eth + Lrn, family = "quasipoisson", 
    data = quine)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-6.808  -3.065  -1.119   1.819   9.909  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.7154     0.2347  11.569  < 2e-16 ***
AgeF1        -0.3339     0.2543  -1.313 0.191413    
AgeF2         0.2578     0.2265   1.138 0.256938    
AgeF3         0.4277     0.2456   1.741 0.083831 .  
SexM          0.1616     0.1543   1.047 0.296914    
EthN         -0.5336     0.1520  -3.511 0.000602 ***
LrnSL         0.3489     0.1888   1.848 0.066760 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 13.16691)

    Null deviance: 2073.5  on 145  degrees of freedom
Residual deviance: 1696.7  on 139  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Como puede ver a pesar de que tenemos una fuerte sobredispersión de 12.21 en este conjunto de datos (por deviance(modp)/modp$df.residual) los coeficientes de regresión (estimaciones puntuales) no cambian en absoluto. Pero observe cómo cambian los errores estándar.

La cuestión del efecto de la sobredispersión en modelos de Poisson penalizados

$\theta$$\beta$

$g(\mu)=x^\top\beta + f(\beta)$

$\beta$$\theta \mu$

glmnet$f(\beta)=0$

> library(glmnet)
> y <- quine[,5]
> x <- model.matrix(~Age+Sex+Eth+Lrn,quine)
> modl <- glmnet(y=y,x=x, lambda=c(0.05,0.02,0.01,0.005), family="poisson")
> coefficients(modl)
8 x 4 sparse Matrix of class "dgCMatrix"
                    s0         s1         s2         s3
(Intercept)  2.7320435  2.7221245  2.7188884  2.7172098
(Intercept)  .          .          .          .        
AgeF1       -0.3325689 -0.3335226 -0.3339580 -0.3340520
AgeF2        0.2496120  0.2544253  0.2559408  0.2567880
AgeF3        0.4079635  0.4197509  0.4236024  0.4255759
SexM         0.1530040  0.1581563  0.1598595  0.1607162
EthN        -0.5275619 -0.5311830 -0.5323936 -0.5329969
LrnSL        0.3336885  0.3428815  0.3459650  0.3474745

Entonces, ¿qué hace OD a los modelos de regresión penalizados? Como ya sabrán, todavía existe cierto debate sobre la forma correcta de calcular los errores estándar para los modelos penalizados (ver, por ejemplo, aquí ) y de glmnettodos modos no se está produciendo, probablemente por esa razón. Es muy posible que el OD influya en la parte de inferencia del modelo, tal como lo hace en el caso no penalizado, pero a menos que se llegue a un consenso sobre la inferencia en este caso, no lo sabremos.

Por otro lado, uno puede dejar todo este desorden si está dispuesto a adoptar una visión bayesiana donde los modelos penalizados son solo modelos estándar con un previo específico.