R-WinBUGS code to fit the 2sSEM

We present the used code to fit Model 4 in the paper. This code can easily be adapted to fit the other models discussed in the paper. For a detailed explanation of the R2WinBUGSpackage the reader is referred to Sturtz, Ligges, and Gelman (2005); and Gelman and Hill (2007).

Preliminaries

  1. Install WinBUGS if it is not already installed on your system. The program, and all the information for installation can be found at
  2. Install R and the R2WinBUGS package if they are not already installed. R and all available packages can be found at
  3. Write the WinBUGS code of the model in the file model4.txt (see below)
  4. Run the R script to call WinBUGS from R (see below)

WinBUGS code

Copy and paste the following WinBUGS code to create the model4.txt file.

model;

{

for(p in 1:P){ #1#

for(s in 1:S) { #2#

logit(pi[p,s,1]) <- eta.per[p,1] + eta.sit[s,1] #3#

y[p,s,1] ~dbern(pi[p,s,1]) #4#

y.rep[p,s,1]~dbern(pi[p,s,1]) #5#

t[p,s,1]<-pow(y[p,s,1]-pi[p,s,1],2)/(pi[p,s,1]*(1-pi[p,s,1])) #6#

t.rep[p,s,1]<-pow(y.rep[p,s,1]-pi[p,s,1],2)/(pi[p,s,1]*(1-pi[p,s,1])) #7#

logit(pi[p,s,2]) <- eta.per[p,2] + eta.sit[s,2] #8#

y[p,s,2] ~ dbern(pi[p,s,2]) #9#

y.rep[p,s,2] ~ dbern(pi[p,s,2]) #10#

t[p,s,2]<-pow(y[p,s,2]-pi[p,s,2],2)/(pi[p,s,2]*(1-pi[p,s,2])) #11#

t.rep[p,s,2]<-pow(y.rep[p,s,2]-pi[p,s,2],2)/(pi[p,s,2]*(1-pi[p,s,2])) #12#

logit(pi[p,s,3]) <- b[1]*eta.per[p,1]+b[2]*eta.per[p,2]+zp3[p]+b[1]*eta.sit[s,1]+b[2]*eta.sit[s,2]+zs3[s] #13#

y[p,s,3] ~ dbern(pi[p,s,3]) #14#

y.rep[p,s,3] ~ dbern(pi[p,s,3]) #15#

t[p,s,3]<-pow(y[p,s,3]-pi[p,s,3],2)/(pi[p,s,3]*(1-pi[p,s,3])) #16#

t.rep[p,s,3]<-pow(y.rep[p,s,3]-pi[p,s,3],2)/(pi[p,s,3]*(1-pi[p,s,3])) #17#

logit(pi[p,s,4]) <- b[3]*eta.per[p,1]+b[4]*eta.per[p,2]+zp4[p]+b[3]*eta.sit[s,1]+b[4]*eta.sit[s,2]+zs4[s] #18#

y[p,s,4] ~ dbern(pi[p,s,4]) #19#

y.rep[p,s,4] ~ dbern(pi[p,s,4]) #20#

t[p,s,4]<-pow(y[p,s,4]-pi[p,s,4],2)/(pi[p,s,4]*(1-pi[p,s,4])) #21#

t.rep[p,s,4]<-pow(y.rep[p,s,4]-pi[p,s,4],2)/(pi[p,s,4]*(1-pi[p,s,4])) #22#

} #23#

eta.per[p,1:2] ~ dmnorm(mu.eta.per[],R[,]) #24#

zp3[p]~dnorm(0,tt) #25#

zp4[p]~dnorm(0,tt) #26#

} #27#

tt ~ dgamma(0.0001,0.0001) #28#

vt<-1/tt #29#

mu.eta.per[1]<-0 #30#

mu.eta.per[2]<-0 #31#

R[1:2,1:2]~dwish(O[,],2) #32#

O[1,1]<-0.001 #33#

O[2,1]<-0 #34#

O[1,2]<-0 #35#

O[2,2]<-0.001 #36#

for( s in 1 : S) { #37#

eta.sit[s,1:2] ~ dmnorm(mu.eta.sit[],E[,]) #38#

zs3[s] ~ dnorm(tau[3],tb) #39#

zs4[s] ~ dnorm(tau[4],tb) #40#

} #41#

tb ~ dgamma(0.0001,0.0001) #42#

vb<-1/tb #43#

mu.eta.sit[1]<-tau[1] #44#

mu.eta.sit[2]<-tau[2] #45#

E[1:2,1:2]~dwish(G[,],2) #46#

G[1,1]<-0.001 #47#

G[2,1]<-0 #48#

G[1,2]<-0 #49#

G[2,2]<-0.001 #50#

for(r in 1:4){ #51#

tau[r]~ dnorm(0,0.0001) #52#

b[r]~ dnorm(0,0.0001) #53#

} #54#

Psi.per12[1:2,1:2]<-inverse(R[1:2,1:2]) #55#

Psi.sit12[1:2,1:2]<-inverse(E[1:2,1:2]) #56#

Corr.per12<-Psi.per12[1,2]/sqrt(Psi.per12 [1,1]*Psi.per12[2,2]) #57#

Corr.sit12<- Psi.sit12[1,2]/sqrt(Psi.sit12[1,1]*Psi.sit12[2,2]) #58#

tau3<-b[1]*tau[1]+b[2]*tau[2]+tau[3] #59#

tau4<-b[3]*tau[1]+b[4]*tau[2]+tau[4] #60#

Tobs<-sum(t[,,]) #61#

Trep<-sum(t.rep[,,]) #62#

p.value <- step(Trep-Tobs) #63#

}

