将 glmer(二项式)翻译成 jags 以包含相关的随机效应(时间)

Translating glmer (binomial) into jags to include a correlated random effect (time)

上下文:

我有一个 12 项风险评估,其中每个人的评分为 0-4(4 为最高风险)。可以对每个人进行多次风险评估(最多 = 19,但大多数只有不到 5 次测量)。 风险的基线水平因人而异,因此我正在寻找随机截距模型,但也需要反映风险的动态性质,即添加 'time' 作为随机系数。

结果是二进制的:

最终我想要做的是根据其他人(具有相同特征的人)的评估历史、背景因素和可能随时间变化的因素来预测一个人将来是否会冒犯。

目标:

我希望通过添加时变(级别 1)和时不变(级别 2)预测变量来添加到我的 'basic' 模型中:

问题:

基本上是翻译之一。这是我的 'basic' 基于时间的模型和使用 lme4 的风险评估中的 12 个项目:

Basic_Model1 <- glmer(BinaryResponse ~ item1 + item2 + item3 + ... + item12 + time + (1+time|individual), data=data, family=binomial)

这是我将其转化为 BUGS 模型的尝试:

# the number of Risk Assessments = 552
N <-nrow(data)                                                            

# number of Individuals (individual previously specified) = 88
J <- length(unique(Individual))                                           

# the 12 items (previously specified)
Z <- cbind(item1, item2, item3, item4, ... item12)                        

# number of columns = number of predictors, will increase as model enhanced
K <- ncol(Z)                                                              

## Store all data needed for the model in a list

jags.data1 <- list(y = FO.bin, Individual =Individual, time=time, Z=Z, N=N, J=J, K=K)                    

model1 <- function() {
    for (i in 1:N) {
    y[i] ~ dbern(p[i])
    logit(p[i]) <- a[Individual[i]] + b*time[i] 
  }
  
  for (j in 1:J) {
    a[j] ~ dnorm(a.hat[j],tau.a)
    a.hat[j]<-mu.a + inprod(g[],Z[j,])
  }
  b ~ dnorm(0,.0001)
  tau.a<-pow(sigma.a,-2)
  sigma.a ~ dunif(0,100)
  
  mu.a ~ dnorm (0,.0001)
  for(k in 1:K) {
    g[k]~dnorm(0,.0001)
  }
}

write.model(model1, "Model_1.bug") 

查看输出,我的直觉是我没有为时间添加可变系数,到目前为止我所做的只是相当于

Basic_Model2 <- glmer(BinaryResponse ~ item1 + item2 + item3 + ... + item12 + time + (1|individual), data=data, family=binomial)

如何调整我的 BUGS 代码以将时间反映为可变系数,即 Basic_Model1?

根据我设法找到的例子,我知道我需要在J循环中做一个额外的规范,以便我可以监视U[j],并且需要更改第二部分涉及时间的 logit 语句,但它已经达到了我见树不见林的地步!

我希望比我更专业的人能给我指明正确的方向。最终,我希望通过添加额外的 1 级和 2 级预测变量来扩展模型。使用 lme4 查看这些内容后,我预计必须指定交互跨级交互,因此我正在寻找一种足够灵活的方法来以这种方式扩展。我是编码新手,所以请多多关照!

对于这种情况,您可以暂时使用自回归高斯模型 (CAR)。由于您的标签是 winbugs(或 openbugs),您可以按如下方式使用函数 car.normal。此代码需要适应您的数据集!

数据

y 应该是一个矩阵,观察值在行中,时间在列中。如果每个 i 的次数不同,只需输入 NA 个值。
您还需要定义时间过程的参数。这是具有权重的邻域矩阵。对不起,但我不完全记得如何创建它......对于一阶自回归,这应该是这样的:

jags.data1 <- list(
    # Number of time points
    sumNumNeigh.tm = 14, 
    # Adjacency but I do not remember how it works
    adj.tm = c(2, 1, 3, 2, 4, 3, 5, 4, 6, 5, 7, 6, 8, 7), 
    # Number of neighbours to look at for order 1
    num.tm = c(1, 2, 2, 2, 2, 2, 2, 1),
    # Matrix of data ind / time
    y = FO.bin, 
    # You other parameters
    Individual =Individual, Z=Z, N=N, J=J, K=K)     

型号

model1 <- function() {
    for (i in 1:N) {
      for (t in 1:T) {
    y[i,t] ~ dbern(p[i,t])
    # logit(p[i]) <- a[Individual[i]] + b*time[i] 
    logit(p[i,t]) <- a[Individual[i]] + b*U[t] 
  }}

  # intrinsic CAR prior on temporal random effects
  U[1:T] ~ car.normal(adj.tm[], weights.tm[], num.tm[], prec.nu)
  for(k in 1:sumNumNeigh.tm) {weights.tm[k] <- 1 }
  # prior on precison of temporal random effects
  prec.nu ~ dgamma(0.5, 0.0005)       
  # conditional variance of temporal random effects
  sigma2.nu <- 1/prec.nu                                


  for (j in 1:J) {
    a[j] ~ dnorm(a.hat[j],tau.a)
    a.hat[j]<-mu.a + inprod(g[],Z[j,])
  }
  b ~ dnorm(0,.0001)
  tau.a<-pow(sigma.a,-2)
  sigma.a ~ dunif(0,100)

  mu.a ~ dnorm (0,.0001)
  for(k in 1:K) {
    g[k]~dnorm(0,.0001)
  }
}

供您参考,使用 JAGS,您需要使用 dmnorm.

自己编写 CAR 模型