Como o glmnet lida com superdispersão?

9

Eu tenho uma pergunta sobre como modelar texto sobre dados de contagem, em particular como eu poderia usar a lassotécnica para reduzir recursos.

Digamos que eu tenha N artigos on-line e a contagem de visualizações de página para cada artigo. Eu extraí 1 grama e 2 gramas para cada artigo e eu queria fazer uma regressão sobre os 1,2 gramas. Como os recursos (1,2 gramas) são muito mais que o número de observações, o laço seria um bom método para reduzir o número de recursos. Além disso, eu achei glmnetrealmente útil para executar a análise de laço.

No entanto, o número de contagens de exibições de página é superdisperso (variação> média), mas glmnetnão oferece quasipoisson(explicitamente) ou negative binomialapenas poissonpara dados de contagem. A solução que eu pensei é log transformnos dados de contagem (um método comumente usado entre cientistas sociais) e fazer com que a variável de resposta siga aproximadamente uma distribuição normal. Como tal, eu poderia modelar os dados com a família gaussiana usando glmnet.

Então, minha pergunta é: é apropriado fazê-lo? Ou devo usar apenas poisson para identificadores de glmnetcaso ? Ou existem outros pacotes R que lidam com essa situação?glmnetquasipoisson

Muito obrigado!

Sonya S.
fonte

Respostas:

14

Resposta curta

A superdispersão não importa ao estimar um vetor de coeficientes de regressão para a média condicional em um modelo de quase / poisson! Você ficará bem se esquecer a sobredispersão aqui, use o glmnet com a família poisson e se concentre apenas em saber se o erro de previsão com validação cruzada é baixo.

A qualificação segue abaixo.


Poisson, Quasi-Poisson e funções de estimativa:

Digo o acima exposto porque a super-dispersão (DO) em um modelo de poisson ou quase-poisson influencia qualquer coisa relacionada à dispersão (ou variação ou escala ou heterogeneidade ou propagação ou o que você quiser chamar) e, como tal, afeta o padrão erros e intervalos de confiança, mas deixa as estimativas para a média condicional de (chamada ) intocadas. Isso se aplica particularmente a decomposições lineares da média, comoyμxβ .

Isso deriva do fato de que as equações de estimativa para os coeficientes da média condicional são praticamente as mesmas para os modelos de poisson e quase-poisson. Quase-poisson especifica a função de variação em termos da média e um parâmetro adicional (digamos ) como (com para Poisson = 1), mas o não em ser relevante ao otimizar a equação de estimativa. Assim, o não desempenha nenhum papel na estimativa do quando a média e a variância condicionais são proporcionais. Portanto, as estimativas pontuais são idênticas para os modelos quase e poisson!θVar(y)=θμθθθββ^

Deixe-me ilustrar com um exemplo (observe que é preciso rolar para ver o código e a saída completos):

> 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 você pode ver, apesar de termos uma super-dispersão forte de 12,21 neste conjunto de dados (por deviance(modp)/modp$df.residual), os coeficientes de regressão (estimativas pontuais) não mudam. Mas observe como os erros padrão mudam.

A questão do efeito da super-dispersão em modelos de poisson penalizados

Modelos penalizados são usados ​​principalmente para previsão e seleção de variáveis ​​e não (ainda) para inferência. Portanto, as pessoas que usam esses modelos estão interessadas nos parâmetros de regressão para a média condicional, apenas encolheram para zero. Se a penalização for a mesma, as equações de estimativa para as médias condicionais derivadas da probabilidade penalizada (quase) também não dependem de e, portanto, a super - dispersão não importa para as estimativas de em um modelo do tipo:θβ

g(μ)=xβ+f(β)

como é estimado da mesma maneira para qualquer função de variação da forma , novamente para todos os modelos em que a média condicional e a variação são proporcionais. βθμÉ como nos modelos de poisson / quase-pontocom não-penalizados.

Se você não quiser aceitar isso pelo valor de face e evitar a matemática, poderá encontrar suporte empírico no fato de que glmnet, se definir o parâmetro de regularização para 0 (e, portanto, ), você acaba praticamente onde os modelos de poisson e quasipoisson pousam (veja a última coluna abaixo, onde lambda é 0,005).f(β)=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

Então, o que OD faz para modelos de regressão penalizados? Como você deve saber, ainda há algum debate sobre a maneira correta de calcular erros padrão para modelos penalizados (veja, por exemplo, aqui ) e glmnetnão produz nenhum resultado, provavelmente por esse motivo. Pode muito bem ser que o OD influencie a parte de inferência do modelo, exatamente como no caso não penalizado, mas, a menos que seja alcançado algum consenso sobre a inferência nesse caso, não saberemos.

Como um aparte, pode-se deixar toda essa bagunça para trás se estiver disposto a adotar uma visão bayesiana em que modelos penalizados são apenas modelos padrão com um prior específico.

Momo
fonte
@ Mono, obrigado por sua explicação muito detalhada! Aqui está meu entendimento e, por favor, corrija-me se estiver errado: poissone as quasipoissonregressões estimam os coeficientes da mesma maneira e o que eles diferem é como estimam erros padrão e, portanto, a significância. No entanto, para o método do laço, como calcular erros padrão ainda não chegou a um consenso, e, portanto, seu uso atual reside principalmente na seleção de variáveis ​​em vez de inferência. Como tal, não importa se usamos glmnetcom poisson ou quasipoisson, mas o que faz é que esse erro de validação cruzada deve ser minimizado.
Sonya S.
@ Mono, outra observação, eu summary(modqp)mesmo corri e vi que ele tinha exatamente as mesmas estimativas de coeficiente. Acredito que sua resposta beneficiará mais pessoas sobre esse assunto porque eu não encontrei nenhuma, por isso sugiro que você adicione a saída do resumo (modqp) para um exemplo ilustrado ainda melhor. Mais uma vez, muito obrigado!
Sonya S.
11
@ Sonya Yours é um bom resumo. A chave é que, ao estimar os parâmetros para a média condicional, as funções de estimativa (digamos, a função de pontuação) para poisson e quasipoisson são as mesmas! Portanto, não importa para esses parâmetros se há uma penalização ou não, desde que seja a mesma penalização. Eu deixo isso mais claro acima. Obrigado também pelo ponteiro sobre o resumo (modq), mas que já está lá, ele apenas é "encaixotado" em uma tela normal, então é preciso rolar para baixo.
Momo
Ainda me pergunto se é possível que menos variáveis ​​sejam reduzidas em Poisson do que se houvesse uma especificação quase-Poisson, que é mais correta, e provavelmente levaria a uma melhor precisão preditiva do que o modelo de Poisson, porque seu modelo de amostragem é mais correto.
Equilíbrio Brash
Na mesma nota, também pode ser que mais variáveis ​​sejam reduzidas em Poisson do que deveriam ser reduzidas em casos de dispersão UNDER (como quando você está usando um modelo robusto de Poisson para estimar relações de risco relativas para dados 0/1).
Equilíbrio Brash