We assume that 50 individuals were sampled and that yearly measurements were taken on these individuals during 5 years. We simulate some data where the response variable is potentially influenced by 4 covariates, which we assume do not change through time because study participants only answered our questionnaire in the beginning of the study. These data are stored in a matrix where the rows correspond to different individuals and the columns to different measurements made on these individuals. To simulate these data, we will need to generate samples from a multivariate normal distribution using the package “mvtnorm”:

rm(list=ls(all=TRUE)) 
library(mvtnorm)
set.seed(18)

#set up covariates and basic settings
nind=50
ntime=5
ncov=4
cov1=matrix(rnorm(nind*ncov),nind,ncov)

#parameters
betas=rnorm(ncov+1); betas[c(2,4)]=0
true.betas=betas
media=betas[1]+betas[2]*cov1[,1]+betas[3]*cov1[,2]+betas[4]*cov1[,3]+betas[5]*cov1[,4]
true.sig2=sig2=0.5
true.rho=rho=0.9

#create AR1 covariance matrix
H <- abs(outer(1:ntime, 1:ntime, "-"))
Sigma <- sig2 * rho^H

#create data
dat=matrix(NA,nind,ntime)
for (i in 1:nind){
  dat[i,]=rmvnorm(1,mean=rep(media[i],ntime),sigma=Sigma)
}
colnames(dat)=paste('measur',1:ntime,sep='')
rownames(dat)=paste('person',1:nind,sep='')

#export covariates, response variables and slope coefficients
setwd('U:\\uf\\courses\\bayesian course\\rmarkdown\\AR1 model')
write.csv(dat,'fake data.csv')
write.csv(cov1,'fake covariates.csv',row.names=F)
write.csv(true.betas,'fake betas.csv',row.names=F)

Next, we specify our JAGS model. A key detail here is that, similar to the univariate normal distribution, the multivariate normal distribution is parameterized in JAGS using the precision matrix (i.e., the inverse of the covariance matrix). The JAGS model provided below is stored in the file “AR1_JAGS.R”:

model{
  #calculate covariance matrix
  for (i in 1:ntime){
    for (j in 1:ntime){
    Sigma[i,j]<- sig2*(rho^H[i,j])
    }
  }
  invSigma[1:ntime,1:ntime]<-inverse(Sigma[1:ntime,1:ntime])
  
  #likelihood
  for (i in 1:nind){
    tmp[i]=inprod(betas,cov1[i,])  
    for (j in 1:ntime){
      media[i,j]<-tmp[i]
    }
    y[i,1:ntime] ~ dmnorm(media[i,],invSigma[,])  
  }
  
  #priors
  rho~dunif(0,1)
  sig2 <- 1/prec
  prec~dgamma(0.1,0.1)
  for (i in 1:(ncov+1)){
    betas[i]~dnorm(0,1/10)
  }
}

Note that we use the inprod() function to tell JAGS to calculate the inner product between two vectors “betas” and “cov[i,]”. If you know a little bit of linear algebra, you will know that this inner product is equal to:

\[\begin{pmatrix} \beta_0 & \beta_1 & ... & \beta_p \end{pmatrix} \times \begin{pmatrix} 1 \\ x_{i1} \\ ... \\ x_{ip} \end{pmatrix} = \beta_0 + \beta_1 x_{i1} + ... + \beta_p x_{ip}\]

We use the function inprod() here to avoid explicitly writing \(\beta_0 + \beta_1 x_{i1} + ... + \beta_p x_{ip}\). You will appreciate the convenience of this if you have many covariates. Next, we fit this model and store its output to compare it to those from a simple linear regression using the lm function:

rm(list=ls(all=TRUE)) 
library(jagsUI)

# MCMC settings 
ni <- 10000 #number of iterations
nt <- 5    #interval to thin 
nb <- ni/2  #number of iterations to discard as burn-in
nc <- 3     #number of chains

#set parameters to monitor
params=c("sig2","rho",'betas')

#run model
setwd('U:\\uf\\courses\\bayesian course\\rmarkdown\\AR1 model')
dat=read.csv('fake data.csv',as.is=T)[,-1]
cov1=read.csv('fake covariates.csv',as.is=T)

nind=nrow(dat)
ntime=ncol(dat)
ncov=ncol(cov1)

#auxiliary matrix
H <- abs(outer(1:ntime, 1:ntime, "-"))

#create function that sets the initial values
inits=function(){
  list(betas=rnorm(ncov+1,0,1), prec=rgamma(1,1,1), rho=runif(1,0,1))
}

#fit the model
model2=jags(model.file="AR1_JAGS.R",
            parameters.to.save=params,inits=inits,
            data=list(cov1=cbind(1,cov1),y=dat,nind=nind,ntime=ntime,H=H,ncov=ncov),
            n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE)

# summary(model2)
betas=model2$sims.list$betas
colnames(betas)=paste('betas',1:ncol(betas),sep='')
write.csv(betas,'jags ar1 results.csv',row.names=F)