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]]