##Example code for simulating and fitting a bivariate multilevel autoregressive model, with random means and Phi,
##using the ML Wishart specification. After fitting the model, standardized parameters (BP, WP and grand) are calculated.
##The File "Biv_IW_databased.txt" contains the winbugs modelcode. Put this file in the working directory for this R script before you run this code.


####Load Packages####
###first install these packages if not yet installed###
library(MASS)
library(R2WinBUGS)
library(lme4)
library(DataCombine)

####Set dimensions####

tt= 200 ## number of time points
nt=tt ## number of time points
np= 50 #number of subjects

ny=2    ##number of indicators
ne=2    ##number of latent states
nr=6    ##number of random parameters

####Random effects#####

d=matrix(0,np,ny,byrow=T) ####storage for random means 
H=matrix(0,ne*np,ne,byrow=T)####storage for random Phi matrix

parmu=c(.4,.2,.2,.4,3,2) ####the fixed effects for the Phi matrix, followed by those for the random means
parcov= matrix(c(.01, 0, 0, 0, 0, 0,
                 0,.01,0, 0,0,0,
                 0,0, .01, 0 ,0,0,
                 0, 0, 0,0.01, 0, 0,
                 0, 0, 0, 0, 1, 0,
                 0, 0, 0, 0, 0, .5),6,6,byrow=TRUE) ###the covariance matrix for Phi and the random means

truephis=c(0)

for (j in 1:np) { 
  muphi=mvrnorm(1,mu=parmu,Sigma=parcov)
  
  
  H[(j+j-1),1]= muphi[1] ##autoregression factor 1
  H[(j+j),2]=muphi[4] ##autoregression factor 2
  H[(j+j-1),2]=muphi[2] ##crossregression 1
  H[(j+j),1]=muphi[3] ##crossregression factor 2
  
  d[j,1]= muphi[5]
  d[(j),2]= muphi[6]
} ###Sample the random parameters and store Phi for each person in H and the means for each person in d.


###Fixed effects###

Q=matrix(c(
  .5,0,
  0,1
),ne,ne,byrow=T) ###the innovation covariance matrix.

##########################################################
###calculate the covariance matrices of the process###
Id = matrix(c(1,0,0,0, ##identity matrix
              0,1,0,0,
              0,0,1,0,
              0,0,0,1),4,4)
truesig1=c(rep(0,np))
truesig2=c(rep(0,np))
truevar1=c(rep(0,np))
truevar2=c(rep(0,np))
truecov=c(rep(0,np))
covmat=matrix(0,ne*np,ne,byrow=T)
cormat=matrix(0,ne*np,ne,byrow=T)

for (j in 1:np) { 
  cvvec= solve(Id-(H[(j+j-1):(j+j),1:ne]%x%H[(j+j-1):(j+j),1:ne])) %*% c(Q) ###calculate the covmatrix
  truevar1[j]=cvvec[1] 
  truecov[j] = cvvec[2]
  truevar2[j]=cvvec[4] 
  truesig1[j]=sqrt(cvvec[1])
  truesig2[j]=sqrt(cvvec[4])
  cvmat= matrix(cvvec,2,2)
  crmat = cov2cor(cvmat)
  covmat[(j+j-1):(j+j),1:ne]= cvmat ##covariance matrices for each person
  cormat[(j+j-1):(j+j),1:ne]= crmat ##correlation matrices for each person
  
  
}


######Simulate data######
#########################

zdata=matrix(0,tt*np,ne) ## storage for z
ydata=matrix(0,tt*np,ny) #### storage for y
ppnum=matrix(0,tt*np,1) ##storage for participant numbers


## simulate data for t=1 ##
for (j in 1:np) { 
  ydata[(j-1)*tt+1,1:ny]=mvrnorm(1,d[j,1:ne],Sigma=covmat[(j+j-1):(j+j),1:ne]) ##sample from normal distribution based on personal mean and covariance matrix
  zdata[(j-1)*tt+1,1:ne]=ydata[(j-1)*tt+1,1:ny]-d[j,1:ne]
  ppnum[(j-1)*tt+1,1]=j
}

