[R] R2WinBUGS Error

[Liu, J.]

Hi,

I'm trying to run R2WinBUGS using the R code below (Thinkpad Yoga 260,
Win8, system x86_64, mingw32, R version 3.3.1). It worked fine for several
times but then one error began to pop up in every run: command #Bugs:set
cannot be executed (is greyed out). I've been trying for more than one week
but still can't figure out where is the problem. It would be great if
someone could help me with this. Thanks in advance!

Kind regards,
JLiu

Here's the code:

sink("mod1.txt")
cat("
    model {

    for( k in 1 : n ) {
    for( i in 1 : n - 1 ) {
    for( j in i + 1 : n ) {
    y[k , i , j] ~ dbern(p[k , i , j])
    y[k , j , i] ~ dbern(p[k , j , i])
    logit(p[k , i , j]) <- mu + a[i] + b[j] + g[k] + ab[i , j] + ag[i , k]
+ ag[j , k]
    logit(p[k , j , i]) <- mu + a[j] + b[i] + g[k] + ab[j , i] + ag[j , k]
+ ag[i , k]
    }
    }
    }

    # Compute difference ab[1,2] - ab[2,1] for Figure 3 of the paper
    dif12 <- ab[1 , 2] - ab[2 , 1]

    # Prior for the overall mean effect mu
    mu ~ dnorm(0, 1)

    # Tri-normal prior for actor, partner, rater effects (a, b, g)
    # with zero-sum constraints
    for( i in 1 : n ) {
    a[i] <- a1[i , 1] - mean(a1[ , 1])
    b[i] <- a1[i , 2] - mean(a1[ , 2])
    g[i] <- a1[i , 3] - mean(a1[ , 3])

    # without zero-sum constraints will make the code run faster
    # for( i in 1 : n ) {
    # a[i] <- a1[i,1]
    # b[i] <- a1[i,2]
    # g[i] <- a1[i,3]

    a1[i , 1:3] ~ dmnorm(zero[1:3], S1[ , ])
    }

    # Wishart prior for precision matrix S1 = Sigma_1-inverse
    # degrees of freedom nu = 3
    # Om1 = nu*Identity Matrix is provided in the data list
    S1[1:3 , 1:3] ~ dwish(Om1[1:3 , 1:3], 3)
    Sig1[1:3 , 1:3] <- inverse(S1[1:3 , 1:3])

    # Compute the correlations and the variances for the main effects
    rho1 <- Sig1[1 , 2] / sqrt(Sig1[1 , 1] * Sig1[2 , 2])
    rho2 <- Sig1[1 , 3] / sqrt(Sig1[1 , 1] * Sig1[3 , 3])
    rho3 <- Sig1[2 , 3] / sqrt(Sig1[3 , 3] * Sig1[2 , 2])
    sig.a <- Sig1[1 , 1]
    sig.b <- Sig1[2 , 2]
    sig.g <- Sig1[3 , 3]

    # Standard deviation values are reported in the paper (Table 2)
    sd.a <- sqrt(sig.a)
    sd.b <- sqrt(sig.b)
    sd.g <- sqrt(sig.g)
    # A vector of zero's is needed for the mean values of some vectors
    for( i in 1 : 4 ) {
    zero[i] <- 0
    }

    # Priors for interaction effects ab_ij and ag_ij in equation (3) of the
paper
    for( i in 1 : n ) {
    # alpha_beta[i,i] is not defined in the model, just put  =0
    ab[i , i] <- 0

    # del[i] = personal bias of subject i in reporting his tendency to
establish friendship ties
    # The posterior distribution of del[i]'s is shown as boxplots in Figure
3 of the paper.
    # Side-by-side boxplots are created using Inference -> Compare menu in
WinBUGS

    del[i] <- ag[i , i] - ag.mean[i]
    ag.mean[i] <- (sum(ag[i , ]) - ag[i , i]) / (n - 1)
    ag[i , i] ~ dnorm(0, tau.ag)
    }

    # Generate 4 pair-wise interaction parameters as Normal_4 with
precision matrix S2
    for( i in 1 : n - 1 ) {
    for( j in i + 1 : n ) {
    ab[i , j] <- a2[i , j , 1]
    ab[j , i] <- a2[i , j , 2]
    ag[i , j] <- a2[i , j , 3]
    ag[j , i] <- a2[i , j , 4]
    a2[i , j , 1:4] ~ dmnorm(zero[1:4], S2[ , ])
    }
    }

    # Get the precision matrix of interaction parameters from the var-cov
matrix Sig_2 in the paper
    # Note that we don't use Inv-Wishart prior for this precision matrix
    S2[1:4 , 1:4] <- inverse(Sig2[1:4 , 1:4])

    # Now build the var-cov matrix Sig2 from variances and correlations
(phi's) in equation (7)

    phi1 ~ dunif(-0.99, 0.99)
    phi2 ~ dunif(-0.99, 0.99)
    phi3 ~ dunif(-0.99, 0.99)
    phi4 ~ dunif(-0.99, 0.99)

    # Compute two detrminants in equations (9-10) in the paper
    det1 <-  1 -  phi1 * phi1 - phi2 * phi2 - phi3 * phi3 + 2 * phi1 * phi2
* phi3
    det2 <-  1 -  phi1 * phi1 - 2 * phi2 * phi2 - 2 * phi3 * phi3 - phi4 *
phi4 + 4 * phi1 * phi2 * phi3
    + 4 * phi2 * phi3 * phi4 + phi1 * phi1 * phi4 * phi4 - 2 * phi2 * phi2
* phi3 * phi3
    - 2 * phi1 * phi3 * phi3 * phi4 - 2 * phi1 * phi2 * phi2 * phi4 + phi2
* phi2 * phi2 * phi2
    + phi3 * phi3 * phi3 * phi3

    # Check if determinants are positive
    cons1 <- step(det1)
    cons2 <- step(det2)

    # Zeros trick
    O1 <- 0
    O2 <- 0
    q1 <- 1 - cons1
    q2 <- 1 - cons2
    O1 ~ dbern(q1)
    O2 ~ dbern(q2)

    # Inverse Gamma priors for variance parameters sig.ab, sig.ag
    tau.ab ~ dgamma(100, 10)
    tau.ag ~ dgamma(100, 10)
    sig.ab <- 1 / tau.ab
    sig.ag <- 1 / tau.ag

    # Standard deviations sd.ab and sd.ag are reported in the paper (Table
2)
    sd.ab <- sqrt(sig.ab)
    sd.ag <- sqrt(sig.ag)
    Sig2[1 , 1] <- sig.ab
    Sig2[3 , 3] <- sig.ag
    Sig2[2 , 2] <- Sig2[1,1]
    Sig2[4 , 4] <- Sig2[3 , 3]

    # Create the off-diagonal entries of Sig2
    # If the generated set of phi1, phi2, phi3 and phi4 values don't have
both det1 > 0
    # and det2 > 0, the off-diagonal elements of Sig2 are all set equal to
zero
    Sig2[1,2] <- cons1*cons2*sig.ab*phi1
    Sig2[2 , 1] <- Sig2[1,2]
    Sig2[1 , 3] <- cons1*cons2*sqrt(sig.ab * sig.ag) * phi2
    Sig2[3 , 1] <- Sig2[1 , 3]
    Sig2[1 , 4] <- cons1*cons2*sqrt(sig.ab * sig.ag) * phi3
    Sig2[4 , 1] <- Sig2[1 , 4]
    Sig2[2 , 3] <- cons1*cons2*sqrt(sig.ab * sig.ag) * phi3
    Sig2[3 , 2] <- Sig2[2 , 3]
    Sig2[2 , 4] <- cons1*cons2*sqrt(sig.ab * sig.ag) * phi2
    Sig2[4 , 2] <-  Sig2[2 , 4]
    Sig2[3 , 4] <- cons1*cons2*sig.ag * phi4
    Sig2[4 , 3] <- Sig2[3 , 4]
    #############################################################

    # Predict the network for studying the model fit using residuals
    # for Model Selection (Section 3 of the paper)
    for( k in 1 : n ) {
    for( i in 1 : n - 1 ) {
    for( j in i + 1 : n ) {
    y.pred[k , i , j] ~ dbern(p[k , i , j])
    y.pred[k , j , i] ~ dbern(p[k , j , i])

    # How often we get wrong predictions? Compute residual e[k,i,j]
    e[k , i , j] <- abs(y[k , i , j] - y.pred[k , i , j])
    e[k , j , i] <- abs(y[k , j , i] - y.pred[k , j , i])
    }
    }
    }

    # Put diagonal residuals equal to zero. This is needed for computing
the sum() function
    for( k in 1 : n ) {
    for( i in 1 : n ) {
    e[k , i , i] <- 0
    }
    }
    # Compute epd = the total number of incorrect predictions, equation
(13) of the paper
    epd <- sum(e[,,])

    } # End of model


    ", fill = TRUE)
sink()

# read data as a whole three-dimension matrix
nr <- 21 # netowrk F
X <- "Kr"
network <- array(NA, dim=c(nr, nr, nr))
for(n in 1:nr){
  network[,,n] = as.matrix(read.table(paste0("C:/R/internship/network",X,
"/Sheet", n, ".txt")))
}

y = structure(.Data=network)
n = 21
Om1=structure(.Data=c(3,0,0,0,3,0,0,0,3), .Dim = c(3,3))
data = list(n=n, Om1=Om1, y=y)


params = c("mu", "a1", "S1", "Sig1", "a2", "S2", "sig2", "phi1", "phi2",
"phi3", "phi4","tau.ab", "tau.ag")

inits1 <- function () {list(mu = 0, tau.ab=1,tau.ag=1,
                            phi1=0.1,phi2=0.3,phi3=-0.4,phi4=0.1,
                            S1=structure(.Data=c(1,0,0,0,1,0,0,0,1), .Dim =
c(3,3)))
}

inits2 <- function () {list(mu = -5.0, tau.ab=10,tau.ag=10,
                           phi1=0,phi2=0,phi3=0.2,phi4=-0.3,
                           S1=structure(.Data=c(10,0,0,0,10,0,0,0,10), .Dim
= c(3,3)))
  }
inits3 <- function () {list(mu = 2.0, tau.ab=5,tau.ag=10,
                            phi1=0.3,phi2=-0.4,phi3=0.4,phi4=0,
                            S1=structure(.Data=c(5,0,0,0,5,0,0,0,5), .Dim =
c(3,3)))
  }

inits <- function () {list(inits1, inits2, inits3)}

nc <- 1      #number of MCMC chains to run
ni <- 4000  #number of samples for each chain
nb <- 2000   #number of samples to discard as burn-in
nt <- 1      #thinning rate, increase this to reduce autocorrelation

bugs.out <- bugs(data=data, inits=inits1, parameters.to.save=params,
model.file="mod1.txt",
                 n.chains=nc, n.iter=ni, n.burnin=nb, n.thin=nt,
debug=FALSE, DIC=TRUE,
                 bugs.directory = "C:/winBUGS/WinBUGS14",
working.directory=getwd())

	[[alternative HTML version deleted]]