Model does not contain any samples, network meta-analysis

[Richard]

Hello,

I'm starting to learn Stan and I thought I would begin by trying to convert WinBUGS code from here (http://www.nicedsu.org.uk/evidence-synthesis-tsd-series(2391675).htm) into Stan code.  I started with the TSD2-1 files and have run into problems with 1c. The WinBUGS code is this (taken from the .odc file I've attached):

# Binomial likelihood, logit link

# Random effects model for multi-arm trials

model{                               # *** PROGRAM STARTS

for(i in 1:ns){                      # LOOP THROUGH STUDIES

    w[i,1] <- 0    # adjustment for multi-arm trials is zero for control arm

    delta[i,1] <- 0             # treatment effect is zero for control arm

    mu[i] ~ dnorm(0,.0001)           # vague priors for all trial baselines

    for (k in 1:na[i]) {             # LOOP THROUGH ARMS

        r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood

        logit(p[i,k]) <- mu[i] + delta[i,k]  # model for linear predictor

        rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

#Deviance contribution

        dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k]))  

            +  (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k])))         }

#  summed residual deviance contribution for this trial

    resdev[i] <- sum(dev[i,1:na[i]])       

    for (k in 2:na[i]) {             # LOOP THROUGH ARMS

# trial-specific LOR distributions

        delta[i,k] ~ dnorm(md[i,k],taud[i,k])

# mean of LOR distributions (with multi-arm trial correction)

        md[i,k] <-  d[t[i,k]] - d[t[i,1]] + sw[i,k]

# precision of LOR distributions (with multi-arm trial correction)

        taud[i,k] <- tau *2*(k-1)/k

# adjustment for multi-arm RCTs

        w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]])

# cumulative adjustment for multi-arm trials

        sw[i,k] <- sum(w[i,1:k-1])/(k-1)

      }

  }   

totresdev <- sum(resdev[])           # Total Residual Deviance

d[1]<-0       # treatment effect is zero for reference treatment

# vague priors for treatment effects

for (k in 2:nt){  d[k] ~ dnorm(0,.0001) }

sd ~ dunif(0,5)     # vague prior for between-trial SD

tau <- pow(sd,-2)   # between-trial precision = (1/between-trial variance)

# Provide estimates of treatment effects T[k] on the natural (probability) scale

# Given a Mean Effect, meanA, for 'standard' treatment A, 

# with precision (1/variance) precA

A ~ dnorm(meanA,precA)

for (k in 1:nt) { logit(T[k]) <- A + d[k]  }

}                                    # *** PROGRAM ENDS                          

My attempted Stan code is this: 

data {

  int<lower=0> ND; // Number of data points

  int<lower=0> NT; // Number of treatments

  int<lower=0> NS; // Number of studies

  int<lower=0> sampleSize[ND]; // Number of patients in study and treatment arm 

  int<lower=0> responders[ND]; // Number of responders in study and treatment arm

  int<lower=1> treatment[ND]; // Treatment used in study and treatment arm

  int<lower=1> study[ND]; // Study for data point

  int<lower=1> treatment_base[ND]; // The first or control treatment per study

} 

parameters {

  vector[NT] d;

  vector[NS] mu;

  real delta[NS,NT];

  real<lower=0> sigmasq_delta; 

  real delta_new;

}

transformed parameters {

  real<lower=0> sigma_delta;

  matrix[NS,NT] w;

  matrix[NS,NT] sw;

//  matrix[NS,NT] md;

//  matrix[NS,NT] vard;

//  matrix[NS,NT] sqrt_vard; 

  sigma_delta <- sqrt(sigmasq_delta);

  for (i in 1:ND) {

    if (treatment[i] == 1) {

 w[study[i],treatment[i]] <- 0.0; // adjustment for multi-arm trials is zero for control arm

    } else {

 // adjustment for multi-arm RCTs

      w[study[i],treatment[i]] <- delta[study[i],treatment[i]] - d[treatment[i]] + d[treatment_base[i]];

 // cumulative adjustment for multi-arm RCTs

      sw[study[i],treatment[i]] <- sum(w[study[i],1:(treatment[i]-1)])/(treatment[i]-1);

    }

//  md[study[i],treatment[i]] <- d[treatment[i]] - d[treatment_base[i]] + sw[study[i],treatment[i]];

//        vard[study[i],treatment[i]] <- (pow(sigmasq_delta, 2) * treatment[i])./(2 * (treatment[i]-1));

//  sqrt_vard[study[i],treatment[i]] <- sqrt(vard[study[i],treatment[i]]);

  }

}  