## simulate data for the remaining time points ##
for (j in 1:np) { 
  for (i in ((j-1)*tt+2):(j*tt)) {
    inno=mvrnorm(1,mu=rep(0,ne),Sigma=Q)  ###sample innovations 
    
    zdata[i,1:ne]=H[(j+j-1):(j+j),1:ne]%*%zdata[i-1,1:ne]+inno  ###latent autoregressive process
    ydata[i,1:ny]=zdata[i,1:ne]+d[j,1:ne]      ###observed data with mean d
    
    ppnum[i]=j ###participant number for each observation
    
  }
}


###Estimate univariate models with ML using lme4##


#####Make dataset for ML with person-centered, lag1, predictors#####
###See Hamaker, E. L., & Grasman, R. P. (2014). To center or not to center? Investigating inertia with a multilevel autoregressive model. Frontiers in psychology, 5. ###

means=matrix(NA,np,2)
for (j in 1:np)     {
  means[j,]=colMeans(ydata[((j-1)*tt+1):(j*tt),],na.rm=TRUE)
}  

datacent=ydata
for (j in 1:np) { 
  datacent[((j-1)*tt+1):(j*tt),1] = ydata[((j-1)*tt+1):(j*tt),1] - means[j,1]
  datacent[((j-1)*tt+1):(j*tt),2] = ydata[((j-1)*tt+1):(j*tt),2] - means[j,2]
}


lagdata=as.data.frame(cbind(ppnum,datacent)) ###prepare dataframe for ML analyses
colnames(lagdata)=c("ppn","Y1","Y2")

new=slide(lagdata, Var = "Y1", GroupVar = "ppn", NewVar="Y1lag", slideBy = -1)
lagdata=slide(new, Var = "Y2", GroupVar = "ppn", NewVar="Y2lag", slideBy = -1) ###make the lagged, centered predictor variables
lagdata[,c("Y1", "Y2")]=ydata  #uncentered data for the outcome variables 
mldata=lagdata

##fit data ML in lme4###

out1=lmer(Y1 ~  Y1lag + Y2lag + (1+Y1lag + Y2lag|ppn), data=mldata) ###estimate first univariate model in lme4
out2=lmer(Y2 ~ Y1lag + Y2lag + (1+Y1lag + Y2lag|ppn), data=mldata) ###estimate second univariate model in lme4

####Calculate ML prior input####
###Read out the estimated variances from the lme4 output and put it into matrix W, for the prior specification of the 
###covariance matrix of the random parameters###
###see Schuurman, Grasman and Hamaker, in press###

wmat1=VarCorr(out1)
wmat2=VarCorr(out2)

wvars1= c(wmat1$ppn[1],wmat1$ppn[5],wmat1$ppn[9])
wvars2= c(wmat2$ppn[1],wmat2$ppn[5],wmat2$ppn[9])



df=6 ###degrees of freedom for the Wishart in WinBUGS 
W=matrix(0,6,6) #W is the scale matrix for the Wishart in WinBUGS
W[1,1]= df*wvars1[1]
W[2,2]= df*wvars2[1]
W[3,3]= df*wvars1[3] 
W[4,4]= df*wvars2[2] 
W[5,5]= df*wvars1[2]
W[6,6]= df*wvars2[3] 



#######Call Winbugs and fit the full model#####
f=ydata
nr=6 ##number of random parameters
ns=20000 ## number of samples (excluding burnin)
burnin = 10000 ##number of burnin samples
nit = ns + burnin ###number of iterations

dat=list("np","f","tt","W")###data for winbugs

###specify some starting values, for example: ####
stbpre1 = matrix(c(.01, 0, 0, 0, 0, 0,
                   0,.01,0, 0,0,0,
                   0,0, .01, 0 ,0,0,
                   0, 0, 0,0.01, 0, 0,
                   0, 0, 0, 0, .01, 0,
                   0, 0, 0, 0, 0, .01),6,6,byrow=TRUE)
