语境:
我有一个 12 项风险评估,其中个人的评分为 0-4(4 是最高风险)。可以对每个人进行多次风险评估(最多 = 19,但大多数只有少于 5 次测量)。风险的基线水平因人而异,因此我正在寻找一个随机截距模型,但还需要反映风险的动态性质,即添加“时间”作为随机系数。
结果是二进制的:
- 在测量级别发生的进一步犯罪 (FO.bin),这意味着我基本上是在查看 12 项中的一项或多项发生了哪些动态变化,以及它们如何导致个人在此期间再次犯罪测量之间
最终,我主要想做的是根据其他人(具有相同特征的人)的评估历史、背景因素和可能随时间变化的因素来预测个人将来是否会犯罪。
目标:
我希望通过添加时变(级别 1)和时不变(级别 2)预测器来添加到我的“基本”模型中:
时变包括围绕刑事司法程序的虚拟变量,例如不合规、上法庭和在押时间。这些被反映为在两次评估之间发生的“事件”
时不变包括虚拟变量,例如女性、非白人,以及连续预测变量,例如首次进攻时的年龄2 个预测变量,包括存在交互和交叉交互的位置。然而,增强模型的复杂性正在引发各种警告消息,包括关于未能收敛的警告消息。因此,我觉得使用 Rjags 切换到贝叶斯框架是合适的,这样我就可以对我的发现更有信心。
问题:
基本上它是一种翻译。这是我基于时间和使用 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 查看了这些内容后,我预计必须指定交互跨级别交互,因此我正在寻找一种足够灵活的方法来以这种方式扩展。我对编码很陌生,所以请对我温柔一点!