Context:
I have a 12 item risk assessment where individuals are given a rating from 0-4 (4 being the highest risk). The risk assessment can be done multiple times for each individual (max = 19, but most only have less than 5 measurements). The baseline level of risk varies by individual so I am looking for a random intercepts model, but also need to reflect the dynamic nature of the risk ie adding 'time' as a random coefficient.
The outcome is binary:
Ultimately what I am essentially looking to do is to predict whether an individual will offend in the future, based on other's (who share the same characteristics) assessment history, contextual factors and factors which may change over time.
Goal:
I wish to add to my 'basic' model by adding time-varying (level 1) and time-invariant (level 2) predictors:
Time varying include dummy variables around the criminal justice process such as non-compliance, going to court and spending time in custody. These are reflected as being an 'event' which has occurred in the period between assessments
Time invariant includes dummy variables such as being female, being non-White, and continuous predictors such as age at time of first offence I've managed to set this up OK using lmer4 and have some potentially interesting results from adding the level 1 and level 2 predictors including where there are interactions and cross-interactions. However, the complexity of the enhanced models is throwing up all kinds of warning messages including ones about failing to converge. I therefore feel that it would be appropriate to switch to a Bayesian framework using Rjags so that I can feel more confident about my findings.
The Problem:
Basically it is one of translation. This is my 'basic' model based on time and the 12 items in the risk assessment using lme4:
Basic_Model1 <- glmer(BinaryResponse ~ item1 + item2 + item3 + ... + item12 + time + (1+time|individual), data=data, family=binomial)
This is my attempt to translate this into a BUGS model:
# 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")
Looking at the output, my gut feeling is that I've not added the varying coefficient for time and that what I have done so far is only the equivalent of
Basic_Model2 <- glmer(BinaryResponse ~ item1 + item2 + item3 + ... + item12 + time + (1|individual), data=data, family=binomial)
How do I tweak my BUGS code to reflect time as a varying co-efficient ie Basic_Model1 ?
Based on the examples I have managed to find, I know that I need to make an additional specification in the J loop so that I can monitor the U[j], and there is a need to change the second part of the logit statement involving time, but its got to the point where I can't see the wood for the trees!
I'm hoping that someone with a lot more expertise than me can point me in the right direction. Ultimately I am looking to expand the model by adding additional level 1 and level 2 predictors. Having looked at these using lme4, I anticipate having to specify interactions cross-level interactions, so I am looking for an approach which is flexible enough to expand in this way. I'm very new to coding so please be gentle with me!
For that kind of case you can use an autoregressive gaussian model (CAR) for the time. As your tag is winbugs (or openbugs), you can use function car.normal as follows. This code needs to be adapted to your dataset !
y should be a matrix with observations in line and time in columns. If you do not have same number of time for each i, just put NA values.
You also need to define the parameters of the temporal process. This is the matrix of neighborhood with the weights. I am sorry, but I do not totally remember how to create it... For autoregressive of order one, this should be something like:
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)
}
}
For your information, with JAGS, you would need to code yourself the CAR model using dmnorm.