Network meta analysis model

Modeling
1

Dear Stan developers team,
I’m trying to implement a network meta analysis model in Stan introducing the data as sparse matrices. The following is a simplified form of the model. Is there a clear reason why the sampler cannot initialize? Any suggestions for fixing the issue or alternative parametrization?
Thanks,
Shirin

this is data in Bugs implementation format

t = cbind(c(1,1,1,3,3,4,4), c(3,2,2,4,4,5,5), c(0,0,4,0,0,0,0))
y = cbind(c(-1.22,-0.7,-0.3,-0.24,-0.73,-2.2,-1.8), c(-1.53,-2.4,-2.6,-0.59,-0.18,-2.5,-2.1),
c(0,0,-1.2,0,0,0,0))
se = cbind(c(0.504,0.282,0.505,0.265,0.335,0.197,0.200), c(0.439,0.258,0.510,0.354,0.442,0.190,0.250),
c(0,0,.478,0,0,0,0))

the data is converted into matrix format

nmax = ncol(t)
ns = nrow(t)
nt = max(t)
y1 = matrix(0, ns, nt)
se1 = matrix(0, ns, nt)
for (i in 1:ns) {
for (j in 1:nmax) {
y1[i,t[i,j]] = y[i,j]
se1[i,t[i,j]] = se[i,j]
}
}
b = t[,1]
nz = length(which(y1!=0))

stan implementation

data {
int ns; //number of studies
int nt; // number of treatments
int nz; // number of nonzero elements of y
matrix[ns,nt] y; // matrix of estimates
matrix[ns,nt] se; // matrix of se’s for the estimates

}
transformed data {
int v[nz]; // vector of column indices for non-zero values in y
v = csr_extract_v(y);
}
parameters {
vector[nt] theta; // effect sizes
}
model {
theta ~ normal(0, 10);
for (i in 1:ns) {
for (k in 1:nz) {
y[i,v[k]] ~ normal(theta[v[k]], se[i,v[k]]); // likelihood for treatment v[k] in study i
}
}
}

2

The sampler won’t be able to initialize if any of the elements of se you use aren’t greater than zero. Make sure and double check that is not the case.

3

Whoops I see you included your data sorry! Giving this a run…

4

Yeah, the error I get is:

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Scale parameter is 0, but must be > 0!

I think it’s how the first index is handled here:

for (i in 1:ns) {
  for (k in 1:nz) {
    y[i,v[k]] ~ normal(theta[v[k]], se[i,v[k]]); // likelihood for treatment v[k] in study i
  }
}

There are only nz nonzeros, so you only want one loop over all the nonzero entries in the matrix. Right now you’re looping over ns * nz things, so that’s probly gonna break.

Since the non-zeros of se match the nonzeros of y, you could just extract them both with csr_extract_w like:

// In transformed data
vector[nz] y_w = csr_extract_w(y);
vector[nz] se_w = csr_extract_w(se);

// In model
for (k in 1:nz)
  y_w[k] ~ normal(theta[v[k]], se_w[k]); // likelihood for treatment v[k] in study i

But also since this is transformed data, you could write your own functions that do the extraction and be verbose about things if that helps you avoid mistakes (performance isn’t really critical in the transformed data blocks – and vanilla Stan code is pretty fast). Working with sparse formats is tricky.

Hope that helps!

5

Yes, thanks a lot!

6

Shirin_Golchi,

did your NMA conduct successfully?

Thanks a lot.