I use R2jags and CODA to get some diagnostic statistics for my MCMC networks, but I'm having problems. I want to run MCMC as follows:
modelfit <- jags(data=jags.data, inits=jags.inits, model.params, n.iter = 100000, model.file=jags.model, model.params)
The error is as follows:
Error in vector("list", n.chains) : invalid 'length' argument
I am using RStudio 0.97.551 and R version 3.0.0 (2013-04-03).
I appreciate any help!
This is my R script:
setwd("C:\\") y1 <-c(1538727, 1444672, 1206999, 1002960, 744597, 390301, 1640130, 1472255, 1383947, 1109395, 984775, 697701, 1769569, 1573498, 1489025, 1351284, 1111397, 935166, 1747764, 1790841, 1626852, 1407388, 1284583, 995236, 1676555, 1787181, 1655400, 1527122, 1421772, 1309989, 1561922, 1643467, 1598855, 1570645, 1495999, 1319439, 1456258, 1561892, 1567872, 1555237, 1551579, 1532222, 1243436, 1387943, 1436659, 1511134, 1549578, 1539580) y2 <- c(2634569, 3031916, 3138776, 2875868, 2495888, 1886174, 2148776, 2567507, 2747428, 2696199, 2593985, 2138303, 1662296, 2224336, 2489723, 2698322, 2655746, 2450716, 1304387, 1734318, 2180203, 2396749, 2629088, 2555934, 1087351, 1380119, 1616309, 2109287, 2408800, 2369855, 821642, 1041702, 1221283, 1661647, 2098345, 2426842, 708327, 873092, 952245, 1237084, 1628334, 2123709, 549763, 666699, 774205, 981393, 1243888, 1538431) y3 <- c(1245931, 1664176, 1659375, 2313647, 3850196, 4254634, 825634, 1293382, 1454776, 1736181, 2596719, 3655532, 554953, 901957, 1186747, 1490664, 2083400, 2738988, 335824, 630232, 847486, 1239538, 1702256, 2296941, 218213, 373786, 555286, 907876, 1397221, 2005940, 143202, 237344, 344229, 594993, 1012777, 1510283, 121187, 151070, 219731, 351040, 650930, 1157146, 87211, 120279, 140551, 226530, 393887, 733699) n <- c(5862309, 6673625, 6534802, 6942747, 8329067, 8152696, 5049199, 5913474, 6268113, 6253757, 7298375, 8260640, 4319559, 5245545, 5840408, 6306245, 6785242, 7492958, 3588778, 4553684, 5259609, 5813653, 6517271, 7001560, 3105173, 3797508, 4271831, 5180290, 6086716, 7002991, 2591140, 3063506, 3428373, 4305319, 5326889, 6217360, 2329398, 2661972, 2886111, 3418403, 4327922, 5565798, 1906676, 2224544, 2444586, 2864892, 3473404, 4362648) A <- c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8) P <- c(1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6) C <- c(8, 9, 10, 11, 12, 13, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 5, 6, 7, 8, 9, 10, 4, 5, 6, 7, 8, 9, 3, 4, 5, 6, 7, 8, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6) W <- 48 I <- 8 J <- 6 JJ <- 8 L <- 13 LL <- 15 jags.model<- function() { ## Loop over all observations for(w in 1: W){ m[w] <- n[w]-y1[w] h[w] <- n[w]-y1[w]-y2[w] y1[w] ~ dbin(delta[w],n[w]) y2[w] ~ dbin(theta[w],m[w]) y3[w] ~ dbin(eta[w],h[w]) y4[w] <- n[w]-y1[w]-y2[w]-y3[w] ## Define model logit(delta[w]) <- mu1+theta1[A[w]]+phi1[P[w]]+psi1[C[w]] logit(theta[w]) <- mu2+theta2[A[w]]+phi2[P[w]]+psi2[C[w]] logit(eta[w]) <- mu3+theta3[A[w]]+phi3[P[w]]+psi3[C[w]] ## Define probabilities p1[w] <- delta[w] p2[w] <- theta[w]*(1-delta[w]) p3[w] <- eta[w]*(1-delta[w])*(1-theta[w]) p4[w] <- (1-delta[w])*(1-theta[w])*(1-eta[w]) L1[w] <- log(p1[w]/(p2[w]+p3[w]+p4[w])) L2[w] <- log(p2[w]/(p3[w]+p4[w])) L3[w] <- log(p3[w]/p4[w]) ## Likelihood log.like[w] <- y1[w]*log(delta[w]) + (n[w]-y1[w])*log(1-delta[w]) + y2[w]*log(theta[w]) + (n[w]-y1[w]-y2[w])*log(1-theta[w]) + y3[w]*log(eta[w]) + (n[w]-y1[w]-y2[w]-y3[w])*log(1-eta[w]) + logfact(n[w]) - logfact(y1[w]) - logfact(y2[w]) - logfact(y3[w]) - logfact(n[w]-y1[w]-y2[w]-y3[w]) } ## Prior model for global mean effects mu1~ dnorm(0,1000) mu2~ dnorm(0,1000) mu3~ dnorm(0,1000) ## Autoregressive prior model for P effects phi1mean[1] <- 0.0 phi1prec[1] <- tauphi1*1.0E-6 phi1mean[2] <- 0.0 phi1prec[2] <- tauphi1*1.0E-6 phi2mean[1] <- 0.0 phi2prec[1] <- tauphi2*1.0E-6 phi2mean[2] <- 0.0 phi2prec[2] <- tauphi2*1.0E-6 phi3mean[1] <- 0.0 phi3prec[1] <- tauphi3*1.0E-6 phi3mean[2] <- 0.0 phi3prec[2] <- tauphi3*1.0E-6 phi4mean[1] <- 0.0 phi4prec[1] <- tauphi4*1.0E-6 phi4mean[2] <- 0.0 phi4prec[2] <- tauphi4*1.0E-6 for (j in 3:JJ) { phi1mean[j] <- 2*phi1[j-1]-phi1[j-2] phi1prec[j] <- tauphi1 phi2mean[j] <- 2*phi2[j-1]-phi2[j-2] phi2prec[j] <- tauphi2 phi3mean[j] <- 2*phi3[j-1]-phi3[j-2] phi3prec[j] <- tauphi3 phi4mean[j] <- 2*phi4[j-1]-phi4[j-2] phi4prec[j] <- tauphi4 } # Sampling P effects and subtracting mean for observed P (J) for (j in 1:JJ) { phi1[j] ~ dnorm(phi1mean[j],phi1prec[j]) phi2[j] ~ dnorm(phi2mean[j],phi2prec[j]) phi3[j] ~ dnorm(phi3mean[j],phi3prec[j]) phi4[j] ~ dnorm(phi4mean[j],phi4prec[j]) phi1c[j] <- phi1[j]-mean(phi1[1:J]) phi2c[j] <- phi2[j]-mean(phi2[1:J]) phi3c[j] <- phi3[j]-mean(phi3[1:J]) phi4c[j] <- phi4[j]-mean(phi4[1:J]) } # Hyperpriors for the precision parameters tauphi1 ~ dgamma(1.0E-1,1.0E-1) tauphi2 ~ dgamma(1.0E-1,1.0E-1) tauphi3 ~ dgamma(1.0E-1,1.0E-1) tauphi4 ~ dgamma(1.0E-1,1.0E-1) sigmaphi1 <- 1/sqrt(tauphi1) sigmaphi2 <- 1/sqrt(tauphi2) sigmaphi3 <- 1/sqrt(tauphi3) sigmaphi4 <- 1/sqrt(tauphi4) ## Autoregressive prior model for C effects psi1mean[1] <- 0.0 psi1prec[1] <- taupsi1*1.0E-6 psi1mean[2] <- 0.0 psi1prec[2] <- taupsi1*1.0E-6 psi2mean[1] <- 0.0 psi2prec[1] <- taupsi2*1.0E-6 psi2mean[2] <- 0.0 psi2prec[2] <- taupsi2*1.0E-6 psi3mean[1] <- 0.0 psi3prec[1] <- taupsi3*1.0E-6 psi3mean[2] <- 0.0 psi3prec[2] <- taupsi3*1.0E-6 psi4mean[1] <- 0.0 psi4prec[1] <- taupsi4*1.0E-6 psi4mean[2] <- 0.0 psi4prec[2] <- taupsi4*1.0E-6 for (l in 3:LL) { psi1mean[l] <- 2*psi1[l-1]-psi1[l-2] psi1prec[l] <- taupsi1 psi2mean[l] <- 2*psi2[l-1]-psi2[l-2] psi2prec[l] <- taupsi2 psi3mean[l] <- 2*psi3[l-1]-psi3[l-2] psi3prec[l] <- taupsi3 psi4mean[l] <- 2*psi4[l-1]-psi4[l-2] psi4prec[l] <- taupsi4 } # Sampling C effects and subtracting mean for observed C (L) for (l in 1:LL) { psi1[l] ~ dnorm(psi1mean[l],psi1prec[l]) psi2[l] ~ dnorm(psi2mean[l],psi2prec[l]) psi3[l] ~ dnorm(psi3mean[l],psi3prec[l]) psi4[l] ~ dnorm(psi4mean[l],psi4prec[l]) psi1c[l] <- psi1[l]-mean(psi1[1:L]) psi2c[l] <- psi2[l]-mean(psi2[1:L]) psi3c[l] <- psi3[l]-mean(psi3[1:L]) psi4c[l] <- psi4[l]-mean(psi4[1:L]) } # HyPpriors for the precision parameters taupsi1 ~ dgamma(1.0E-1,1.0E-1) taupsi2 ~ dgamma(1.0E-1,1.0E-1) taupsi3 ~ dgamma(1.0E-1,1.0E-1) taupsi4 ~ dgamma(1.0E-1,1.0E-1) sigmapsi1 <- 1/sqrt(taupsi1) sigmapsi2 <- 1/sqrt(taupsi2) sigmapsi3 <- 1/sqrt(taupsi3) sigmapsi4 <- 1/sqrt(taupsi4) ## Autoregressive prior model for A effects theta1mean[1] <- 0.0 theta1prec[1] <- tautheta1*1.0E-6 theta1mean[2] <- 0.0 theta1prec[2] <- tautheta1*1.0E-6 theta2mean[1] <- 0.0 theta2prec[1] <- tautheta2*1.0E-6 theta2mean[2] <- 0.0 theta2prec[2] <- tautheta2*1.0E-6 theta3mean[1] <- 0.0 theta3prec[1] <- tautheta3*1.0E-6 theta3mean[2] <- 0.0 theta3prec[2] <- tautheta3*1.0E-6 theta4mean[1] <- 0.0 theta4prec[1] <- tautheta4*1.0E-6 theta4mean[2] <- 0.0 theta4prec[2] <- tautheta4*1.0E-6 for (i in 3:I) { theta1mean[i] <- 2*theta1[i-1]-theta1[i-2] theta1prec[i] <- tautheta1 theta2mean[i] <- 2*theta2[i-1]-theta2[i-2] theta2prec[i] <- tautheta2 theta3mean[i] <- 2*theta3[i-1]-theta3[i-2] theta3prec[i] <- tautheta3 theta4mean[i] <- 2*theta4[i-1]-theta4[i-2] theta4prec[i] <- tautheta4 } # Sampling A effects for (i in 1:I) { theta1[i] ~ dnorm(theta1mean[i],theta1prec[i]) theta2[i] ~ dnorm(theta2mean[i],theta2prec[i]) theta3[i] ~ dnorm(theta3mean[i],theta3prec[i]) theta4[i] ~ dnorm(theta4mean[i],theta4prec[i]) } # Hyperpriors for the precision parameters tautheta1 ~ dgamma(1.0E-1,1.0E-1) tautheta2 ~ dgamma(1.0E-1,1.0E-1) tautheta3 ~ dgamma(1.0E-1,1.0E-1) tautheta4 ~ dgamma(1.0E-1,1.0E-1) sigmatheta1 <- 1/sqrt(tautheta1) sigmatheta2 <- 1/sqrt(tautheta2) sigmatheta3 <- 1/sqrt(tautheta3) sigmatheta4 <- 1/sqrt(tautheta4) # Removing linear trends from A effects for (i in 1:I) { ivec1[i] <- i-(I+1)/2 aivec1[i] <- ivec1[i]*theta1[i] theta1c[i] <- theta1[i]-ivec1[i]*sum(aivec1[])/(I*(I+1)*(I-1)/12) ivec2[i] <- i-(I+1)/2 aivec2[i] <- ivec2[i]*theta2[i] theta2c[i] <- theta2[i]-ivec2[i]*sum(aivec2[])/(I*(I+1)*(I-1)/12) ivec3[i] <- i-(I+1)/2 aivec3[i] <- ivec3[i]*theta3[i] theta3c[i] <- theta3[i]-ivec3[i]*sum(aivec3[])/(I*(I+1)*(I-1)/12) ivec4[i] <- i-(I+1)/2 aivec4[i] <- ivec4[i]*theta4[i] theta4c[i] <- theta4[i]-ivec4[i]*sum(aivec4[])/(I*(I+1)*(I-1)/12) } ## Loop over A and P ## Computing fitted and projected probabilities for (i in 1:I) { for (j in 1:JJ) { deltapred[i,j] <- exp(mu1+theta1[i]+phi1[j]+psi1[I+ji])/(1+exp(mu1+theta1[i]+phi1[j]+psi1[I+ji])) thetapred[i,j] <- exp(mu2+theta2[i]+phi2[j]+psi2[I+ji])/(1+exp(mu2+theta2[i]+phi2[j]+psi2[I+ji])) etapred[i,j] <- exp(mu3+theta3[i]+phi3[j]+psi3[I+ji])/(1+exp(mu3+theta3[i]+phi3[j]+psi3[I+ji])) p1p[i,j] <- deltapred[i,j] p2p[i,j] <- thetapred[i,j]*(1-deltapred[i,j]) p3p[i,j] <- etapred[i,j]*(1-deltapred[i,j])*(1-thetapred[i,j]) p4p[i,j] <- (1-deltapred[i,j])*(1-thetapred[i,j])*(1-etapred[i,j]) } } } jags.data<-list("y1","y2","y3","n","A","P","C","W","I","J","JJ","L","LL") jags.inits <- function(){ list("tauphi1" <-1,"tauphi2" <-1,"tauphi3" <-1,"tauphi4" <-1,"taupsi1" <-1, "taupsi2" <-1,"taupsi3" <-1,"taupsi4" <-1,"tautheta1" <-1,"tautheta2" <-1, "tautheta3" <-1,"tautheta4" <-1,"mu1" <-0,"mu2" <-0,"mu3" <-0, "theta1" <-c(0, 0, 0, 0, 0, 0, 0, 0), "theta2" <-c(0, 0, 0, 0, 0, 0, 0, 0), "theta3" <-c(0, 0, 0, 0, 0, 0, 0, 0), "theta4" <-c(0, 0, 0, 0, 0, 0, 0, 0), "phi1" <-c(0, 0, 0, 0, 0, 0, 0, 0), "phi2" <-c(0, 0, 0, 0, 0, 0, 0, 0), "phi3" <-c(0, 0, 0, 0, 0, 0, 0, 0), "phi4" <-c(0, 0, 0, 0, 0, 0, 0, 0), "psi1" <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), "psi2" <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), "psi3" <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), "psi4" <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) } library(R2jags) library(coda) model.params <- c('mu1', 'mu2', 'mu3', 'theta1', 'theta2', 'theta3', 'phi1', 'phi2', 'phi3', 'psi1', 'psi2', 'psi3', 'p1p', 'p2p', 'p3p', 'p4p')