stan-network.R

library("qgraph")
library("igraph")
library("IsingSampler")
library("IsingFit")
library('rstan')

set.seed(1337)
Kappa <- as.matrix(get.adjacency(watts.strogatz.game(1,10,1,0)))
Kappa[1,5] <- Kappa[5,1] <- Kappa[3,9] <- Kappa[9,3] <-Kappa[7,1] <-Kappa[1,7] <- Kappa[5,9] <-Kappa[9,5] <-Kappa[7,3] <-Kappa[3,7] <-.5

true_network <- qgraph(Kappa, layout="circle")

data <- IsingSampler(10000, Kappa, thresholds=c(rep(0,10)),
  beta = 1, nIter = 100, method = "MH")

### Stan model
ndim=ncol(data)
nsub=nrow(data)
II=diag(ncol(data))

sm<-'
data{
int ndim;
int nsub;
matrix[ndim,ndim] II;
int dat[nsub,ndim];
vector[ndim] dat2[nsub];
}

parameters{
vector[ndim*ndim-ndim] Avec;
}

transformed parameters{
matrix[ndim,ndim] A;

for(i in 1:ndim){
for(j in 1:ndim){
if(i==j) A[i,j]=0;
if(j>i) A[i,j]=Avec[(i-1)*ndim+j-i];
if(j<i) A[i,j]=Avec[(i-1)*ndim+j-i+1];
}
}

}

model{
Avec~normal(0,1);

for(rowi in 1:nsub){
dat[rowi] ~ bernoulli_logit(A * dat2[rowi]);
}

}
'

sf=stan(model_code = sm,data = list(dat=data,II=II,ndim=ndim,nsub=nsub,dat2=data),chains=3,cores = 3,iter=200)
sum=summary(sf)
A=diag(ndim)
A[row(A)!=col(A)]=sum$summary[,'50%'][1:(ndim^2-ndim)]

par(mfrow=c(1,3))
qgraph(Kappa, layout="circle")
estimated_network <- IsingFit(data, family='binomial', progressbar=T, gamma=0, layout="circle", plot= FALSE)
qgraph(estimated_network$weiadj, layout=true_network$layout)
qgraph(A, layout="circle")

### for windows stan temp folder bug - run if repeating fit
# tmpdir=tempdir()
# tmpdir=gsub('\\','/',tmpdir,fixed=T)
# system(paste0("rm ",tmpdir,'/*.*'))