stbpre2 = matrix(c(.015, 0, 0, 0, 0, 0,
                   0,.015,0, 0,0,0,
                   0,0, .015, 0 ,0,0,
                   0, 0, 0,0.015, 0, 0,
                   0, 0, 0, 0, .015, 0,
                   0, 0, 0, 0, 0, .015),6,6,byrow=TRUE)
stbpre3 = matrix(c(.005, 0, 0, 0, 0, 0,
                   0,.005,0, 0,0,0,
                   0,0, .005, 0 ,0,0,
                   0, 0, 0,0.005, 0, 0,
                   0, 0, 0, 0, .005, 0,
                   0, 0, 0, 0, 0, .005),6,6,byrow=TRUE)

phistart = matrix(0,np,nr)
phistart2 = phistart + .1
phistart3 = phistart + .2

inits1 <- list( bmu=c(0,0,0,0,1,1),  b=phistart, bpre=stbpre1) ###for chain1
inits2 <- list( bmu=c(0.1,0.1,0.1,0.1,2,2),  b=phistart2,  bpre=stbpre2) ###for chain 2
inits3 <- list( bmu=c(0.2,0.2,0.2,0.2,1.5,1.5),  b=phistart3,  bpre=stbpre3) ###for chain 3

inits <- list(inits1)#,inits2,inits3) ###inits that are provided to winbugs, list all three inits if running 3 chains

params= c("Qvar","Qpre","b", "bmu","bpre", "bcov", "bcor") ##parameters for which we want to see the estimates

modelfile= "Biv_IW_databased.txt" ###the model file for winbugs. Please save this file in your R working directory.


####Calling R2Winbugs, for more info, run: ?bugs ####

results= bugs(data=dat,
              inits=inits,
              parameters.to.save=params,
              model.file=modelfile,
              n.chains=1,
              n.iter=nit,
              n.burnin=burnin,
              n.thin=1,
              debug=FALSE,
              DIC=TRUE,
              digits=5,
              codaPkg=TRUE,
              bugs.directory="C:/Program Files/WinBUGS14",
              program=c("WinBUGS"),
              bugs.seed=NULL,
              summary.only=TRUE,
              save.history=FALSE,
              over.relax = FALSE)


###tip for potential warning "cannot create file ..... permission denied...": Run R as an administrator.


####Make matrix Phi and vector Sigma for calculating the estimated covariance matrix of the data###
###################################################################################################
codaobject<-read.bugs(results) ###load the MCMC samples
mcmcobject = as.mcmc(codaobject) 
mcmcmatrix = as.matrix(mcmcobject)


H<-array(NA,dim=c(2,2,ns,np)) ###matrix for storing Phi
Q<-array(NA,dim=c(4,ns))  ##vector for storing Sigma

for (pp in 1:np){
  
  for (sample in 1:ns) {
    H[1,1,sample,pp] = mcmcmatrix[sample,paste("b[",pp,",5]", sep = "")]
    H[1,2,sample,pp] = mcmcmatrix[sample,paste("b[",pp,",4]", sep = "")]
    H[2,1,sample,pp] = mcmcmatrix[sample,paste("b[",pp,",3]", sep = "")]
    H[2,2,sample,pp] = mcmcmatrix[sample,paste("b[",pp,",6]", sep = "")]
    
    
    cat("subject", pp,"sample", sample, " \n")
    flush.console()
    
  }  ###read out the relevant coefficients and put into matrix.
  
}
save(H, file="Coda_H.RData")

for (sample in 1:ns) {
  Q[1,sample] = mcmcmatrix[sample,"Qvar[1,1]"]
  Q[2,sample] = mcmcmatrix[sample,"Qvar[1,2]"]
  Q[3,sample] = mcmcmatrix[sample,"Qvar[2,1]"]
  Q[4,sample] = mcmcmatrix[sample,"Qvar[2,2]"]
  
  
  cat("iteration", sample, " \n")
  flush.console()
}

save(Q, file="Coda_Q.RData") ###read out the relevant coefficients and put into vector.


#######################################################################################





