Bayesian Estimation of Stereotype Logit Model

[jun...@gmail.com]

Dear Stan (rstan) user,

I am trying to estimate a stereotype logit model using rstan, but run into problems. I would appreciate it if anyone could point me to the right direction. I would be more than happy to share the codes once I work it out. Following are two documents that discuss this model. Basically it is a reduced rank multinomial logit model, with the slopes varying across equations by a set of scalar factors that are constrained between 0 and 1 (corner constraints) to retain ordinality and obtain identification.

http://psfaculty.ucdavis.edu/bsjjones/stereotypelogit.pdf

http://deepblue.lib.umich.edu/bitstream/handle/2027.42/64310/3693_ftp.pdf?sequence=1
(see page 3149 for how some of the key parameters are given priors using Dirichlet distribution)

Jun Xu, PhD
Associate Professor
Department of Sociology
Ball State University

Here are some excerpts of the error/warning messages that I got from running the following codes:
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Invalid value of phi: validate transformed params(N4stan5agrad3varE): phi[4] is -0.211538:0, but must be greater than or equal to 0
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Rejecting proposed initial value with zero density.

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Invalid value of phi: validate transformed params(N4stan5agrad3varE)error occurred during calling the sampler; sampling not done
> print(stanFit)
Stan model 'modelString' does not contain samples.

Below are my r codes (see the attached files for the r script and data)
Codes
###

rm(list=ls(all=TRUE))
require(rstan)
require(MASS)
require(VGAM)
require(foreign)


mydta <- read.dta("gssCum7212Teach.dta", convert.factor=F)

usevar <- c("health", "age", "educ", "impinc", "female", "year")
mydta <- mydta[complete.cases(mydta[usevar]),]
              
usedta <- subset(mydta, select=c(health, age, educ, impinc, female))

# set seed to get the same random sample each time
set.seed(47306)
usedta <- usedta[sample(1:nrow(usedta), size = 1000, replace=T),]

# here is the DATA used for estimation
save(usedta, file = "stanData.RData")

myologit <- polr(as.ordered(usedta$health) ~ usedta$age + usedta$educ + usedta$impinc +usedta$female,
                 method = c("logistic"))
summary(myologit)

# stereotype logit model
stereotypeLogitR1 <- rrvglm(as.factor(health) ~ age + educ + impinc + female,
                            multinomial,
                            usedta)
summary(stereotypeLogitR1)
# slotNames(stereotypeLogitR1)

#------------------------------------------------------------------------------
# THE MODEL.

modelString = "
data {
    int<lower=2> M; // number of response levels
    int<lower=1> N; // number of observations
    int<lower=1> K; // number of predictors
    int<lower=1, upper=M> y[N]; // y variables
    row_vector[K] x[N]; // x vectors
    vector<lower=0>[M-1] alpha;
}

parameters {
    vector[K] b;
    vector[M-1] theta;
    vector[M-1] gamma;
}

transformed parameters{
    vector<lower=0, upper=1>[M] phi;
    phi[1] <- 0;
    phi[2] <- gamma[1];
    phi[3] <- gamma[1] + gamma[2];
    phi[4] <- gamma[1] + gamma[2] + gamma[3];
}

model {
    vector[M] prob;
    vector[M] emu;
   
    // specify priors
    for (m in 1:M-1)
        gamma ~ dirichlet(alpha);
    for (k in 1:K)
        b[k] ~ normal(0, 100);
    for (m  in 1:M-1)
        theta[m] ~ normal(0, 100);
    for (n in 1:N) {
            emu[1] <- exp(phi[1]*x[n]*b);
            for (m in 2:M)
                emu[m] <- exp(theta[m] + phi[m]*x[n]*b);
            for (m in 1:M)
                prob[m] <- emu[m]/sum(emu);
            y[n] ~ categorical(prob);
    }
}

# generated quantitaties {
#
# }

" # close quote for modelstring

# translate model to C++ and compile to DSO:
stanDSO <- stan_model(model_code=modelString)

# prepare data
dataMat = as.matrix(usedta)

# Rename for use by generic processing later:
nData = NROW(dataMat)
x = dataMat[,-1]
predictorNames = colnames(dataMat)[-1]
nPredictors = NCOL(x)
y = as.matrix(dataMat[,1])
predictedName = colnames(dataMat)[1]
nYlevels = max(y)
distAlpha = rep(1,3)

dataList = list(
    x = x ,
    y = as.vector( y ) , # BUGS does not treat 1-column mat as vector
    K = nPredictors ,
    N = nData ,
    M = nYlevels,
    alpha = distAlpha
)

# get initial values for b
bInit = - coefficients(stereotypeLogitR1)[-c(1:3)]
thetaInit = - coefficients(stereotypeLogitR1)[c(1:3)]

