How to monitor intermediate variables?

[Thao Tran]

Hi everybody,
I have one question: How to monitor an intermediate variables?
I have some model in Jags and want to translate it into Stan. However, I faced the difficulty to moninor an intermediate variables:
The model in Jags has some parts like:

model {
  # Initial condition: At the lowest age
  S[1] <- S1;  I[1] <- I1;  R[1] <- R1
  rateRecover     <- 1-exp(-sigma)
  # From the 2nd age-group
  for (i in 2:84)
  {
    # lambda[i] = age-dependent-FOI.
    lambda[i] <-  beta1 + 2*beta2*age[i]
    # Newly infecteds, recovers, susceptible again.
    rate[i]   <- 1-exp(-lambda[i])
    Inew[i]   ~ dbin(rate[i], S[i-1])                            # New infected
    Rnew[i]   ~ dbin(rateRecover, I[i-1])                        # New recovers
    # Transition process
    S[i]  <- S[i-1]  - Inew[i] ;    I[i]  <- I[i-1]  + Inew[i] - Rnew[i];    R[i]  <- R[i-1]  + Rnew[i]
    # Probability to be infected
    p[i]  <- (I[i] + R[i])/ (I[i] + R[i] + S[i])
    # Put in a logistic regression
    Pos[i] ~ dbin(p[i], Tot[i])
  }
   # Low information prior for parameters
   beta1 ~ dnorm(0,.0001);   beta2 ~ dnorm(0,.0001);   sigma ~ dnorm(8,0.0001)                                        
}

The important part that I could not do the same thing in Stan is:
  Inew[i]   ~ dbin(rate[i], S[i-1])                            # New infected
  Rnew[i]   ~ dbin(rateRecover, I[i-1])                        # New recovers

I was thinking to do something like this for the model in Stan:

model {
  real rateRecover;
  real lambda[n];
  real rate[n];
  int S[n]; int Inew[n];
  int I[n]; int Rnew[n];
  int R[n]; real p[n];
  real prob[n];
 
  rateRecover <- 1-exp(-sigma);
  S[1] <- S1; Inew[1] <- 0;
  I[1] <- I1; Rnew[1] <- 0;
  R[1] <- R1; p[1]    <- (I1+R1)/(I1+R1+S1);
  lambda[1] <- p1+p2*0.5;
  rate[1]   <- 1-exp(-lambda[1]);
  for (i in 2:n) {
    lambda[i] <- p1+p2*age2[i];
    rate[i]   <- 1 - exp(-lambda[i]);
    Inew[i]   ~ binomial(S[i-1] , rate[i])  ;                  # New infected
    Rnew[i]   ~ binomial(I[i-1] , rateRecover) ;               # New recovers# Transition process
    S[i]  <- S[i-1]  - Inew[i];
    I[i]  <- I[i-1]  + Inew[i] - Rnew[i];
    R[i]  <- R[i-1]  + Rnew[i]           ;                       # Probability to be infected
    p[i]  <- (I[i] + R[i])/ (I[i] + R[i] + S[i]);                # Put in a logistic regression
  }
 
 
  Pos2 ~ binomial(Tot2, p);
  # Low information prior for parameters
  p1 ~ normal(0,1000);
  p2 ~ normal(0,1000);
  sigma ~ normal(8,10)   ;                                      # Prior for recovery rates
  # Duration in infected class
}

However, when I ran this model in Stan, I ran into the following error:

[459] "  Error evaluating the log probability at the initial value."                                                                          
[460] "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"                        
[461] "Exception thrown at line 39: stan::math::binomial_log: Probability parameter is -0.0480671, but must be between (0, 1)"                
[462] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
[463] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."                              
[464] "Rejecting initial value:"                                                                                                              
[465] "  Error evaluating the log probability at the initial value."                                                                          
[466] "Rejecting initial value:"                                                                                                              
[467] "  Log probability evaluates to log(0), i.e. negative infinity."                                                                        
[468] "  Stan can't start sampling from this initial value."                                                                                  
[469] "Rejecting initial value:"                                                                                                              
[470] "  Log probability evaluates to log(0), i.e. negative infinity."                                                                        
[471] "  Stan can't start sampling from this initial value."                                                                                  
[472] "Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:"   

