Multivariate Ordered Logit Model

[punyabeet]

Hi Prof. Hess,

I have the following doubts on which I need your help.

1. Below attached is the code and output of the Multi-variate Ordered Logit Model. Is the model specification and output correct?

2. How does apollo_combineModels work ?? Can you please explain in brief and please suggest some papers regarding the methodology based on which apollo_combineModels code is developed. This would help me in understanding as well as citing it.

3. For a 5 episode activity, instead of taking starting tau values as -2, -1, 1, 2 can I specify -5, -1, 1, 5 ?? Will these changes impact the coefficients of explanatory variables ??

Thank you for your time and help.

[stephanehess]

Hi

1. Your model is simply combining a number of different ordered logit models. You are not making any links between these components, so I wouldn't really refer to it as multivariate. If you want to create links between the models, you could incorporate a latent variable, for example.

2. The function apollo_combineModels simply creates a joint likelihood by multiplying the likelihoods for the individual model components together. This is explained at the start of Chapter 7 of the manual.

3. Regarding tau, these are starting values, so they will change during estimation anyway and the model, if well behaved, should converge to the same solution as long as you use reasonable starting values. If you use too wide a range, extreme values for the dependent variable will have probabilities near or equal to zero at the starting values, likely leading to a failure of the code.

Stephane

[jjanmaat]

I've got a SP off RP style model that has a multivariate ordered RP component. It is an SP off RP style as the 'RP' component is survey participants choosing a future they believe the most likely from a set of possible futures, and then the subsequent SP experiment is pivoted off of that chosen future. The possible futures are described by three attributes, the same three attributes that are used, together with a payment vehicle, in the SP.

I've run a multivariate ordered model using mvord in R, and find that several individual specific variables are important in predicting the likelihood of an individual choosing a particular future. The probability distribution across the SP alternatives is therefore conditional on these predictors. So, I am thinking I need to combine models in the spirit of Train and Wilson's approach, but rather than having an RP that is based on choices driven by the attribute levels, mine has the 'SF' (stated future) based on individual specific variables.

Any suggestions? Relevant references as well as tips for running something like this in apollo would be most welcome.

Thanks.

John.

[stephanehess]

John

sorry, slightly unclear what you are asking. Could you maybe explain the model specification mathematically?

Best wishes

Stephane

[jjanmaat]

Hello Stephane,

Thanks for reading my post, and apologies for the lack of clarity.

I have attached a document setting out my current theoretical description of the modeling problem. I hope this is a bit more clear.

After you look at the attachment, my question might be more clear. Specifically:

I have apollo performing the conditional estimation, a multinomial logit based on the participants stated future. I don't know if I can model the participants choice of the stated future in apollo.

Thanks,

John.

[dpalma]

Hi John,

Sorry for the belated response. If I understand correctly, what you want to do is a choice experiment in two stages. In the first stage the respondent chooses an alternative, say one between A1, A2 or A3. Let’s imagine the respondent chooses A3. In the second stage, the same respondent faces a series of choices where A3 is always an alternative. For example, they may face a choice between B1, B2 and A3, and then another choice between C1, C2 and A3, and so forth.

As you mention in the uploaded document, the issue here is that the error term in the utility of the chosen alternative in the first stage (A3 in the example) is present in all subsequent choices, therefore introducing correlation between them. Sadly, there is no pre-coded function to do this in Apollo.

The simpler way I can think of implementing this using Apollo is using a multivariate probit model. You would need to figure out the correlation matrix between choices, and then code the likelihood function of a multivariate probit using a multivariate normal cdf (you can use the package mvtnorm in R).

Including the pre-coded likelihood of multivariate probit models might be an interesting addition to Apollo, but I am afraid it is not among our priorities right now, sorry.

I hope this helps a bit
David

[jjanmaat]

Hello David,

Given that it took me till today to check the forum, no apology necessary for the belated response. And thanks for your reply.

My tentative way forward is to use the R package 'mvord' to estimate a multivariate ordinal model of the choice among futures, and then use the predicted probabilities as weights in the multinomial regression in Apollo. If I that the utility impact of the choice relative to the stated future is independent of the choice of stated future, then this should not be an issue.

We are considering a hybrid model, where we have a couple of latent variables correlated with some attitude measures. There is suggestive evidence that these latent variables are predictive of both the chosen future and the choice conditional on the chosen future. This will be our big challenge!!!

Thanks again for your reply.

John.

[jjanmaat]

PS:

I have been trying to use Apollo to reproduce a multivariate ordered probit, using a Choleski representation of the covariance matrix between the components. I am successful in running the model, but it chokes on generating the summary. I have set the bottom break points for the ordered probit to zero. Not sure if this is what is leading to the error:

Error in varcovar[activePar, activePar] <- solve(-hessian(object)[activePar, :
number of items to replace is not a multiple of replacement length


Code components for null model - intercepts only. Aim is to enrich it with exogenous predictors, and then combine it with a multinomial probit/logit for the SP of SF choices. The coefficients make sense, relative to those generated by the mvord package, excepting some scaling.

Comments most welcome.

John.

#
# Set up starting values for parameters to estimate
apollo_beta <- c(
# break points for ordered variables
tau_Wea_1 = 0,
tau_Wea_2 = 1,
tau_Wea_3 = 2,
tau_Spe_1 = 0,
tau_Spe_2 = 1,
tau_Spe_3 = 2,
tau_Nat_1 = 0,
tau_Nat_2 = 1,
tau_Nat_3 = 2,
tau_Nat_4 = 3,
# parameters for predictive elements
gamma_Wea_0 = 0,
gamma_Spe_0 = 0,
gamma_Nat_0 = 0,
# Components of choleski decomposition
Ch11 = 1, Ch12 = 0, Ch13 = 0,
Ch22 = 1, Ch23 = 0, Ch33 = 1)
#
# Fixed parameters, to force identification.
apollo_fixed <- c("tau_Wea_1", "tau_Spe_1", "tau_Nat_1")
#
# Set up the form of the maximum simulated likelihood
apollo_draws <- list(
interDrawsType = "halton",
interNDraws = 100,
interNormDraws = c("eta_Wea", "eta_Spe", "eta_Nat")
)
#
# Set up random coefficients definitions
apollo_randCoeff <- function(apollo_beta, apollo_inputs) {
randcoef <- list()
randcoef[['LV_Wea']] <-
gamma_Wea_0 + eta_Wea*Ch11 + eta_Spe*Ch12 + eta_Nat*Ch13
randcoef[['LV_Spe']] <-
gamma_Spe_0 + eta_Spe*Ch22 + eta_Nat*Ch23
randcoef[['LV_Nat']] <-
gamma_Nat_0 + eta_Nat*Ch33
return(randcoef)
}

....

#
# Set up the probabilities for the choices.
apollo_probabilities <- function(apollo_beta,
apollo_inputs,
functionality ="estimate") {
### Attach inputs and detach after function exit
apollo_attach (apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
### Create list for probabilities
P <- list()
### Likelihood of choosing future
op_settings_Wea <- list(outcomeOrdered = WeaExp,
V = LV_Wea,
tau = c(tau_Wea_1, tau_Wea_2,
tau_Wea_3),
rows = (Task==1),
coding = 1:4,
componentName = "Ord_Wea")
op_settings_Spe <- list(outcomeOrdered = SpeExp,
V = LV_Spe,
tau = c(tau_Spe_1, tau_Spe_2,
tau_Spe_3),
rows = (Task==1),
coding = 1:4,
componentName = "Ord_Spe")
op_settings_Nat <- list(outcomeOrdered = NatExp,
V = LV_Nat,
tau = c(tau_Nat_1, tau_Nat_2,
tau_Nat_3, tau_Nat_4),
rows = (Task==1),
coding = 1:5,
componentName = "Ord_Nat")
#
P[['WeaLV']] <- apollo_op(op_settings_Wea, functionality)
P[['SpeLV']] <- apollo_op(op_settings_Spe, functionality)
P[['NatLV']] <- apollo_op(op_settings_Nat, functionality)
#
### Likelihood of the whole model
P <- apollo_combineModels(P, apollo_inputs, functionality)
### Take product across observation for same individual
P <- apollo_panelProd(P, apollo_inputs, functionality)
### Average across inter-individual draws
P <- apollo_avgInterDraws(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P <- apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}

[dpalma]

Hi John,

I am not familiar with the mvord package, nor with the weighting strategy to correct for the correlation between alternatives, but I am glad you found an approach that works and serves your problem.

About the code for the multiple ordered models with correlated utilities, I do not know where the problem is. I generated a synthetic database and run your code, and it run without issues. Results are obviously spurious as the data is synthetic, but the model does not fail. The code is at the bottom of this post.

Maybe you are using an old version of Apollo? Current version is 0.2.3. Also, as you will note in the code below, it is no longer necessary to include componentName among the settings of a model component, Apollo will use the name assigned in P.

Hope updating Apollo solves the issue.

Code: Select all

### Initialise
rm(list = ls())
library(apollo)
apollo_initialise()
apollo_control = list(
  modelName  ="OP_test",
  indivID    ="id", 
  mixing=TRUE)

### Load data
N <- 200
M <- 5
database <- data.frame(WeaExp= rep(sample(1:4, N, replace=TRUE), each=M), 
                       SpeExp= rep(sample(1:4, N, replace=TRUE), each=M), 
                       NatExp= rep(sample(1:5, N, replace=TRUE), each=M), 
                       Task  = rep(1:M, N), 
                       id    = rep(1:N, each=M))
rm(N, M)

### Set up starting values for parameters to estimate
apollo_beta <- c(
  # break points for ordered variables
  tau_Wea_1 = 0,
  tau_Wea_2 = 1,
  tau_Wea_3 = 2,
  tau_Spe_1 = 0,
  tau_Spe_2 = 1,
  tau_Spe_3 = 2,
  tau_Nat_1 = 0,
  tau_Nat_2 = 1,
  tau_Nat_3 = 2,
  tau_Nat_4 = 3,
  # parameters for predictive elements
  gamma_Wea_0 = 0,
  gamma_Spe_0 = 0,
  gamma_Nat_0 = 0,
  # Components of choleski decomposition
  Ch11 = 1, Ch12 = 0, Ch13 = 0,
  Ch22 = 1, Ch23 = 0, Ch33 = 1)
# Fixed parameters, to force identification.
apollo_fixed <- c("tau_Wea_1", "tau_Spe_1", "tau_Nat_1")

# Set up the form of the maximum simulated likelihood
apollo_draws <- list(
  interDrawsType = "halton",
  interNDraws = 100,
  interNormDraws = c("eta_Wea", "eta_Spe", "eta_Nat")
)

# Set up random coefficients definitions
apollo_randCoeff <- function(apollo_beta, apollo_inputs) {
  randcoef <- list()
  randcoef[['LV_Wea']] <- gamma_Wea_0 + eta_Wea*Ch11 + eta_Spe*Ch12 + eta_Nat*Ch13
  randcoef[['LV_Spe']] <- gamma_Spe_0 +                eta_Spe*Ch22 + eta_Nat*Ch23
  randcoef[['LV_Nat']] <- gamma_Nat_0 +                               eta_Nat*Ch33
  return(randcoef)
}

### Group and validate input
apollo_inputs <- apollo_validateInputs()


### Set up the probabilities for the choices.
apollo_probabilities <- function(apollo_beta,
                                 apollo_inputs,
                                 functionality ="estimate") {
  ### Attach inputs and detach after function exit
  apollo_attach (apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  ### Create list for probabilities
  P <- list()
  ### Likelihood of choosing future
  op_settings_Wea <- list(outcomeOrdered = WeaExp,
                          V = LV_Wea,
                          tau = c(tau_Wea_1, tau_Wea_2, tau_Wea_3),
                          rows = (Task==1),
                          coding = 1:4)
  op_settings_Spe <- list(outcomeOrdered = SpeExp,
                          V = LV_Spe,
                          tau = c(tau_Spe_1, tau_Spe_2, tau_Spe_3),
                          rows = (Task==1),
                          coding = 1:4)
  op_settings_Nat <- list(outcomeOrdered = NatExp,
                          V = LV_Nat,
                          tau = c(tau_Nat_1, tau_Nat_2, tau_Nat_3, tau_Nat_4),
                          rows = (Task==1),
                          coding = 1:5)
  #
  P[['WeaLV']] <- apollo_op(op_settings_Wea, functionality)
  P[['SpeLV']] <- apollo_op(op_settings_Spe, functionality)
  P[['NatLV']] <- apollo_op(op_settings_Nat, functionality)
  #
  ### Likelihood of the whole model
  P <- apollo_combineModels(P, apollo_inputs, functionality)
  ### Take product across observation for same individual
  P <- apollo_panelProd(P, apollo_inputs, functionality)
  ### Average across inter-individual draws
  P <- apollo_avgInterDraws(P, apollo_inputs, functionality)
  ### Prepare and return outputs of function
  P <- apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}

L <- apollo_probabilities(apollo_beta, apollo_inputs)

model = apollo_estimate(apollo_beta, apollo_fixed, 
                        apollo_probabilities, apollo_inputs,
                        estimate_settings=list(writeIter=FALSE))

apollo_modelOutput(model)

[jjanmaat]

Thanks for the example. It works on my computer. However, the error still occurs with my data.

To check on things, installed the latest version of R, updated apollo and all the other packages, and ran the analysis without loading any other packages. Also removed everything else from the database that was not used in the estimation, and made sure that the id values were a sequence with no gaps. None of this stopped the error.

I have attached the script I am using, and the stripped down dataset.

Any suggestions would be welcome.

John.