initsList =  list(
    b = bInit,
    theta = thetaInit)

# specify parameter values for MCMC sampling
parameters = c("b", "theta", "gamma")  # The parameter(s) to be monitored.
adaptSteps = 1000              # Number of steps to "tune" the samplers.
burnInSteps = 1000            # Number of steps to "burn-in" the samplers.
nChains = 3 
seedNum = 60115
# Number of chains to run.
initsChains <- list()
# initials for three chains
for (i in 1:nChains) {
    initsChains[[i]] <- initsList
}
#}
numSavedSteps=10000           # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.

# Sampling MCMC chain
cat( "Sampling final MCMC chain...\n" )

# monitor time
time.used <- proc.time()

stanFit <- sampling(object=stanDSO,
                data=dataList,
                seed = seedNum,
                pars=parameters,
                chains=nChains,
                thin=thinSteps,
                  # init=initsChains
                iter=nPerChain,
                warmup=burnInSteps )
print(stanFit)
###

[Bob Carpenter]

May I also ask you what you read of the Stan doc before
trying this? We need to somehow get the message about
satisfying constraints across, because most of the questions
we're getting lately are from people who didn't realize
that all parameters must be constrained so that any values
satisfying the constraints has a positive (non-zero) density (or
finite log density).

I don't have time to read about a new model, but the
issue you're getting from Stan is exactly what it says it
is. You assigned phi to a value that does not satisfy
its declared constraints. The source of this problem is here:

> parameters {
> vector[K] b;
> vector[M-1] theta;
> vector[M-1] gamma;
> }
>
> transformed parameters{
> vector<lower=0, upper=1>[M] phi;
> phi[1] <- 0;
> phi[2] <- gamma[1];
> phi[3] <- gamma[1] + gamma[2];
> phi[4] <- gamma[1] + gamma[2] + gamma[3];
> }

Because gamma is unconstrained, it's valid for it to have values like -1
or 127, neither of which are going to allow phi to satisfy
that lower=0, upper = 1 constraint. Maybe you want to inv-logit transform
those gammas as is typical in ordered logit models?

Alternative, you can define

real<lower=0, upper=1> gamma1;
real<lower=0, upper= 1 - gamma1> gamma2;
real<lower=0, upper = 1 - gamma2> gamma3;

We have built-in operations to do a lot of the things you're
implementing from scratch, like softmax.

- Bob

> On Oct 27, 2015, at 9:44 AM, jun...@gmail.com wrote:
>
> Dear Stan (rstan) user,
>
> I am trying to estimate a stereotype logit model using rstan, but run into problems. I would appreciate it if anyone could point me to the right direction. I would be more than happy to share the codes once I work it out. Following are two documents that discuss this model. Basically it is a reduced rank multinomial logit model, with the slopes varying across equations by a set of scalar factors that are constrained between 0 and 1 (corner constraints) to retain ordinality and obtain identification.
>
> http://psfaculty.ucdavis.edu/bsjjones/stereotypelogit.pdf
>
> http://deepblue.lib.umich.edu/bitstream/handle/2027.42/64310/3693_ftp.pdf?sequence=1
> (see page 3149 for how some of the key parameters are given priors using Dirichlet distribution)
>

> Here are some excerpts of the error/warning messages that I got from running the following codes:
> Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

> --
> 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.
> <bayesOrdRegStereotypeStan.r><stanData.RData>

[Ben Goodrich]

On Tuesday, October 27, 2015 at 1:00:18 PM UTC-4, Bob Carpenter wrote:

Alternative, you can define

  real<lower=0, upper=1> gamma1;
  real<lower=0, upper= 1 - gamma1> gamma2;
  real<lower=0, upper = 1 - gamma2> gamma3;

Another parameterization would be to make gamma a simplex vector and make phi be the cumulative sum of gamma that is scaled by a free parameter that is constrained to the (0,1) interval.

Ben

[Bob Carpenter]

Yeah, making gamma a simplex is a much better idea, though
that won't give you the gamma0 = 0 thing.

- Bob

[Ben Goodrich]

On Tuesday, October 27, 2015 at 1:41:11 PM UTC-4, Bob Carpenter wrote:

Yeah, making gamma a simplex is a much better idea, though
that won't give you the gamma0 = 0 thing.

Right. Just need to make the simplex vector one element shorter and take gamma0 = 0 to be fixed.

Ben

[jun...@gmail.com]

Bob,

Thanks a lot for your suggestion. In my model block of codes, I specify that

    for (m in 1:M-1)
        gamma ~ dirichlet(alpha);

So here I want gamma to follow a dirichlet distribution so that all phi's are constrained between 0 and 1 and of course, phi[1] and phi[4] are fixed as 0 and one.  Yes, you are right that I need to put constraints on both gamma and alpha (?). I will revise my codes, and update what I got.