######################
#######functions######
Id = matrix(c(1,0,0,0, ##identity matrix
              0,1,0,0,
              0,0,1,0,
              0,0,0,1),4,4)

VAR <- function(ID, Phi,Res){ ##function to obtain person-specific covariance matrices for each iteration of the mcmc sampler
  VECSIGMA = solve(ID-(Phi%x%Phi)) %*% Res
  SIGMA = matrix(VECSIGMA,2,2)
  return(SIGMA)
}

GETS1S2 <- function(ID, Phi,Res){ ##function to obtain person-specific variances for each iteration of the mcmc sampler
  VECSIGMA = solve(ID-(Phi%x%Phi)) %*% Res
  S1S2 =c(VECSIGMA[1],VECSIGMA[4])
  return(S1S2)
}

STAN <- function(b, Sy,Sx){ ##function to obtain person-specific standardized coefficients for each iteration of the mcmc sampler
  beta = b*Sx/Sy
  names(beta)="standardized coefficient"
  return(beta)
}
####################################################################
##########calculate within person variances and standard deviations#####
SIGMAS <-array(0,dim=c(2,2,ns,np)) ###storage array for all the covariance matrices
for(pp in 1:np){ ##calculate person-specific covariance matrices
  for (sample in 1:ns){
    SIGMAS[,,sample,pp] = VAR(Id,H[,,sample,pp],Q[,sample])
  }
}

vars <-array(0,dim=c(ns,2,np)) ##storage array for the variances
colnames(vars)=c("S1","S2") 
for(pp in 1:np){ ##calculate person-specific variances
  for (sample in 1:ns){
    vars[sample,1:2,pp] = GETS1S2(Id,H[,,sample,pp],Q[,sample])
  }
}

sigs <- sqrt(vars) ##calculate person-specific standard deviations

#############################################################
#####calculate within-person standardized coefficients#######
Sphi <-array(0,dim=c(ns,2,np)) ##storage for the within-person standardized coefficients
colnames(Sphi)=c("sphi12","sphi21")
for(pp in 1:np){##calculate the within-person standardized coefficients
  for (sample in 1:ns){
    Sphi[sample,1,pp] = STAN(b=H[1,2,sample,pp],Sy=sigs[sample,1,pp],Sx=sigs[sample,2,pp])
    Sphi[sample,2,pp] = STAN(b=H[2,1,sample,pp],Sy=sigs[sample,2,pp],Sx=sigs[sample,1,pp])
    
  }
}

######Summary results for within-person standardized coefficients########

###person-specific coefficients###
StanResultsSphi = array(NA,dim=c(np,5,2))
colnames(StanResultsSphi)=c("mean","standard dev", "lb 95CI", "median","ub 95CI")
StanResultsSphi[,1,1:2] = t(colMeans(Sphi[,,],na.rm=TRUE)) ###means
for(pp in 1:np){
  StanResultsSphi[pp,2,] = apply(Sphi[,,pp],2,sd, na.rm=TRUE)  ##posterior standard deviations
  sub1 = Sphi[which(is.na(Sphi[,1,pp])==FALSE),1,pp]  ##exclude any na's that are the result of MCMC samples that do not comply with the stationarity restrictions (if these restrictions do not apply then the equation for the person-specific variances is not applicable). Be sure to check if there are very many violations - if so, you may have chosen simulation conditions that are at odds with the stationarity assumption.
  sub2 = Sphi[which(is.na(Sphi[,2,pp])==FALSE),2,pp] ##exclude any na's that are the result of MCMC samples that do not comply with the stationarity restrictions (if these restrictions do not apply then the equation for the person-specific variances is not applicable). Be sure to check if there are many violations.
  xx1 = summary(as.mcmc(sub1),quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975)) ###get the quantiles for the credible intervals
  StanResultsSphi[pp,3,1] = xx1$quantiles[[1]]
  StanResultsSphi[pp,4,1] = xx1$quantiles[[3]]
  StanResultsSphi[pp,5,1] = xx1$quantiles[[5]]
  xx2 = summary(as.mcmc(sub2), quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975))
  StanResultsSphi[pp,3,2] = xx2$quantiles[[1]]
  StanResultsSphi[pp,4,2] = xx2$quantiles[[3]]
  StanResultsSphi[pp,5,2] = xx2$quantiles[[5]]
}