I tried to specify "good" starting values but the sampling could not be initialized.
I also tried move verything else into the transformed parameters block like:

transformed parameters {
  real<lower = 0, upper = 1> rateRecover;
  real<lower = 0> lambda[n];
  real<lower = 0, upper = 1> rate[n];
  real<lower = 0, upper = 1> p[n];
  rateRecover <- 1-exp(-sigma);
  lambda[1] <- p1+p2*0.5;
  rate[1]   <- 1-exp(-lambda[1]);

  {int S[n]; int Inew[n];
  int I[n]; int Rnew[n];
  int R[n];
  S[1] <- S1; Inew[1] <- 0;
  I[1] <- I1; Rnew[1] <- 0;
  R[1] <- R1; p[1]    <- (I1+R1)/(I1+R1+S1);
  for (i in 2:n) {
    lambda[i] <- p1+p2*age2[i];
    rate[i]   <- 1 - exp(-lambda[i]);
    Inew[i]   <- binomial_rng(S[i-1], rate[i]) ;                 # New infected
    Rnew[i]   <-  binomial_rng(I[i-1] , rateRecover);            # New recovers# Transition process
    S[i]  <- S[i-1]  - Inew[i];
    I[i]  <- I[i-1]  + Inew[i] - Rnew[i];
    R[i]  <- R[i-1]  + Rnew[i]           ;                       # Probability to be infected
    p[i]  <- (I[i] + R[i])/ (I[i] + R[i] + S[i]);                # Put in a logistic regression
  }

}
}

But now I got this error:

ERROR: random number generators only allowed in tranformed data block, generated quantities block or user-defined functions with names ending in _rng; found function=binomial_rng in block=transformed parameter

ERROR at line 35

 33:        lambda[i] <- p1+p2*age2[i];
 34:        rate[i]   <- 1 - exp(-lambda[i]);
 35:        Inew[i]   <- binomial_rng(S[i-1], rate[i]) ;                 # New infected
                                                      ^

Does any one know how I can improve what I was doing?
In the attachment is the model I have run.
Thank you a lot.

Thank you.

[Bob Carpenter]

It's impossible to tell how to interpret a JAGs
program without your telling us which variables are
given as data and which are defined by the model.

But I'm also not sure what you mean by monitor. You
mean save in the output for RStan or one of the other
interfaces?

If Inew and Rnew are discrete parameters (not given as
data), then you can't translate this directly into Stan.
You have to marginalize out Inew and Rnew (which might be
computationally intractable if the counts are high). Or if nothing
that depends on Inew or Rnew is used in the rest of the model,
you can move everything that does depend on them into the
generated quantities block and make draws from discrete distributions
like this using binomial_rng.

- Bob

> I was thinking to do something like this for the model in Stan like this:

> --
> 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.
> <BP19_SampleInModel_Nov16.txt><BP19_SampleInTransformedParameter_Nov16.txt>

[Thao Tran]

Thank you Bob.
The data are: Age2, Pos2 and Tot2.
And I think I could not move everything into the generated quantities block because Inew, Rnew are intermediate variables that are involved in the calculation of p (is the probability of being infected: Pos2 ~ binomial(p, Tot2)).
And the idea of marginalization out the discrete paramter is quite complicated in this case???

[Bob Carpenter]

Sorry for not getting back until now. I'm trying to clear my inbox
and found a bunch of old old messages.

This is an unusual way to do a regression a rate; we usually
use logistic links:

> rate[i] <- 1-exp(-lambda[i])

If the counts are low enough, you can marginalize out. The rule
is just

p(a) = SUM_n p(a, n)

where n ranges over all values for b in p(a, b). You just need
to write the density directly in Stan.

- Bob

[Thao Tran]

Thank you Bob for giving the feedback.

I will give it a try.

Have a nice afternoon.

Thao.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
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/Lj51QxQbZUk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Trần Mai Phương Thảo
Master Student - Master of Statistics
Hasselt University - Belgium.
Email: Thao...@gmail.com / maiphuong...@student.uhasselt.be
Phone number: + 84 979 397 410+ 84 979 397 410 / 0032 488 0358430032 488 035843