Regressão generalizada ponderada em BUGS, JAGS

10

No Rque pudermos "peso antes" uma glmregressão via pesos parâmetro. Por exemplo:

glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(), weights=w)

Como isso pode ser feito em um JAGSou BUGSmodelo?

Encontrei algum artigo discutindo isso, mas nenhum deles fornece um exemplo. Estou interessado principalmente em exemplos de Poisson e de regressão logística.

user28937
fonte
+1 muito boa pergunta! Perguntei a alguns especialistas bayesianos e eles apenas dizem que em alguns casos (pesos de acordo com a covariável categórica), é possível calcular a distribuição posterior dos parâmetros para cada categoria e combiná-los em uma média ponderada. Eles não me deram uma solução geral, então eu estaria realmente interessado se ele existe ou não!
Curioso

Respostas:

7

Pode ser tarde ... mas,

Observe duas coisas:

  • Não é recomendável adicionar pontos de dados, pois alteraria os graus de liberdade. As estimativas médias de efeito fixo podem ser bem estimadas, mas toda inferência deve ser evitada com esses modelos. É difícil "deixar os dados falarem" se você os alterar.
  • Obviamente, ele só funciona com pesos com valor inteiro (você não pode duplicar 0,5 ponto de dados), o que não é o que é feito na regressão mais ponderada (lm). Em geral, uma pesagem é criada com base na variabilidade local estimada a partir de réplicas (por exemplo, 1 / s ou 1 / s ^ 2 em um dado 'x') ou com base na altura da resposta (por exemplo, 1 / Y ou 1 / Y ^ 2, em um determinado 'x').

Em Jags, Bugs, Stan, proc MCMC, ou em Bayesiano em geral, a probabilidade não é diferente do que em frequentist lm ou glm (ou qualquer modelo), é a mesma coisa !! Basta criar uma nova coluna "peso" para sua resposta e escrever a probabilidade como

y [i] ~ dnorm (mu [i], tau / peso [i])

Ou um poisson ponderado:

y [i] ~ dpois (lambda [i] * peso [i])

Este código de Bugs / Jags seria simplesmente o truque. Você vai conseguir tudo correto. Não se esqueça de continuar multiplicando a parte posterior da tau pelo peso, por exemplo, ao fazer intervalos de previsão e confiança / previsão.

Pierre Lebrun
fonte
Se fizermos como indicado, alteramos a média e a variância. Por que não é y [i] * peso [i] ~ dpois (lambda [i] * peso [i])? Isso ajustaria apenas a variação. O problema aqui é que y [i] * peso [i] pode ser do tipo real.
precisa saber é o seguinte
de fato, a regressão ponderada muda de média (porque a pesagem leva a regressão a se aproximar dos pontos que têm muitos pesos!) e a variação agora é uma função dos pesos (portanto, não é um modelo homosquástico). A variação (ou precisão) tau não tem mais significado, mas tau / peso [i] pode ser interpretado exatamente como a precisão do modelo (para um dado "x"). Eu não recomendaria a multiplicação dos dados (y) pelos pesos ... espere se isso é precisamente algo que você deseja fazer, mas não entendo o seu modelo neste caso ...
Pierre Lebrun
Eu concordo com você que não altera a média no exemplo normal: y [i] ~ dnorm (mu [i], tau / peso [i]), mas muda no segundo, desde lambda [i] * peso [ i] se torna o "novo" lambda para dpois e isso não corresponde mais a y [i]. Eu tenho que me corrigir, deve ser: ty [i] * exp (peso [i]) ~ dpois (lambda [i] * peso [i]). A idéia com a multiplicação no caso de Poisson é que queremos ajustar a variação, mas também ajustar a média, por isso não precisamos corrigir a média?
precisa saber é o seguinte
Se você precisar ajustar a variação de forma independente, talvez um modelo Binomial Negativo possa ser útil em vez de um Poisson? Ele adiciona um parâmetro de inflação / deflação de variação ao Poisson. Exceto que é muito parecido.
Pierre Lebrun 29/01
Pierre boa ideia. Eu também pensei sobre a representação contínua da distribuição de Poisson definido na corrediça 6/12 em linkd
user28937
4

Primeiro, vale ressaltar que glmnão realiza regressão bayesiana. O parâmetro 'pesos' é basicamente uma mão curta para "proporção de observações", que pode ser substituída pela amostragem adequada de seu conjunto de dados. Por exemplo:

x=1:10
y=jitter(10*x)
w=sample(x,10)

augmented.x=NULL
augmented.y=NULL    
for(i in 1:length(x)){
    augmented.x=c(augmented.x, rep(x[i],w[i]))
    augmented.y=c(augmented.y, rep(y[i],w[i]))
}

# These are both basically the same thing
m.1=lm(y~x, weights=w)
m.2=lm(augmented.y~augmented.x)

Portanto, para adicionar peso aos pontos no JAGS ou BUGS, você pode aumentar seu conjunto de dados de maneira semelhante à anterior.

David Marx
fonte
2
isso não vai funcionar, pesos são usualy números reais, não inteiros
Curious
Isso não impede que você os aproxime com números inteiros. Minha solução não é perfeita, mas funciona aproximadamente. Por exemplo, considerando os pesos (1/3, 2/3, 1), você pode aumentar a amostra da segunda classe por um fator de dois e a terceira classe por um fator de três.
David Marx
0

Tentei adicionar ao comentário acima, mas meu representante está muito baixo.

Devemos

y[i] ~ dnorm(mu[i], tau / weight[i])

não ser

y[i] ~ dnorm(mu[i], tau * weight[i])

em JAGS? Estou executando alguns testes comparando resultados desse método no JAGS com resultados de uma regressão ponderada via lm () e só posso encontrar concordância usando o último. Aqui está um exemplo simples:

aggregated <- 
  data.frame(x=1:5) %>%
  mutate( y = round(2 * x + 2 + rnorm(length(x)) ),
          freq = as.numeric(table(sample(1:5, 100, 
                 replace=TRUE, prob=c(.3, .4, .5, .4, .3)))))
x <- aggregated$x
y <- aggregated$y
weight <- aggregated$freq
N <- length(y)

# via lm()
lm(y ~ x, data = aggregated, weight = freq)

e compare com

lin_wt_mod <- function() {

  for (i in 1:N) {
    y[i] ~ dnorm(mu[i], tau*weight[i])
    mu[i] <- beta[1] + beta[2] * x[i]
  }

  for(j in 1:2){
    beta[j] ~ dnorm(0,0.0001)
  }

  tau   ~ dgamma(0.001, 0.001)
  sigma     <- 1/sqrt(tau)
}

dat <- list("N","x","y","weight")
params <- c("beta","tau","sigma")

library(R2jags)
fit_wt_lm1 <- jags.parallel(data = dat, parameters.to.save = params,
              model.file = lin_wt_mod, n.iter = 3000, n.burnin = 1000)
fit_wt_lm1$BUGSoutput$summary
metrikrn
fonte
Independentemente da reputação, os comentários não devem ser dados como respostas.
Michael R. Chernick 28/08