Intervalos de previsão pragmáticos de um GLM de quase-probabilidade por @ ellis2013nz

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


O blog de hoje vem com duas lições: uma estatística e outra sobre solução de problemas.

Depuração de pato de borracha

Lição de solução de problemas primeiro. O problema descrito abaixo me causa pesar há vários dias. Tenho um caso de uso importante e de tempo crítico em funcionamento, no qual preciso atribuir muitos dados ausentes complexos usando o tipo de modelo que estou prestes a descrever. Mas, devido ao cansaço que me levou a não pensar direito, e ao problema de ser incorporado a um conjunto muito maior de programas de computador complexos, eu estava lutando para obter a resposta. Depois de um tempo, percebi corretamente que precisava abstrair o problema do caso de uso e fazer um exemplo com dados simulados para esclarecer minhas idéias.

Depois de naquela não funcionou apesar de algumas horas de esforço embaraçosas, no fim de semana em que decidi assumir a incrível responsabilidade de escrever uma nova pergunta sobre Validação cruzada. Eu estava relutante em fazer isso porque eu sou claramente pretendia saber a resposta para esse problema (descrito abaixo), portanto parece uma admissão de falha profissional ter que perguntar publicamente. Só porque eu estava genuinamente confuso, eu iria tão longe; parecia claro que o que eu pensava entender sobre modelos lineares generalizados implementados em R estava errado em algum aspecto importante.

Passei cerca de 30 minutos preparando cuidadosamente um exemplo mínimo reproduzível, explicando claramente o que eu precisava fazer e o que tentara, qual era meu pensamento e por que isso não parecia funcionar, relendo meticulosamente toda a documentação. O que me fez pensar: “oh, talvez isso signifique isso …”. Sem muita esperança (pois era a terceira vez no processo, durante vários dias, de repente, pensei ter visto a luz), fiz um exemplo para mostrar que havia tentado isso e por que não funcionou. Mas (quem teria previsto isso) fez trabalhos.

Portanto, nunca subestime o poder da depuração de pato de borracha. Ou seja, a superpotência de solução de problemas que está escrevendo um exemplo reproduzível e explicando cuidadosamente exatamente qual é o seu problema de uma maneira que você acha que outra pessoa pode responder por você. Se você explicar com bastante clareza, é provável que outra pessoa seja você.

Leia Também  Crescimento de cartões na Europa diminui para 2,7%, com pagamentos digitais começando a crescer

Na verdade, possuo quatro patos de borracha expressamente com o objetivo de me ajudar a resolver problemas difíceis que “certamente” devem estar em meus próprios recursos para resolver – um para cada um dos meus postos de trabalho em casa, um para o trabalho (agora como reposição em casa ) e um quarto apenas para viagem. Mas eu não segui meu próprio conselho de explicar o problema a um dos patos. Então, lição aprendida.

Intervalos de previsão para variáveis ​​com relações exponenciais e variação crescente

Para o problema estatístico substantivo. Eu tenho dados que, quando transformados em log, têm um relacionamento linear com outros dados. E minha variável principal claramente tem uma variação crescente à medida que seu valor esperado aumenta. No meu exemplo do mundo real, tenho muitas dessas variáveis ​​e a motivação é imputar valores ausentes.

Quero modelar os dados de uma maneira que não descarte valores negativos individuais, portanto, realizar uma transformação de log é sub-ideal. Em vez disso, quero usar um modelo linear generalizado com a função de link do logaritmo, o que significa que o valor esperado da resposta nunca ficará abaixo de zero, mas os valores individuais podem (os dados originais incluem variáveis ​​financeiras que seguem esse padrão). Eu preciso da variação da variável de resposta para aumentar com o valor esperado, então não posso usar family = gaussian(link = log), que assume variação constante.

A variação da resposta no nível individual é crucial, porque vou usar o modelo para simular valores individuais em uma implementação caseira de imputação múltipla com equações encadeadas (por várias razões complicadas, não posso simplesmente usar o mice pacote, que teria sido a minha preferência e é mais fácil rodar manualmente meu próprio sistema). Estou preparado para fingir que os dados são normalmente distribuídos (condicional a um valor esperado que muda com as outras variáveis ​​do sistema), desde que eu tenha uma variação crescente proporcional à média. Portanto, um resultado final importante de tudo isso é um vetor que eu posso usar para o sd = argumento de rnorm().

As duas estratégias de modelagem em discussão são

  • transformar log e transformar quadrados mínimos ordinários; ou
  • um GLM de quase probabilidade com uma função de link de log e variação proporcional à média da resposta. (Como discutido acima, esse é o objetivo que eu quero usar. No final da pista, eu poderia usar “ao quadrado da média da resposta”, mas isso não é importante para fins ilustrativos.)
Leia Também  Lançada a nova empresa consultiva Avonhurst

Essas estratégias fornecem estimativas de coeficiente semelhantes. Vejamos alguns dados simulados (retidos da minha pergunta interrompida para validação cruzada)

library(tidyverse)
library(scales)

n  50
set.seed(123)

d  tibble(x = rnorm(n, 4, 1)) %>%
  mutate(mean = 2 + 0.25 * x,
		 sd = sqrt(mean) / 20,
		 y = exp(mean + rnorm(n, , sd))
  )

mod1  lm(log(y) ~ x, data = d)
mod2  glm(y ~ x, family = quasi(link = log, variance = mu), data = d)

# similar coefficients (one is minimising squares on the log scale, 
# the other is minimising weighted squares on the original)
rbind(coef(mod1), coef(mod2))

O que nos dá:

         (Intercept)         x
    [1,]    2.021225 0.2478562
    [2,]    2.016842 0.2496386