model {

  for (i in 1:ND) {

    if (treatment[i] == 1) {

   responders[i] ~ binomial_logit(sampleSize[i], mu[study[i]]);

} else {

      responders[i] ~ binomial_logit(sampleSize[i], mu[study[i]] + delta[study[i],treatment[i]]);

}

delta[study[i],treatment[i]] ~ student_t(4,d[treatment[i]],sigma_delta);

//  delta[study[i],treatment[i]] ~ normal(md[study[i],treatment[i]],sqrt_vard[study[i],treatment[i]]);

mu[study[i]] ~ normal(0, sqrt(1E5));

    d[treatment[i]] ~ normal(0, sqrt(1E5)); 

    sigmasq_delta ~ inv_gamma(1E-2, 1E-2); 

  }

}

The problems seems to center around the calculation of sw for the multi-arm adjustment, so I've commented out the md, vard, sqrt_vard calculation until that gets dealt with.  If I comment out the sw line I can get samples from the model.  If I define sw as just sum(w[study[i],1:(treatment[i]-1)]) or just 1/(treatment[i]-1) I get samples from the model, but as-is the model seems to fail. I've attached the data I use in Rda and txt files.

Can anyone see what would cause Stan to output that "Stan model 'TSD2-1cRE_Bi_logit' does not contain samples."? Am I doing something fundamentally wrong?

And, since I'm still learning, is there anything else that anyone would suggest to improve the code?

Thank you,

Richard

[Bob Carpenter]

It helps immensely if you can report the *full message* you get back.
I believe that error arises in RStan when the initialization isn't
legal (too bad it doesn't say that explicitly) or some other error
arose during sampling (like illegal generated quantities).

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <data for stan.Rda><data for stan.txt><TSD2-1cRE_Bi_logit.odc>

[Richard]

Sorry about that...

here are whatever error messages were returned

[[1]]

Stan model 'TSD2-1cRE_Bi_logit' does not contain samples.

[[2]]

Stan model 'TSD2-1cRE_Bi_logit' does not contain samples.

[[3]]

Stan model 'TSD2-1cRE_Bi_logit' does not contain samples.

[[4]]

Stan model 'TSD2-1cRE_Bi_logit' does not contain samples.

Warning messages:

1: class "C++Object" is defined (with package slot ‘Rcpp’) but no metadata object found to revise subclass information---not exported?  Making a copy in package ‘.GlobalEnv’ 

2: In .local(object, ...) :

  some chains had errors; consider specifying chains = 1 to debug

Thanks,

Richard

[Ben Goodrich]

consider specifying chains = 1 to debug

Also, the warning message about the C++Object does not seem to do any harm, but goes away if you upgrade your R.

Ben

[Richard]

Good grief.  I've seen that line dozens of times and thought it said something else.  Sorry.

"  Stan can't start sampling from this initial value."                                                 

[301] "Initialization between (--2, -2) failed after 100 attempts. "                                         

[302] " Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model."

[1] "error occurred during calling the sampler; sampling not done"

Okay, this helps. Any suggestions as to the first thing to try?

Thanks,

Richard

[Bob Carpenter]

The problem here is with inits here.
Jonah is sitting next to me and just ran the model.

The error message confused all of us, and it looks like a bug
in Stan or RStan's reporting of random init intervals.

- Bob

[Richard]

Okay, thanks. I appreciate the help.  I'll see what I can do about figuring out appropriate initial values.

Richard

[Richard]

I'm not sure if this is of interest to anyone, but maybe it can help someone else avoid errors.

The problem was that I was only populating part of the "w" array and then summing over a whole row/column.  I didn't initialize the array to be full of zeroes, though, so it failed all the time.  I would look at values output from the summing over "w" but that only included the samples that succeeded.  After I added a print statement the errors popped up and things were pretty clear.

I initialized the matrix with a "w <- rep_array(0.0, NS, NT);" statement and it's running fine, but slow.

I'm not sure if I should start a new topic, but I'm trying to vectorize the model now.  I'd like to do something like... "y[1:ND] ~ normal(a[study[1:ND]] + b[study[1:ND], treatment[1:ND]], se[1:ND]);" and "w[study[1:ND], treatment[1:ND]] <- delta[study[1:ND],treatment[1:ND] - d[treatment[1:ND]];" where ND is the number of data points and study and treatment are vectors containing the corresponding study and treatment for each data point, but I suspect that won't quite work... an array like b[study[1:ND], treatment[1:ND]] will not behave like a vector b[1:ND], right?  Any idea if I can do something different to vectorize?

Thanks,

Richard

[Bob Carpenter]

> On Apr 27, 2016, at 11:41 AM, Richard <richa...@gmail.com> wrote:
>
> I'm not sure if this is of interest to anyone, but maybe it can help someone else avoid errors.
>
> The problem was that I was only populating part of the "w" array and then summing over a whole row/column. I didn't initialize the array to be full of zeroes, though, so it failed all the time. I would look at values output from the summing over "w" but that only included the samples that succeeded. After I added a print statement the errors popped up and things were pretty clear.

Glad you found it --- that's what the print() statement's there for :-)