###within-person standardized fixed effects###
Sbmu1=0
Sbmu2=0
for (sample in 1:ns){
  Sbmu1[sample]=sum(Sphi[sample,1,])/np  ###calculate average person-specific standardized coefficient in each sample
  Sbmu2[sample]=sum(Sphi[sample,2,])/np
}
FI_Sbmu1=summary(as.mcmc(Sbmu1[which(is.na(Sbmu1)==FALSE)])) ###obtain summary results
FI_Sbmu2=summary(as.mcmc(Sbmu2[which(is.na(Sbmu2)==FALSE)]))

#########Calculate grand standardized coefficients###########
#############################################################

varmeans=codaobject[,c("bcov[1,1]","bcov[2,2]")]  ###The samples for the variances of the person-specific means, AKA the between-person variances
av_wpmean=apply(vars,1:2,mean)  ### the average within person variances in each iteration of the MCMC sampler.
grandvars=av_wpmean+varmeans[[1]]  ### the grand variances for each iteration of the MCMC sampler
grandsigs=sqrt(grandvars)  ##the grand standard deviations

GSphi <-array(0,dim=c(ns,2,np)) ##storage array for the grand standardized person-specific coefficients
colnames(GSphi)=c("Gsphi12","Gsphi21")


for(pp in 1:np){
  for (sample in 1:ns){ ###calculate grand standardized coefficients
    GSphi[sample,1,pp] = STAN(b=H[1,2,sample,pp],Sy=grandsigs[sample,1],Sx=grandsigs[sample,2])
    GSphi[sample,2,pp] = STAN(b=H[2,1,sample,pp],Sy=grandsigs[sample,2],Sx=grandsigs[sample,1])
    
  }
}

######Summary results for grand standardized coefficients########

###person-specific coefficients###
StanResultsGSphi = array(NA,dim=c(np,5,2)) ### storage array for summary results
colnames(StanResultsGSphi)=c("mean","standard dev", "lb 95CI", "median","ub 95CI")
StanResultsGSphi[,1,1:2] = t(colMeans(GSphi[,,],na.rm=TRUE)) ##means
for(pp in 1:np){
  StanResultsGSphi[pp,2,] = apply(GSphi[,,pp],2,sd, na.rm=TRUE) ##posterior standard deviations
  sub1 = GSphi[which(is.na(GSphi[,1,pp])==FALSE),1,pp]
  sub2 = GSphi[which(is.na(GSphi[,2,pp])==FALSE),2,pp]
  xx1 = summary(as.mcmc(sub1),quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975)) #### quantiles for credible intervals
  StanResultsGSphi[pp,3,1] = xx1$quantiles[[1]]
  StanResultsGSphi[pp,4,1] = xx1$quantiles[[3]]
  StanResultsGSphi[pp,5,1] = xx1$quantiles[[5]]
  xx2 = summary(as.mcmc(sub2), quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975))
  StanResultsGSphi[pp,3,2] = xx2$quantiles[[1]]
  StanResultsGSphi[pp,4,2] = xx2$quantiles[[3]]
  StanResultsGSphi[pp,5,2] = xx2$quantiles[[5]]
}

###grand stardized fixed effects ###
GSbmu1=NA
GSbmu2=NA
for (sample in 1:ns){
  GSbmu1[sample]=sum(GSphi[sample,1,])/np  ###calculate average person-specific grand standardized coefficients
  GSbmu2[sample]=sum(GSphi[sample,2,])/np
}

FI_GSbmu1=summary(as.mcmc(GSbmu1[which(is.na(GSbmu1)==FALSE)]))
FI_GSbmu2=summary(as.mcmc(GSbmu2[which(is.na(GSbmu2)==FALSE)]))

