How to specify Bayesian mixed effects model in BUGS

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

I posted this earlier in the week then retracted the question when I found a good source, not wanting to waste people's time. I haven't made much progress I'm afraid. In trying to be a good citizen here I will make the problem as clear as possible. I suspect there will be few takers.

I have a dataframe in R I wish to analyse in BUGS or R. It is in long format. It consists of multiple observations on 120 individuals, with a total of 885 rows. I am examining the occurrence of a categorical outcome - but that's not really relevant here. The question is about something deeper.

The model I have been using up to here is

  mymodel<-gee(Category ~ Predictor 1 + Predictor 2..family=binomial(link="logit"),  data=mydata,    id=Person)

with a marginal model essentially accounting for the clustering of patients. I then examined

mymodel<-gee(Category ~ Predictor 1 + Predictor 2.. , family=binomial(link="logit"),   corstr = "AR-M",   data=mydata, id=Person)

in order to account for the time ordering of the observations on the individual people.

This didn't change much.

Then I tried to model them using the following set of MCMCPack commands:

mymodel<-MCMCglmm(category~  Predictor1 + Predictor2..,   data=mydata, family=binomial(link="logit"))

An examination of the output was thrilling, showing statistical significance for many predictors. I hailed myself as a newly converted bayesian, until I realised I hadn't accounted for repeated measures within patients.

I understand that I have to account for that. I understand that this may mean fitting a hyperprior for each individual - is that right? What form will this take in BUGS?

Here's a basic log reg model: (kudos to Kruschke, J., Indiana)

    model {  
for( i in 1 : nData ) {    y ~ dbern( mu )    mu <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))  }  
b0 ~ dnorm( 0 , 1.0E-12 )   
for ( j in 1 : nPredictors ) {   
b[j] ~ dnorm( 0 , 1.0E-12 )  }
}

However, no hyperprior here for the individual. Here's my best attempt so far at a within-individual design, accounting for repeated measures within people:

Here's Jackman's model for JAGS

1 model{2
## loop over data for likelihood3 for(i in 1:n){4  y ~ dbern( mu )    mu <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))6 }7 sigma ˜ dunif(0,20)
## prior on standard deviation8 tau <- pow(sigma,-2)
## convert to precision910
## hierarchical model for each state’s intercept & slope11 for(p in 1:50){12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,])
## bivariate normal13 }1415
## means, hyper-parameters16 for(q in 1:2){17 mu[q] ˜ dnorm(0,.0016)

}

Here's my bastard-child model for BUGS

1 model{2
## loop over data for likelihood3 for(i in 1:n){4 mu.y <- alpha + beta[s,1] + beta[s,2]*(j-jbar)5 demVote ˜ dnorm(mu.y,tau)6 }7 sigma ˜ dunif(0,20)
## prior on standard deviation8 tau <- pow(sigma,-2)
## convert to precision910
## hierarchical model for each state’s intercept & slope11 for(p in 1:120){12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,])
## bivariate normal13 }1415
## means, hyper-parameters16 for(q in 1:2){17 mu[q] ˜ dnorm(0,.0016)  }

Can somebody let me know if I'm heading in the right direction. My understanding of this is growing, but slowly. Please be gentle. I'm a medic, not a statistic! I have used R quite a bit, but I'm new to BUGS, and new to Bayes.

Thanks

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Bayesian effects specify Effect model multiple possible question earlier problem

You are (were) almost there. Just a few comments - you don't have to make the prior for the beta[,1:2] parameters a joint MV normal; you can make the prior such that beta[i,1] and beta[i,2] are independent, which simplifies things (for example, no prior covariance need be specified.) Note that doing so doesn't mean they will be independent in the posterior.

Other comments: Since you have a constant term - alpha - in the regression, the components beta[,1] should have zero mean in the prior. Also, you don't have a prior for alpha in the code.

Here's a model with hierarchical intercept and slope terms; I've tried to keep to your priors and notation where possible, given the changes:

model {
  for(i in 1:n){
    mu.y[i] <- alpha + beta0[s[i]] + beta1[s[i]]*(j[i]-jbar)
    demVote[i] ~ dnorm(mu.y[i],tau)
  }

  alpha ~ dnorm(0, 0.001) ## prior on alpha; parameters just made up for illustration
  sigma ~ dunif(0,20) ## prior on standard deviation
  tau <- pow(sigma,-2) ## convert to precision

  ## hierarchical model for each state’s intercept & slope
  for (p in 1:120) {
     beta0[p] ~ dnorm(0, tau0)
     beta1[p] ~ dnorm(mu1, tau1)
  }

  ## Priors on hierarchical components; parameters just made up for illustration
  mu1 ~ dnorm(0, 0.001)
  sigma0 ~ dunif(0,20)
  sigma1 ~ dunif(0,20)
  tau0 <- pow(sigma0,-2)
  tau1 <- pow(sigma1,-2)
}
A very useful resource for hierarchical models, including some "tricks" to speed up convergence, is Gelman and Hill.