10
$\begingroup$

In R we can "prior weight" a glm regression via the weights parameter. For example:

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

How can this be accomplished in a JAGS or BUGS model?

I found some paper discussing this, but none of them provides an example. I'm interested mainly into Poisson and logistic regression examples.

$\endgroup$
  • $\begingroup$ +1 very good question! I was asking this some bayesian specialists and they only say that in some cases (weights according to categorical covariate), you can compute the posterior distribution of parameter(s) for each category and then combine them in a weighted mean. They did not give me a general solution though, so I would be really interested if it exists or not! $\endgroup$ – Curious Nov 19 '13 at 15:38
7
$\begingroup$

It might be late... but,

Please note 2 things:

  • Adding data points is not advised as it would change degrees of freedom. Mean estimates of fixed effect could be well estimated, but all inference should be avoided with such models. It is hard to "let the data speaks" if you change it.
  • Of course it only works with integer-valued weights (you cannot duplicate 0.5 data point), which is not what is done in most weighted (lm) regression. In general, a weighing is created based on local variability estimated from replicates (e.g. 1/s or 1/s^2 at a given 'x') or based on response height (e.g. 1/Y or 1/Y^2, at a given 'x').

In Jags, Bugs, Stan, proc MCMC, or in Bayesian in general, the likelihood is not different than in frequentist lm or glm (or any model), it is just the same !! Just create a new column "weight" for your response, and write the likelihood as

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

Or a weighted poisson:

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

This Bugs/Jags code would simply to the trick. You will get everything correct. Don't forget to continue multiplying the posterior of tau by the weight, for instance when making prediction and confidence/prediction intervals.

$\endgroup$
  • $\begingroup$ If we do it as stated, we change the mean and the variance. Why is it not y[i] * weight[i] ~ dpois(lambda[i] * weight[i])? That would adjust only the variance. The problem here is that y[i] * weight[i] might be of type real. $\endgroup$ – user28937 Jan 21 '15 at 21:19
  • $\begingroup$ indeed, weighted regression does change mean (because the weighing leads the regression to go closer to the points that have a lot of weights !) and the variance is now a function of the weights (hence it is not an homoskedastic model). The variance (or precision) tau has no meaning any longer, but tau/weight[i] can be interpreted exactly as the precision of the model (for a given "x"). I would not advise the multiplication of the data (y) by the weights... expect if this is precisely something you want to do, but I don't understand you model in this case... $\endgroup$ – Pierre Lebrun Jan 27 '15 at 9:57
  • $\begingroup$ I agree with you it doesn't changes the mean in normal example: y[i] ~ dnorm(mu[i], tau / weight[i]), but it does in the second, since lambda[i] * weight[i] becomes the "new" lambda for dpois and this is not going to match to y[i] anymore. I have to correct myself it should be: t y[i] * exp(weight[i]) ~ dpois(lambda[i] * weight[i]). The idea with the multiplication in the poisson case is that we want to adjust the variance, but also adjust the mean, so don't we need to correct the mean? $\endgroup$ – user28937 Jan 28 '15 at 15:35
  • $\begingroup$ If you need to adjust the variance independently, maybe a Negative Binomial model may come in handy instead of a Poisson ? It adds a variance inflation/deflation parameter to the Poisson. Except that it is very similar. $\endgroup$ – Pierre Lebrun Jan 29 '15 at 11:37
  • $\begingroup$ Pierre good idea. I thought also about the continuous representation of the poisson distribution defined in slide 6/12 at linkd $\endgroup$ – user28937 Jan 30 '15 at 15:30
4
$\begingroup$

First, it's worth pointing out thatglm does not perform bayesian regression. The 'weights' parameter is basically a short hand for "proportion of observations," which can be replaced with up-sampling your dataset appropriately. For example:

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)

So to add weight to points in JAGS or BUGS, you can augment your dataset in a similar fashion as above.

$\endgroup$
  • 2
    $\begingroup$ this won't work, weights are usualy real numbers, not integer $\endgroup$ – Curious Nov 19 '13 at 15:31
  • $\begingroup$ This doesn't preclude you from approximating them with integers. My solution isn't perfect, but it works approximately. For example, given weights (1/3, 2/3, 1), you could upsample the second class by a factor of two and the third class by a factor of three. $\endgroup$ – David Marx Nov 19 '13 at 17:59
0
$\begingroup$

Tried adding to comment above, but my rep is too low.

Should

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

not be

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

in JAGS? I'm running some tests comparing results from this method in JAGS to results from a weighted regression via lm() and can only find accordance using the latter. Here's a simple example:

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)

and compare to

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
$\endgroup$
  • $\begingroup$ Regardless of reputation comments should not be given as answers. $\endgroup$ – Michael R. Chernick Aug 28 '18 at 23:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.