Взвешенная обобщенная регрессия в BUGS, JAGS

10

В R«предварительном весе» мы можем glmрегрессировать через параметр весов . Например:

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

Как это можно сделать в JAGSили BUGSмодели?

Я нашел некоторые документы, обсуждающие это, но ни один из них не дает пример. Меня интересуют в основном примеры Пуассона и логистической регрессии.

user28937
источник
+1 очень хороший вопрос! Я спрашивал об этом некоторых байесовских специалистов, и они только говорят, что в некоторых случаях (веса согласно категориальному ковариате) вы можете вычислить апостериорное распределение параметра (ов) для каждой категории и затем объединить их в средневзвешенном значении. Они не дали мне общего решения, поэтому мне было бы очень интересно, существует он или нет!
Любопытно

Ответы:

7

Это может быть поздно ... но,

Обратите внимание на 2 вещи:

  • Добавление точек данных не рекомендуется, так как это изменит степени свободы. Средние оценки фиксированного эффекта могут быть хорошо оценены, но с такими моделями следует избегать любых выводов. Трудно «позволить данным говорить», если вы их измените.
  • Конечно, он работает только с целочисленными весами (вы не можете дублировать точку данных 0,5), что не является тем, что делается в большинстве взвешенных (lm) регрессий. В общем, взвешивание создается на основе локальной изменчивости, оцененной по повторностям (например, 1 / с или 1 / с ^ 2 при заданном «x») или на основе высоты отклика (например, 1 / Y или 1 / Y ^ 2, при данный «х»).

В Jags, Bugs, Stan, proc MCMC или в байесовском в целом вероятность не отличается от частоты в lm или glm (или любой модели), она просто одинакова !! Просто создайте новый столбец «вес» для вашего ответа и напишите вероятность как

y [i] ~ dnorm (му [i], тау / вес [i])

Или взвешенный пуассон:

y [i] ~ dpois (лямбда [i] * вес [i])

Этот код Bugs / Jags был бы просто приколом. Вы получите все правильно. Не забывайте продолжать умножать апостериор тау на вес, например, при прогнозировании и доверительных интервалах.

Пьер Лебрун
источник
Если мы сделаем это, как указано, мы изменим среднее значение и дисперсию. Почему это не у [я] * вес [я] ~ dpois (лямбда [я] * вес [я])? Это будет регулировать только дисперсию. Проблема здесь в том, что y [i] * weight [i] может иметь тип real.
user28937
действительно, взвешенная регрессия действительно меняет среднее (потому что взвешивание приводит к тому, что регрессия приближается к точкам, имеющим много весов!), и дисперсия теперь является функцией весов (следовательно, это не модель гомоскедастика). Дисперсия (или точность) тау больше не имеет значения, но тау / вес [i] можно интерпретировать точно как точность модели (для данного «x»). Я бы не советовал умножать данные (y) на весовые коэффициенты ... ожидайте, если это именно то, что вы хотите сделать, но я не понимаю вашу модель в этом случае ...
Pierre Lebrun
Я согласен с вами, это не меняет среднего значения в обычном примере: y [i] ~ dnorm (mu [i], tau / weight [i]), но это происходит во втором, так как lambda [i] * weight [ i] становится «новой» лямбда-версией для dpois, и это больше не будет соответствовать y [i]. Я должен исправить себя, это должно быть: ty [i] * exp (weight [i]) ~ dpois (lambda [i] * weight [i]). Идея с умножением в случае Пуассона состоит в том, что мы хотим скорректировать дисперсию, но также скорректировать среднее значение, так что нам не нужно корректировать среднее?
user28937
Если вам нужно отрегулировать дисперсию независимо, может быть, вам пригодится модель с отрицательным биномиальным значением вместо пуассоновской? Это добавляет параметр инфляции / дефляции дисперсии к Пуассону. За исключением того, что это очень похоже.
Пьер Лебрун
Пьер хорошая идея. Я также подумал о непрерывном представлении распределения Пуассона, определенного на слайде 6/12 по ссылке
user28937
4

Во-первых, стоит отметить, что glmне выполняет байесовскую регрессию. Параметр «веса» в основном является сокращением от «пропорции наблюдений», который может быть заменен соответствующим повышением частоты вашего набора данных. Например:

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)

Таким образом, чтобы добавить вес к точкам в JAGS или BUGS, вы можете увеличить свой набор данных аналогично тому, как описано выше.

Дэвид Маркс
источник
2
это не сработает, весы - это обычно действительные числа, а не целые числа
Любопытно,
Это не мешает вам аппроксимировать их целыми числами. Мое решение не идеально, но оно работает примерно. Например, с учетом весов (1/3, 2/3, 1) вы можете повысить выборку второго класса в два раза, а третьего класса - в три раза.
Дэвид Маркс
0

Пытался добавить комментарий выше, но моя репутация слишком низкая.

Должен

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

не быть

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

в зазубринах? Я запускаю несколько тестов, сравнивающих результаты этого метода в JAGS с результатами взвешенной регрессии с помощью lm (), и могу найти соответствие только с помощью последнего. Вот простой пример:

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)

и сравнить с

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
источник
Независимо от репутации комментарии не следует давать в качестве ответов.
Майкл Р. Черник