Assim, como um exercício de ajuste de curva, ambos os métodos são muito bons para recuperar os valores reais originais de 2 e 0,25 para a interceptação e a inclinação, respectivamente, no preditor linear do meu processo de geração de dados.

Um intervalo de previsão difere de um intervalo de confiança, pois inclui aleatoriedade em nível individual além de a incerteza nas estimativas de parâmetros do modelo. Por exemplo, geom_smooth() desenha um intervalo de confiança, expressando confiança / incerteza na posição e na forma da linha de melhor ajuste, mas não no intervalo em que se espera que 95% dos novos valores apareçam.

Intervalo de previsão com transformação e mínimos quadrados comuns

Criar um intervalo de previsão a partir do modelo OLS transformado por log é um exercício padrão de Regressão 101. A variação da estimativa para qualquer indivíduo é a soma da variação de nossa curva estimada nesse ponto e a variação residual dos indivíduos. Pegue a raiz quadrada e você obtém o desvio padrão; assuma a normalidade e multiplique por 1,96 para obter um intervalo de previsão de 95%; voltar para a escala original (observe que há complicações de ajuste de viés com as quais poderíamos nos preocupar aqui, mas ignoraremos para os propósitos de hoje).

… Gerado com este código:

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
#------------Prediction interval from the log transform version--------------------
sx  tibble(x = seq(from = 1, to = 7, length.out = 100)) 
pred1  predict(mod1, se.fit = TRUE, newdata = sx)

sx1  sx %>%
  # calculate standard deviation for the prediction interval and create lower and upper bounds:
  mutate(se_pi = sqrt(pred1$se.fit ^ 2 + pred1$residual.scale ^ 2),
         lower_pi = exp(pred1$fit  - 1.96 * se_pi),
         upper_pi =exp(pred1$fit +  1.96 * se_pi))

p1  ggplot(d, aes(x = x)) +
  geom_ribbon(data = sx1, aes(ymin= lower_pi, ymax = upper_pi), fill = "steelblue", alpha = 0.1) +
  geom_point(aes(y = y))+
  labs(title = "Original data with prediction interval from log transform model")

Para fazer algo semelhante ao modelo linear generalizado com a resposta de quase-probabilidade, precisamos entender o residual.scale argumento retornado por predict.glm(). Se usamos predict(..., type = "response") ou predict(..., type = "link"), a residual.scale o valor é o mesmo (0,361 no nosso caso). Mas 0,361 de quê?

Leia Também  X é para scale_x | R-bloggers

residual.scale é descrito no arquivo de ajuda como “um escalar que fornece a raiz quadrada da dispersão usada no cálculo dos erros padrão”. Como usamos isso para fazer um intervalo de previsão? Em nossa especificação de mod2, Eu disse glm que a variação de y seria proporcional à sua média (family = quasi(link = log, variance = mu)) Podemos recuperar a variação estimada em um ponto individual pelo valor previsto de y multiplicado pelo fator de dispersão. Portanto, isso significa que, nesse caso específico, podemos gerar o intervalo de previsão da seguinte maneira. A parte principal da nossa discussão é + pred2$fit * pred2$residual.scale ^ 2) – Porque pred2$fit é a média de y no ponto especificado e pred2$residual_scale ^ 2 é o fator de dispersão que escala a média para se tornar sua variação.

#-----------------prediction interval from the quasi family GLM------------

pred2  predict(mod2, se.fit = TRUE, newdata = sx, type = "response")
sx2  sx %>%
  mutate(se_pi = sqrt(pred2$se.fit ^ 2 + pred2$fit * pred2$residual.scale ^ 2),
         lower_pi = pred2$fit - 1.96 * se_pi,
         upper_pi = pred2$fit + 1.96 * se_pi)

ggplot(d, aes(x = x)) +
  geom_ribbon(data = sx2, aes(ymin= lower_pi, ymax = upper_pi), fill = "steelblue", alpha = 0.1) +
  geom_point(aes(y = y))+
  labs(title = "Original data with prediction interval from quasi likelihood glm with log link")

E aqui vamos nós:

Parece simples quando você sabe.

Existem diferenças sutis, mas importantes, entre os intervalos de previsão de nossos dois métodos. Como escrevi o processo original de geração de dados, sei que o segundo está mais próximo do modelo “real”. O modelo de transformação de log, na verdade, aumenta a variação muito rapidamente à medida que a média aumenta, em comparação com o processo real. Para valores mais altos de y, o intervalo de previsão é maior do que o necessário. Obviamente, se esse será o caso com dados reais, dependerá da situação específica.

Observe que esse método não funcionará para todos os GLMs. Em particular, qualquer coisa com uma resposta não contínua precisaria definitivamente de outro pensamento. E para usar aqueles +/- 1.96 valores para o intervalo de previsão, confiei em estar preparado para fazer uma suposição simplificadora sobre a distribuição normal da variável de resposta para os fins do intervalo de previsão, que eu não fiz na estimativa original do modelo de quase-probabilidade. Não existe um teorema do limite central ou qualquer coisa em que confiar aqui, porque estes são Individual valores de y, e eles terão qualquer forma que y normalmente tenha.

O que mais importa aqui é entender como isso residual.scale estimativa de dispersão funciona, em relação à especificação de variação na chamada original para quasi(). Essa relação com a variação não depende de qualquer suposição de normalidade adicional.

De qualquer forma, use essa abordagem com cuidado. Obviamente, se eu realmente Quando estávamos preparados para dizer que y é normalmente distribuído (mas com a variação aumentando proporcional à média), eu poderia modelá-lo explicitamente como tal, em vez de brincar com quase-probabilidades. Então, o que é descrito aqui é muito sobre pragmatismo.



Se você chegou até aqui, por que não inscreva-se para atualizações do site? Escolha o seu sabor: o email, Twitter, RSS ou facebook …



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br