####Calculate between person standardized coefficients####
varmeans=codaobject[,c("bcov[1,1]","bcov[2,2]")]  ###the Between-person variances (the variances of the person-specific means)
BPVars=varmeans[[1]]
BPsigs=sqrt(BPVars)  ###obtain the bp standard deviations
BSphi <-array(0,dim=c(ns,2,np))
colnames(BSphi)=c("Bsphi12","Bsphi21")

for(pp in 1:np){
  for (sample in 1:ns){ ###calculate the between person standardized person-specific coefficients
    BSphi[sample,1,pp] = STAN(b=H[1,2,sample,pp],Sy=BPsigs[sample,1],Sx=BPsigs[sample,2])
    BSphi[sample,2,pp] = STAN(b=H[2,1,sample,pp],Sy=BPsigs[sample,2],Sx=BPsigs[sample,1])
    
  }
}

####summary results for the BP standardized coefficients
StanResultsBSphi = array(NA,dim=c(np,5,2))
colnames(StanResultsBSphi)=c("mean","standard dev", "lb 95CI", "median","ub 95CI")
StanResultsBSphi[,1,1:2] = t(colMeans(BSphi[,,],na.rm=TRUE)) # mean
for(pp in 1:np){
  StanResultsBSphi[pp,2,] = apply(BSphi[,,pp],2,sd, na.rm=TRUE) #posterior standard deviation
  sub1 = BSphi[which(is.na(BSphi[,1,pp])==FALSE),1,pp]
  sub2 = BSphi[which(is.na(BSphi[,2,pp])==FALSE),2,pp]
  xx1 = summary(as.mcmc(sub1),quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975)) ###credible interval
  StanResultsBSphi[pp,3,1] = xx1$quantiles[[1]]
  StanResultsBSphi[pp,4,1] = xx1$quantiles[[3]]
  StanResultsBSphi[pp,5,1] = xx1$quantiles[[5]]
  xx2 = summary(as.mcmc(sub2), quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975))
  StanResultsBSphi[pp,3,2] = xx2$quantiles[[1]]
  StanResultsBSphi[pp,4,2] = xx2$quantiles[[3]]
  StanResultsBSphi[pp,5,2] = xx2$quantiles[[5]]
}

###results for BP fixed effects###
BSbmu1=0
BSbmu2=0
for (sample in 1:ns){
  BSbmu1[sample]=sum(BSphi[sample,1,])/np
  BSbmu2[sample]=sum(BSphi[sample,2,])/np
}
FI_BSbmu1=summary(as.mcmc(BSbmu1[which(is.na(BSbmu1)==FALSE)]))
FI_BSbmu2=summary(as.mcmc(BSbmu2[which(is.na(BSbmu2)==FALSE)]))


#########All results###########
##unstandardized##
summary(codaobject)

##standardized##
###person-specific###
StanResultsSphi ##WPS
StanResultsGSphi ##GS
StanResultsBSphi ##BPS

###fixed effects###
FI_Sbmu1;FI_Sbmu2  ##WPS
FI_GSbmu1;FI_GSbmu2  ##GS
FI_BSbmu1;FI_BSbmu2 ##BPS

###plots###
plot(StanResultsSphi[,4,1],StanResultsSphi[,4,2],pch=16,col="blue",ylim=c(-.1,.5),xlim=c(-.1,.5))
points(StanResultsGSphi[,4,1],StanResultsGSphi[,4,2],pch=15, col="red")
points(StanResultsBSphi[,4,1],StanResultsBSphi[,4,2],pch=17, col="green")


points(FI_Sbmu1$quantiles[[3]],FI_Sbmu2$quantiles[[3]],pch=10,cex=3,lwd=2)
points(FI_GSbmu1$quantiles[[3]],FI_GSbmu2$quantiles[[3]],pch=12,cex=3,lwd=2)
points(FI_BSbmu1$quantiles[[3]],FI_BSbmu2$quantiles[[3]],pch=6,cex=3,lwd=2)
l=seq(from = -1, to = 1, by=0.01)
k=seq(from = -1, to = 1, by=0.01)
lines(k,l)