I have to admit that I didn't spend a lot of time going over the manual. I just quickly went over the regression section on ordered regression and multi-category logit models. But that's a mistake that I should not have made :(

I haven't quite understood/know how to implement your suggestion here. But I will try and see what happens. Thank you very much.

[jun...@gmail.com]

Now, I declare gamma parameters to be "vector<lower=0, upper=1>[M-1] gamma;"

but I am still getting these out of boundary sampling values. I am a bit confused as to why this is so since I already declare gamma to follow the dirichlet distribution in the model block and I supply value to alpha in the data block and list....

Jun

[Bob Carpenter]

We tried to emphasize the way that sampling statements
work. Anyway, not a big deal as we're planning to drop
the "~" notation anyway because it's just confusing our
users who assume on seeing it that it either (a) works
like BUGS to define a graphical model, or (b) actually
performs a draw as in a forward simulation. It's just
shorthand for incrementing the log density.

The key point is that every set of values that satisfies
the constraints must be in support. If you do

vector<lower=0,upper=1>[K] gamma;

then gamma[1] = gamma[2] = ... = gamma[K] = 0.99 satisfies
the variable's constraint, but is not in support for a Dirichlet.

But then you don't want to have that 0 you had distributed
by dirichlet(alpha) anyway --- it's out of support if you look
at the definition.

If you really just need a simplex, do this:

data {
vector<lower=0>[K] alpha;
...
parameters {
simplex[K] gamma;
...
model {
gamma ~ dirichlet(alpha);
...

- Bob

[jun...@gmail.com]

Bob,

Thank you so much for all your help. You are credited with the successful completion of my attempt. As promised, I am sharing my Stan codes below:

data {
    int<lower=2> M; // number of response levels
    int<lower=1> N; // number of observations
    int<lower=1> K; // number of predictors
    int<lower=1, upper=M> y[N]; // y variables
    row_vector[K] x[N]; // x vectors
    vector<lower=0>[M-1] alpha;
}

parameters {

    vector[K] b; // slopes
    vector[M-1] theta; // intercepts
    simplex[M-1] gamma; // gammas in Anh et al. 2009
}

transformed parameters{
    vector<lower=0, upper=1>[M] phi; // scale factors

    phi[1] <- 0;
    phi[2] <- gamma[1];
    phi[3] <- gamma[1] + gamma[2];
    phi[4] <- gamma[1] + gamma[2] + gamma[3];
}

model {
    vector[M] prob;
    vector[M] emu;
   
    // specify priors

    gamma ~ dirichlet(alpha);
    for (k in 1:K)
        b[k] ~ normal(0, 100);

    for (m  in 1:M-1)

        theta[m] ~ normal(0, 100);
    for (n in 1:N) {
            emu[1] <- exp(phi[1]*x[n]*b);
            for (m in 2:M)

                emu[m] <- exp(theta[m-1] + phi[m]*x[n]*b);

            for (m in 1:M)
                prob[m] <- emu[m]/sum(emu);
            y[n] ~ categorical(prob);
    }
}

generated quantities {
    vector[M] Gprob; // G for generated
    vector[M] Gemu;
    vector[N] log_lik;

    for (n in 1:N) {

        Gemu[1] <- exp(phi[1]*x[n]*b);

        for (m in 2:M)

            Gemu[m] <- exp(theta[m-1] + phi[m]*x[n]*b);

        for (m in 1:M)

            Gprob[m] <- Gemu[m]/sum(Gemu);
        log_lik[n] <- categorical_log(y[n], Gprob);

[Andrew Gelman]

Wow—cool!

Just a quick comment:  I suggest using normal(0,10) or even something more informative such as normal(0,1), rather than normal(0,100) which could create problems if you have weak data.  I’m assuming your x’s have been scaled so that the coefs are on unit scale.

A

[jun...@gmail.com]

Andrew,

Thanks for your comments and suggestions. Is it when we have weak data, we want to have strong priors to get convergence? I didn't scale x's. I know when the scale for some of the x's is large, I need to subtract their means/standardize for convergence/sampling purposes. But if the original scales are ok, do I have to rescale? Thanks a lot!

jun

[jun...@gmail.com]

I forgot to mention that the codes are for four-level response variable, and they should be adapted for variables with a greater number of response levels.

[Bob Carpenter]

Thanks for sharing.

Yes, it helps to rescale all the predictors so that
the coefficients are roughly unit scale. Otherwise, Stan needs
to learn those scales during warmup.

What matters is the posterior, which you can think of as
the combination of prior knowledge plus data. With little
data, you need more prior, as in hierarchical models or
something more constrained on an ordinary prior.

- Bob