Auto-regressive ordered logistic regression

[Gregory]

Hello everyone,

I'm putting together an auto-regressive ordered logistic regression model, where the interpretation relates to the value of the latent variable yStar.

At the latent-variable level, the model is the following: yStar_t ~ logistic(beta * x1 + rho * yStar_{t -1}, 1)

At the observed-variable level: y_t ~ ordered_logistic(beta * x1 + rho * yStar_{t -1}, kappa)

I started with the simplest version of this model: a single cut point and a single covariate.

I set the model up in Stan in two steps:

(1) draws are made from the posterior for the values of yStar using a truncated logistic density: given the observed value of y, the latent variable yStar is known to be within a given range as defined by the cut point parameters

(2) ordered_logistic() is used to estimate the regression of interest for the observed variable

The problem is that I get the following error and I seem to have run out of ideas for where I'm going wrong:

[299] "  Log probability evaluates to log(0), i.e. negative infinity."                                       

[300] "  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."

Any suggestions for would great!

Simulation and model code:

library(rstan)

n_obs <- 2000
year <- 1:n_obs
x1 <- rnorm(n_obs, -1, 1) # Covariate

beta1 <- 0.4 # Parameter on covariate x1
rho <- 0.7 # Auto-regressive parameter

kappa <- c(-Inf, 0, Inf) # Cut points
k <- length(kappa) - 2 # Number of cut points

S <- data.frame(yStar = NA, y = NA, year = year, x1 = x1)

# Simulate first-order auto-regressive ordinal logistic regression data
for(t in 1:nrow(S)) {
 
  if(t == 1) { # Generate latent (linear) variable
    S$yStar[t] <- 0 # Arbitrary set yStar for t = 1
  } else {
    S$yStar[t] <- rho * S$yStar[t - 1] +
                  beta1 * x1[t] +
                  rlogis(1, 0, 1)
  }

  for(i in 1:(length(kappa) - 1)) { # Generate observed (ordered) variable
    if(S$yStar[t] > kappa[i] & S$yStar[t] <= kappa[i + 1]) S$y[t] <- i
  }

}

model_data <- list(n_obs = n_obs, y = S$y, year = S$year, x1 = S$x1)

model_code <- "
data {

  int n_obs;

  int y[n_obs];
  int year[n_obs];
  real x1[n_obs];

}

parameters {
 
  vector[n_obs] yStar; # Latent variable
  real rho; # Auto-regression parameter
 
  vector[1] beta; # Standard covariate parameters
  ordered[1] kappa; # Cut points

}

transformed parameters {
}

model {

  // Priors
  rho ~ normal(0, 5);
 
  beta ~ normal(0, 5);
  kappa ~ normal(0, 5);

  // Prior on first latent variable
  if(y[1] == 1) yStar[1] ~ normal(0, 5) T[, kappa[1]];
  if(y[1] == 2) yStar[1] ~ normal(0, 5) T[kappa[1], ];

  // Draws for yStar
  for(t in 2:n_obs) {

    if(y[t] == 1) {
      yStar[t] ~ logistic(rho * yStar[t - 1] +
                          beta[1] * x1[t],
                          1) T[, kappa[1]];
    }
    if(y[t] == 2) {
      yStar[t] ~ logistic(rho * yStar[t - 1] +
                          beta[1] * x1[t],
                          1) T[kappa[1], ];
    }
  }

  // Likelihood
  for(t in 2:n_obs) {
    y[t] ~ ordered_logistic(rho * yStar[t - 1] + beta[1] * x1[t], kappa);
  }
}

generated quantities {
}"


# Compile to C
model_latent <- stan_model(model_code = model_code)


# Sample from posterior
posterior_latent <- sampling(model_latent, data = model_data,
                             warmup = 1000, iter = 2000,
                             open_progress = FALSE, chains = 1)

[Gregory]

From what I can tell, the problem relates to the draws from truncated distributions: that the support of the truncated distributions need to fully overlap with the constraints on the parameters. The problem is then that yStar is unconstrained, but the draws are from a truncated distribution. I see that Bob Carpenter has noted this a number of times in his responses to others.

So, in the end, I think I just need to constrain the parameters yStar to be between the relevant cut points kappa at each iteration.

By example, in a model with a single cut point:

If y_i = 1, constrain yStar_i to be less than kappa,

If y_i = 2, constrain yStar_i to be greater than kappa

The constraint thus depends on the value(s) of kappa drawn at a given iteration. Can I do this in the transformed parameter or model block? Is there another way to do this?

Thank again!

[Ben Goodrich]

Hi Greg,

On Saturday, March 12, 2016 at 1:01:02 PM UTC-5, Gregory wrote:

The constraint thus depends on the value(s) of kappa drawn at a given iteration. Can I do this in the transformed parameter or model block? Is there another way to do this?

 
kappa proposed, not kappa drawn

It can be done but you have to account for the change-of-variables and such models tend not to mix well.

parameters{
  pos_ordered[J - 1] kappa;
  vector<lower=0,upper=1>[N] u;
  real alpha;     // be sure to estimate this intercept since kappa is positive
  vector[K] beta;
}
model {
  vector[N] yStar;
  for (i in 1:N) {
    if (y[i] == 1) {
      yStar[i] <- u[i] * kappa[1];
      increment_log_prob(log(kappa[1]));
    }
    else if (y[i] < J) {
      real diff;
      diff <- kappa[y + 1] - kappa[y];
      yStar[i] <- u[i] * diff + kappa[y];
      increment_log_prob(log(diff));
    }
    else {
      yStar[i] <- kappa[J] / (1 - u[i]);
      increment_log_prob(log(kappa[j]) + 2 * log1m(u[i]));
    }
}

Ben

[Gregory]

Thanks a lot for this Ben.

I'm going to have to look into how transforms work here, so I suspect I'll have a follow-up question in the next few days.

(Edit: found own mistake)