Brief explanation of the WinBUGS code

From lines 1 to 23 the likelihood of the model in Equations(10) and (11) of the paper is specified. In lines 5, 10, 15, and 20 replicated data are being generated to be used in the posterior predictive check analysis (see section Bayesian model selection and model checking). From lines 24 to 54, theprior distributions for the parameters are specified. From 55 tothe end of the code, covariance matrices, correlations,reparametrized values of τ (see Appendix C, Equation (23))and the p-value are obtained. The step() function is used to calculate the variable “p.value“, which takes the value 1 if Trep-Tobs ≥ 0 and zero otherwise. Hence, the posterior mean of “p.value” is the proportion of iterations for which Trep ≥ Tobs.

R code to call WinBUGS from R via the R2WinBUGS package

library(R2WinBUGS)

dat=read.table("datafile.dat",header=F)

P=nrow(dat)

S=11

y=array(as.matrix(dat),dim=c(P,S,4))

data=list("P","y","S")

inits=function(){

list(eta.per=structure(.Data = c(rnorm(2*P,0,1)),.Dim = c(P, 2)), zp3=c(rnorm(P,0,1)), zp4=c(rnorm(P,0,1)), eta.sit=structure(.Data = c(rnorm(2*S,0,1)),.Dim = c(S, 2)), zs3=c(rnorm(S,0,1)), zs4=c(rnorm(S,0,1)), b=c(rnorm(4,0,1)),tau=c(rnorm(4,0,1)), R=structure(.Data = c(1,0,0,1),.Dim = c(2,2)), E=structure(.Data = c(1,0,0,1),.Dim = c(2,2)),tt=rgamma(1,1),tb=rgamma(1,1))

}

parameters=c("b","tau","Psi.per12","Psi.sit12","Corr.per12","Corr.sit12","vt","vb","tau3","tau4","p.value","Tobs","Trep")

mod4<-bugs (data, inits, parameters, "model4.txt",n.chains=3, n.iter=500,n.thin=1)

The structure of the “datafile.dat” file is a rectangular matrix with 679 persons in the rows and 44 columns. The first 11 columns are the Frustration answers under the 11 situation, the next 11 columns are the Antagonistic action responses under the 11 situations, and so on. In this example only 3 chains with 500 iterations each are run.