2017-06-01 150 views
1

语境:翻译glmer(二项式)到尖齿包括一个相关的随机效应(时间)

我有个人在那里从0-4(4级为最高风险给出等级12项风险评估)。风险评估可以针对每个人进行多次(最多= 19,但大多数只有少于5次)。 风险的基准水平因人而异,因此我正在寻找随机截距模型,但也需要反映风险的动态性,即将“时间”添加为随机系数。

结果是二进制的:

  • 进一步违规(FO.bin)发生在测量水平,将意味着我基本上是在寻找具有一个或多个内发生的动态变化12项内容以及他们如何在测量期间对个人犯下进一步犯罪做出的贡献

最终,我基本上想要做的是预测一个ind基于其他人(具有相同特征)的评估历史,背景因素和随时间可能改变的因素,个人将在未来得罪。

目标:

我希望通过增加时变(1级)和时间不变性添加到我的 '基本' 的模式(等级2)预测:

  • 时间变化的包括围绕刑事司法程序的虚拟变量,如违规,上法庭和花时间羁押。这些被反射为一个“事件”已经发生在评估

    之间的周期
  • 时间不变包括在第一次进攻的时间虚拟变量诸如为女性,是非白色,和连续预测如年龄 我已经设法使用lmer4进行设置,并且添加了1级和2级预测变量,包括有交互作用和交互作用的地方有一些可能有趣的结果。然而,增强模型的复杂性正在抛出各种警告信息,包括未能收敛的信息。因此,我认为使用Rjags切换到贝叶斯框架是适当的,这样我可以对我的发现感到更加自信。

问题:

基本上它是一个平移。这是一个基于时间和风险评估的12个项目利用lme4我的“基本”的模式:

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],并且需要更改第二部分关于时间的逻辑陈述,但是它达到了我无法看到树木的地步!

我希望有比我更多专业知识的人可以指引我走向正确的方向。最终,我希望通过增加额外的1级和2级预测变量来扩展模型。通过使用lme4查看这些内容,我预计不得不指定交互跨级交互,所以我正在寻找一种足够灵活的方法来以这种方式进行扩展。我对编码很陌生,所以请和我一起温柔!

回答

0

对于这种情况,您可以使用自回归高斯模型(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模型。

+0

我会报告我如何继续。你的解释看起来非常清楚,在研究如何编写CAR模型时没有问题,因为这个过程正在帮助我同时改进编码。 – Helen

相关问题