4
Я использую R2jags и CODA для получения некоторых диагностических статистических данных в моих сетях MCMC, но у меня возникают проблемы. Я хочу, чтобы запустить MCMC следующим образом:Ошибка в векторе («list», n.chains): недопустимый аргумент «длина»
modelfit <- jags(data=jags.data, inits=jags.inits, model.params, n.iter = 100000,
model.file=jags.model, model.params)
Ошибка заключается в следующем:
Error in vector("list", n.chains) : invalid 'length' argument
Я использую RStudio 0.97.551 и R версии 3.0.0 (2013-04-03).
Я ценю любую помощь!
Это мой R сценарий:
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+j-i])/(1+exp(mu1+theta1[i]+phi1[j]+psi1[I+j-i]))
thetapred[i,j] <- exp(mu2+theta2[i]+phi2[j]+psi2[I+j-i])/(1+exp(mu2+theta2[i]+phi2[j]+psi2[I+j-i]))
etapred[i,j] <- exp(mu3+theta3[i]+phi3[j]+psi3[I+j-i])/(1+exp(mu3+theta3[i]+phi3[j]+psi3[I+j-i]))
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')