I'm fairly new to using JAGS and I'm trying to get a growth curve model over time with two levels of variable (drug test group + sex). I'm running the following code but end up with : "Error in jags.model("GCM.txt", data = jagsData, n.chains = chains, n.adapt = adaptation, : Error in node betas[1,1:4] Invalid parent values" when creating the jags model object.
# Creating a matrix with only weights:
# We `unlist' all N rows of selected columns from the data set,
# then transform these values into numeric entries of a matrix
data <- as.matrix(Test_wide[,9:124])
# Creating a variable that indicates the number of time points
nrT <- 116
# Saving X1, X2, X3 and X4 as separate variables (same unlisting process as performed above)
grouping <- as.matrix(Test_wide[,5:8])
X1 <- grouping[,1] # 1 for group1
X2 <- grouping[,2] # 1 for group2
X3 <- grouping[,3] # 1 for Male
X4 <- grouping[,4] # 1 for Female
# Creating a time vector for the measurement waves
time <- c(2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,
31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,48,50,51,52,53,54,55,56,57,58,59,60,
61,62,63,64,65,66,67,68,69,70,71,72,73,74,82,83,84,87,
91,92,93,94,95,96,97,98,99,100,101,102,103,104,106,107,118,121,122,123,124,125,126,127,128,129,130,
131,132,133,134,135,136,137,138,139,144,151,152,153,154)
# time vector (T) based on the time values of our measures
# Now we have all the data needed to be passed to JAGS
# Creating a list of all the variables that we created above
jagsData <- list("Y"=data,"X1"=X1,"X2"=X2,"X3"=X3,"X4"=X4,"N"=N,"nrT"=nrT,"time"=time)
## ----Step 2: Specify the growth curve model with linear trend, highlight=FALSE----
# Step 2
# Specifying the growth curve model
LinearGrowthCurve = cat("
model {
# Loop over individuals
for (i in 1:N) {
# Loop over time
for (t in 1:nrT) {
# Specifying likelihood function, corresponding to:
Y[i,t] ~ dnorm(betas[i,1] + betas[i,2] * time[t] + betas[i,3] + betas[i,4] * time[t], 1/sdLevel1Error^2)
}
# Specifying the multivariate distribution of intercepts and slopes
betas[i,1:4] ~ dmnorm.vcov(Level2MeanVector[i,1:4], interpersonCovMatrix[1:4,1:4])
# The mean of the intercept is modeled as a function of positivity group membership
Level2MeanVector[i,1] <- Med11PInt + betaAviPInt*X1[i] + betaConPInt*X2[i]
Level2MeanVector[i,2] <- Med11PSlope + betaAviPSlope*X1[i] + betaConPSlope*X2[i]
Level2MeanVector[i,3] <- Med12PInt + betaMalePInt*X3[i] + betaFemalePInt*X4[i]
Level2MeanVector[i,4] <- Med12PSlope + betaMalePSlope*X3[i] + betaFemalePSlope*X4[i]
}
# Specifying prior distributions variable
##sderror
sdLevel1Error ~ dunif(0,100)
##statut XP
Med11PInt ~ dnorm(0,0.01)
Med11PSlope ~ dnorm(0,0.01)
betaAviPInt ~ dnorm(0,0.01)
betaConPInt ~ dnorm(0,0.01)
betaAviPSlope ~ dnorm(0,0.01)
betaConPSlope ~ dnorm(0,0.01)
sd11Intercept ~ dunif(0,100)
sd11Slope ~ dunif(0,100)
##Sex
Med12PInt ~ dnorm(0,0.01)
Med12PSlope ~ dnorm(0,0.01)
betaMalePInt ~ dnorm(0,0.01)
betaFemalePInt ~ dnorm(0,0.01)
betaMalePSlope ~ dnorm(0,0.01)
betaFemalePSlope ~ dnorm(0,0.01)
sd12Intercept ~ dunif(0,100)
sd12Slope ~ dunif(0,100)
##Correlation matrice
corr11IntSlope ~ dunif(-1,1)
corr21IntSlope ~ dunif(-1,1)
corr12IntSlope ~ dunif(-1,1)
corr22IntSlope ~ dunif(-1,1)
# Transforming model parameters
# 1. Defining the elements of the level-2 covariance matrix
interpersonCovMatrix[1,1] <- sd11Intercept * sd11Intercept
interpersonCovMatrix[1,2] <- corr11IntSlope * sd11Intercept* sd11Slope
interpersonCovMatrix[1,3] <- sd11Intercept * sd12Intercept
interpersonCovMatrix[1,4] <- corr12IntSlope * sd11Intercept * sd12Slope
interpersonCovMatrix[2,1] <- interpersonCovMatrix[1,2]
interpersonCovMatrix[2,2] <- sd11Slope * sd11Slope
interpersonCovMatrix[2,3] <- corr21IntSlope * sd11Slope * sd12Intercept
interpersonCovMatrix[2,4] <- sd11Slope * sd12Slope
interpersonCovMatrix[3,1] <- interpersonCovMatrix[1,3]
interpersonCovMatrix[3,2] <- interpersonCovMatrix[2,3]
interpersonCovMatrix[3,3] <- sd12Intercept * sd12Intercept
interpersonCovMatrix[3,4] <- corr22IntSlope * sd12Intercept* sd12Slope
interpersonCovMatrix[4,1] <- interpersonCovMatrix[1,4]
interpersonCovMatrix[4,2] <- interpersonCovMatrix[2,4]
interpersonCovMatrix[4,3] <- interpersonCovMatrix[3,4]
interpersonCovMatrix[4,4] <- sd12Slope * sd12Slope
# 2. Creating status XP variables representing:
## 2.1. Group1 intercept
AviPInt <- Med11PInt + betaAviPInt
## 2.2. Group2 intercept
ConPInt <- Med11PInt + betaConPInt
## 2.3. Group1 slope
AviPSlope <- Med11PSlope + betaAviPSlope
## 2.4. Group2 slope
ConPSlope <- Med11PSlope + betaConPSlope
#3 Creating sex variables representing:
## 3.1. Male intercept
MalePInt <- Med12PInt + betaMalePInt
## 3.2. Female intercept
FemalePInt <- Med12PInt + betaFemalePInt
## 3.3. Male slope
MalePSlope <- Med12PSlope + betaMalePSlope
## 3.4. Female slope
FemalePSlope <- Med12PSlope + betaFemalePSlope
# 4. Contrast terms between Groups and Sex intercepts and slopes
AviConPInt <- betaAviPInt - betaConPInt
MaleFemalePInt <- betaMalePInt - betaFemalePInt
AviConPSlope <- betaAviPSlope - betaConPSlope
MaleFemalePSlope <- betaMalePSlope - betaFemalePSlope
}
",file = "GCM.txt")
## ----Step 3: Collect model parameters, set sampler, tidy = TRUE----------
# Step 3
# Monitoring parameters
parameters <- c("Med11PInt","betaAviPInt", "betaConPInt", "Med12PInt","betaMalePInt", "betaFemalePInt",
"Med11PSlope", "betaAviPSlope", "betaConPSlope", "Med12PSlope", "betaMalePSlope", "betaFemalePSlope",
"sd11Intercept", "sd11Slope",
"sdLevel1Error","betas",
"sd12Intercept", "sd12Slope",
"corr11IntSlope", "corr12IntSlope", "corr21IntSlope", "corr22IntSlope",
"AviPInt","ConPInt","AviPSlope","ConPSlope",
"MalePInt", "FemalePInt", "MalePSlope", "FemalePSlope",
"AviConPInt","MaleFemalePInt",
"AviConPSlope","MaleFemalePSlope")
# Specifying sampler settings
adaptation <- 200 # Number of steps to "tune" the samplers
chains <- 2 # Re-start the exploration "chains" number of times
burnin <- 1000 # Number of steps to get rid of the influence of initial values
thinning <- 10
# Defining the number of samples drawn from the posterior in each chain
postSamples <- 3000
# Computing the number of posterior samples needed per chain for JAGS
nrOfIter <- ceiling((postSamples*thinning)/chains)
## ----echo=FALSE, eval=TRUE, warning=FALSE, message=FALSE, error=FALSE, include=FALSE----
#, results=hide,>>=
# fixing the random seed for reproducibility
fixedinits<- list(list(.RNG.seed=5,.RNG.name="base::Wichmann-Hill"),list(.RNG.seed=6,.RNG.name="base::Wichmann-Hill"))
## ----Step 4: Call JAGS and fit the model, tidy = TRUE, eval=TRUE---------
# Step 4
# Loading the rjags package
library(r2jags)
# Creating JAGS model object
jagsModel<-jags.model("GCM.txt",data=jagsData,n.chains=chains,n.adapt=adaptation,inits=fixedinits)
# Running burn-in iterations
update(jagsModel,n.iter=burnin)
# Drawing posterior samples
codaSamples<-coda.samples(jagsModel,variable.names=parameters,n.iter=nrOfIter, thin = thinning,seed=5)
The data are here: text.
Any help you can provide is greatly appreciated,
Thanks for your time, Marc Silvestre