> I initialized the matrix with a "w <- rep_array(0.0, NS, NT);" statement and it's running fine, but slow.
>
> I'm not sure if I should start a new topic, but I'm trying to vectorize the model now.

Either way's OK.

> I'd like to do something like... "y[1:ND] ~ normal(a[study[1:ND]] + b[study[1:ND], treatment[1:ND]], se[1:ND]);"

This doesn't make sense as there's only two arguments to normal. And if
"se" is for "standard error", you probably don't mean that. The std err
is a measure of estimator accuracy. We just use sigma or tau.

To answer your question, yes you can do what you want to do:

I'm going to ignore that treatment thing because I don' tknow what it is,
but you can program a model with two varying intercepts as

y ~ normal(a[study] + b[study], sigma);

assuming y, study, and sigma are all the same size N, so that it's equal to

for (n in 1:N)
y[n] ~ normal(a[study[n]] + b[study[n]], sigma[n]);

It'll also work if sigma is a scalar (in which case sigma[n] becomes sigma above).
You don't usually have a varying scale for each observation because you can't
can't be fit with a single observation.

> and "w[study[1:ND], treatment[1:ND]] <- delta[study[1:ND],treatment[1:ND] - d[treatment[1:ND]];" where ND is the number of data points and study and treatment are vectors containing the corresponding study and treatment for each data point, but I suspect that won't quite work... an array like b[study[1:ND], treatment[1:ND]] will not behave like a vector b[1:ND], right? Any idea if I can do something different to vectorize?

See the indexing section of the manual for full details.
If a is a vector and nn is an array of integers, then
a[nn] is a vector of the same size as the array nn. And so
on for other types.

- Bob

[Richard]

Hi Bob,

Thanks, I've gotten that working for my fixed effects model, but the random effects is a little different.  Is there anything equivalent to ...

  for (n in 1:N) 

    y[n] ~ normal(a[study[n]] + b[study[n], value2[n]], sigma[n]); 

I only want n = 1:N values from a and b, but a is a vector and b is a matrix or array.  When I just put them into the model it complains that a vector can't be added to a real[,].

Richard

[Richard Zur]

It seems like something like an array of vectors might work.  I'll take a close look at the manual and try and figure it out myself.

Thanks,

Richard

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/yKOOfY7hGVg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I see. No, there's no good way to vectorize that as a one-liner.

The basic idea is that if you have multiple indexes, Stan's indexing
works as:

b[study, value2][i,j] = b[study[i], value2[j]]

So b[study, value2] is of the same type as b (e.g., real[,] or matrix).

That's the same way that R works, by the way:

> a <- matrix(1:6, 2, 3)

> a
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6

> a[c(1,2), c(2,3)]
[,1] [,2]
[1,] 3 5
[2,] 4 6

Now if it's a square matrix, you could use

diag(b[study, value]),

but it wouldn't be efficient. Just use loops to define the value you want.

real b[N];

...
for (n in 1:N)

b_sv[n] <- b[study[n], value[n]];

... use b_study_value ...

- Bob