Important: Read this before posting to this forum

  1. This forum is for questions related to the use of Apollo. We will answer some general choice modelling questions too, where appropriate, and time permitting. We cannot answer questions about how to estimate choice models with other software packages.
  2. There is a very detailed manual for Apollo available at http://www.ApolloChoiceModelling.com/manual.html. This contains detailed descriptions of the various Apollo functions, and numerous examples are available at http://www.ApolloChoiceModelling.com/examples.html. In addition, help files are available for all functions, using e.g. ?apollo_mnl
  3. Before asking a question on the forum, users are kindly requested to follow these steps:
    1. Check that the same issue has not already been addressed in the forum - there is a search tool.
    2. Ensure that the correct syntax has been used. For any function, detailed instructions are available directly in Apollo, e.g. by using ?apollo_mnl for apollo_mnl
    3. Check the frequently asked questions section on the Apollo website, which discusses some common issues/failures. Please see http://www.apollochoicemodelling.com/faq.html
    4. Make sure that R is using the latest official release of Apollo.
  4. If the above steps do not resolve the issue, then users should follow these steps when posting a question:
    1. provide full details on the issue, including the entire code and output, including any error messages
    2. posts will not immediately appear on the forum, but will be checked by a moderator first. This may take a day or two at busy times. There is no need to submit the post multiple times.

Scaling in joint model estimation with non-discrete dual response

Ask general questions about model specification and estimation that are not Apollo specific but relevant to Apollo users.
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Scaling in joint model estimation with non-discrete dual response

Post by cybey »

Hello, everyone,

I have a data set where respondents first had to choose between three alternatives (demand-side management programs; so-called forced choices) and then were asked whether they would actually participate in the alternative chosen before (using a Likert scale from 1 to 7; so-called free choices).

The list of utilities looks like this:

Code: Select all

### CBC (Forced Choices)
  V[['alt1']] = b_asc_1_value +
    b_Quota_Count_value * Quota_Count.1  +
    b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
    b_Quota_Compensation_value * Quota_Compensation.1
  
  V[['alt2']] = b_asc_2_value +
    b_Quota_Count_value * Quota_Count.2  +
    b_Quota_Max_Duration_value * Quota_Max_Duration.2  +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.2 + b_Quota_Period_afternoon_value * Quota_Period3.2 + b_Quota_Period_evening_value * Quota_Period4.2 + b_Quota_Period_night_value * Quota_Period5.2 +
    b_Quota_Compensation_value * Quota_Compensation.2
  
  V[['alt3']] = b_asc_3_value +
    b_Quota_Count_value * Quota_Count.3  +
    b_Quota_Max_Duration_value * Quota_Max_Duration.3  +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.3 + b_Quota_Period_afternoon_value * Quota_Period3.3 + b_Quota_Period_evening_value * Quota_Period4.3 + b_Quota_Period_night_value * Quota_Period5.3 +
    b_Quota_Compensation_value * Quota_Compensation.3
  
  V[['alt4']] = b_DR_asc_0_value
  
  V[['alt5']] = b_DR_asc_1_value +
    b_Quota_Count_value * Quota_Count.1  +
    b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
    b_Quota_Compensation_value * Quota_Compensation.1
  
  # Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = Choice,
    V             = lapply(V, "*", (1+mu_CBC)),
    rows          = (CBC_type=="R"),
    componentName = "CBC"
  )
  
  # Compute probabilities using MNL model
  P[['CBC']] = apollo_mnl(mnl_settings, functionality)
  
  
  ### DR (Free Choices)
  
  # Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt4=0, alt5=1),
    avail         = list(alt4=1, alt5=1),
    choiceVar     = Choice,
    V             = lapply(V, "*", (1+mu_DR)),
    rows          = (CBC_type=="RD"),
    componentName = "CBC_DR"
Now I was asking myself whether it makes sense to use scale parameters (mu_CBC and mu_DR) when merging the free and forced choice data? On the one hand, it makes sense, as one of the two choice processes could be more/less deterministic, resulting in different error terms. On the other hand, the error term is directly influenced by the cutoff value of the Likert scale. If, for example, I assume that a respondent chooses an alternative only at a value of 7, which is rather a conservative assumption, then I get a different error term than at a cutoff value of 6.

I look forward to your opinions!

Nico
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Scaling in joint model estimation with non-discrete dual response

Post by stephanehess »

Nico

why would you use the binary model for the free choice instead of an OL model?

A scale difference would indeed make sense

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Re: Scaling in joint model estimation with non-discrete dual response

Post by cybey »

A binary variable (participation/no participation) is easier to communicate than an ordinal variable (1 = "participation very unlikely" to 7 = "participation very likely"). The point is that demand-side management programs are something very new in practice and, to that extent, it is of course debatable whether someone who ticks a 5 on a Likert scale would actually participate in them. But yes, you're right, I've thought about an ordinal model too, as it is methodically more sound.

Is it possible to include the utility coefficients directly (with a scale parameter mu_DR) in the ordinal model or should a latent variable be used?

Approach 1: without a latent variable

Code: Select all

gamma_Quota_Count_WillingnessToParticipate = (1+mu_DR) * b_Quota_Count
  gamma_Quota_Max_Duration_WillingnessToParticipate = (1+mu_DR) * b_Quota_Max_Duration
  gamma_Quota_Period_morning_WillingnessToParticipate = (1+mu_DR) * b_Quota_Period_morning
  gamma_Quota_Period_midmorning_WillingnessToParticipate = (1+mu_DR) * b_Quota_Period_midmorning
  gamma_Quota_Period_afternoon_WillingnessToParticipate = (1+mu_DR) * b_Quota_Period_afternoon
  gamma_Quota_Period_evening_WillingnessToParticipate = (1+mu_DR) * b_Quota_Period_evening
  gamma_Quota_Period_night_WillingnessToParticipate = (1+mu_DR) * b_Quota_Period_night
  gamma_Quota_Compensation_WillingnessToParticipate = (1+mu_DR) * b_Quota_Compensation
  
  [...]
  
  ol_settings1 = list(outcomeOrdered=Choice,
                      V=### Attributes
                        gamma_Quota_Count_WillingnessToParticipate * Quota_Count.1 +
                        gamma_Quota_Max_Duration_WillingnessToParticipate * Quota_Max_Duration.1 +
                        gamma_Quota_Period_morning_WillingnessToParticipate * 0 +
                        gamma_Quota_Period_midmorning_WillingnessToParticipate * Quota_Period2.1 +
                        gamma_Quota_Period_afternoon_WillingnessToParticipate * Quota_Period3.1 +
                        gamma_Quota_Period_evening_WillingnessToParticipate * Quota_Period4.1 +
                        gamma_Quota_Period_night_WillingnessToParticipate * Quota_Period5.1 +
                        gamma_Quota_Compensation_WillingnessToParticipate * Quota_Compensation.1 +
                        
                        ### Sociodemographics
                        gamma_Gender_WillingnessToParticipate * COV_Gender +
                        gamma_Age_WillingnessToParticipate * COV_Age_centered +
                        gamma_Education_WillingnessToParticipate * COV_Education_centered +
                        gamma_Income_WillingnessToParticipate * COV_IncomeClasses_centered +
                        
                        ### Current conditions
                        gamma_ApartmentOwner_WillingnessToParticipate * COV_ApartmentOwner +
                        gamma_ExistingApps_WillingnessToParticipate * COV_ExistingApps +
                        gamma_CurrentMixs_WillingnessToParticipate * COV_CurrentMix,
                      tau=list(tau_WillingnessToParticipate_1, tau_WillingnessToParticipate_2, tau_WillingnessToParticipate_3, tau_WillingnessToParticipate_4, tau_WillingnessToParticipate_5, tau_WillingnessToParticipate_6),
                      rows=(CBC_type=="RD"),
                      componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_ol(ol_settings1, functionality)
  

Approach 2: with a latent variable

Code: Select all

# ----------------------------------------------------------------- #
  #---- Latent Variables
  # ----------------------------------------------------------------- #
  randcoeff[["LV_WillingnessToParticipate"]] =
    ### Attributes
    b_Quota_Count * Quota_Count.1 +
    b_Quota_Max_Duration * Quota_Max_Duration.1 +
    b_Quota_Period_morning * 0 +
    b_Quota_Period_midmorning * Quota_Period2.1 +
    b_Quota_Period_afternoon * Quota_Period3.1 +
    b_Quota_Period_evening * Quota_Period4.1 +
    b_Quota_Period_night * Quota_Period5.1 +
    b_Quota_Compensation * Quota_Compensation.1 +
    
    ### Sociodemographics
    gamma_Gender_WillingnessToParticipate * COV_Gender +
    gamma_Age_WillingnessToParticipate * COV_Age_centered +
    gamma_Education_WillingnessToParticipate * COV_Education_centered +
    gamma_Income_WillingnessToParticipate * COV_IncomeClasses_centered +
    
    ### Current conditions
    gamma_ApartmentOwner_WillingnessToParticipate * COV_ApartmentOwner +
    gamma_ExistingApps_WillingnessToParticipate * COV_ExistingApps +
    gamma_CurrentMixs_WillingnessToParticipate * COV_CurrentMix +
    
    ### Error term
    sigma_WillingnessToParticipate * eta_WillingnessToParticipate
    
    [...]
    
    ol_settings1 = list(outcomeOrdered=Choice,
                      V=zeta_WillingnessToParticipate_b_DR_asc_1 * LV_WillingnessToParticipate,
                      tau=list(tau_WillingnessToParticipate_1, tau_WillingnessToParticipate_2, tau_WillingnessToParticipate_3, tau_WillingnessToParticipate_4, tau_WillingnessToParticipate_5, tau_WillingnessToParticipate_6),
                      rows=(CBC_type=="RD"),
                      componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_ol(ol_settings1, functionality)
    
zeta_WillingnessToParticipate_b_DR_asc_1 is fixed at 1.

I think that approach 1 does not work because the individual-specific error term is missing, leading to a decrease in the beta coefficients, which now capture the unexplainable part (=error) of the measurement model. Or did is miss something?


Nico
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Scaling in joint model estimation with non-discrete dual response

Post by stephanehess »

Nico

I see your point, but if the survey used a Likert scale, then turning this into binary involves lots of approximations.

You don't need a LV for this, nor do you need to rewrite your utilities that way. You should just be able to reuse the relevant V from the MNL, but multiplied by a scale term

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Re: Scaling in joint model estimation with non-discrete dual response

Post by cybey »

Hello Stephane,

thank you for your feedback!
I tried both, an ordinal and a linear model. Could you please take a quick look at both specifications and tell me if they are correct?

Ordinal:

Code: Select all

b_DR_asc_1_value =
    ### Sociodemographics
    b_Gender_DR_asc_1 * COV_Gender +
    b_Age_DR_asc_1 * COV_Age_centered +
    b_Education_DR_asc_1 * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_DR_asc_1 * COV_IncomeClasses_centered +
    
    ### Current conditions
    b_ApartmentOwner_DR_asc_1 * COV_ApartmentOwner +
    b_ExistingApps_DR_asc_1 * COV_ExistingApps +
    b_CurrentMix_DR_asc_1 * COV_CurrentMix

Code: Select all

### DR (Free Choices)
  ol_settings1 = list(outcomeOrdered=Choice,
                      V = (1+mu_DR) * (b_Quota_Count_value * Quota_Count.1  +
                                         b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
                                         b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
                                         b_Quota_Compensation_value * Quota_Compensation.1 +
                                         b_DR_asc_1_value),
                      tau=list(tau_DR_asc_1_1, tau_DR_asc_1_2, tau_DR_asc_1_3, tau_DR_asc_1_4, tau_DR_asc_1_5, tau_DR_asc_1_6),
                      rows=(CBC_type=="RD"),
                      componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_ol(ol_settings1, functionality)

Linear:

Code: Select all

b_DR_asc_1_value = b_DR_asc_1 +
    ### Sociodemographics
    b_Gender_DR_asc_1 * COV_Gender +
    b_Age_DR_asc_1 * COV_Age_centered +
    b_Education_DR_asc_1 * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_DR_asc_1 * COV_IncomeClasses_centered +
    
    ### Current conditions
    b_ApartmentOwner_DR_asc_1 * COV_ApartmentOwner +
    b_ExistingApps_DR_asc_1 * COV_ExistingApps +
    b_CurrentMix_DR_asc_1 * COV_CurrentMix

Code: Select all

### DR (Free Choices)
  normalDensity_settings1 = list(outcomeNormal = Choice, 
                                 xNormal       = (1+mu_DR) * (b_Quota_Count_value * Quota_Count.1  +
                                                                b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
                                                                b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
                                                                b_Quota_Compensation_value * Quota_Compensation.1 +
                                                                b_DR_asc_1_value),
                                 mu            = 0, 
                                 sigma         = sigma_DR_asc_1,
                                 rows=(CBC_type=="RD"),
                                 componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_normalDensity(normalDensity_settings1, functionality)
Please note that in the linear specification, I used b_DR_asc_1 (intercept) to capture the generel tendency of respondents to participate in DSM programs. In addition, I mean-centered the responses to the Likert-scale questions.

Kind regards
Nico
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Scaling in joint model estimation with non-discrete dual response

Post by stephanehess »

Nico

can you confirm what is used in the first stage of the model, i.e. is it

Code: Select all

b_Quota_Count_value * Quota_Count.1  +
                                         b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
                                         b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
                                         b_Quota_Compensation_value * Quota_Compensation.1 +
                                         b_DR_asc_1_value
maybe show us the whole code?

Thanks

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Re: Scaling in joint model estimation with non-discrete dual response

Post by cybey »

Thank you for the fast response. In the following, as requested the code of the two MNL models.

Ordinal:

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)
library(alr4) # for deltaMethod

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName = "flexQgrid_MNL_OL1Ordinal",
  modelDescr = "MNL model in utility space (classical estimation)",
  indivID ="sys_RespNum",
  nCores = 3, # parallel :: detectCores()
  panelData = TRUE,
  mixing = FALSE
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #
input.file <- "./input/CBC_results_merged.wide.rds"
database <- readRDS(input.file)

### Dual response
# database[which(database$CBC_type == "RD" & database$Choice <6), "Choice"] <- 0
# database[which(database$CBC_type == "RD" & database$Choice >=6), "Choice"] <- 1
# database$Choice[database$CBC_type=="RD"] <- scale(database$Choice[database$CBC_type=="RD"], center=TRUE, scale=FALSE)

### Dummy for first row of a respondent
database[, "sys_RespNum.first"] <- rep(0, times=nrow(database))
database[match( unique(database$sys_RespNum), database$sys_RespNum), "sys_RespNum.first"] <- 1

### Mean centering of indicators
NonAttendance_colIndices <- grep("NonAttendance", colnames(database))
database[,NonAttendance_colIndices] <- apply(database[,NonAttendance_colIndices], 2, scale, center=TRUE, scale=FALSE)

### Rescale data
# old <- c(1,2,3,4,5)
# new <- c(30,60,90,120,150)
# 
# database$Quota_Compensation.1[database$Quota_Compensation.1 %in% old] <- new[match(database$Quota_Compensation.1, old, nomatch = 0)]
# database$Quota_Compensation.2[database$Quota_Compensation.2 %in% old] <- new[match(database$Quota_Compensation.2, old, nomatch = 0)]
# database$Quota_Compensation.3[database$Quota_Compensation.3 %in% old] <- new[match(database$Quota_Compensation.3, old, nomatch = 0)]

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation

## Dummy Coding + Constants
apollo_beta = c(
  # ----------------------------------------------------------------- #
  #---- Scaling parameter
  # ----------------------------------------------------------------- #
  mu_CBC = 0,
  mu_DR = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Alternative specific constants (ASC)
  # ----------------------------------------------------------------- #
  b_asc_1 = 0,
  b_asc_2 = 0,
  b_asc_3 = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_0
  # ----------------------------------------------------------------- #
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_1
  # ----------------------------------------------------------------- #
  ### Ordered measurement model
  tau_DR_asc_1_1 = -3,
  tau_DR_asc_1_2 = -2,
  tau_DR_asc_1_3 = -1,
  tau_DR_asc_1_4 = 1,
  tau_DR_asc_1_5 = 2,
  tau_DR_asc_1_6 = 3,
  
  # b_DR_asc_1 = 0,
  
  ### Sociodemographics
  b_Gender_DR_asc_1 = 0,
  b_Age_DR_asc_1 = 0,
  b_Education_DR_asc_1 = 0,
  b_Income_DR_asc_1 = 0,

  ### Current conditions
  b_ApartmentOwner_DR_asc_1 = 0,
  b_ExistingApps_DR_asc_1 = 0,
  b_CurrentMix_DR_asc_1 = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Count
  # ----------------------------------------------------------------- #
  b_Quota_Count = 0,
  
  b_Gender_Quota_Count = 0,
  b_Age_Quota_Count = 0,
  b_Education_Quota_Count = 0,
  b_ExistingApps_Quota_Count = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Max_Duration
  # ----------------------------------------------------------------- #
  b_Quota_Max_Duration = 0,
  
  b_Gender_Quota_Duration = 0,
  b_Age_Quota_Duration = 0,
  b_Education_Quota_Duration = 0,
  b_ExistingApps_Quota_Duration = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_morning
  # ----------------------------------------------------------------- #
  b_Quota_Period_morning = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_midmorning
  # ----------------------------------------------------------------- #
  b_Quota_Period_midmorning = 0,
  
  b_Gender_Quota_Period_midmorning = 0,
  b_Age_Quota_Period_midmorning = 0,
  b_Education_Quota_Period_midmorning = 0,
  b_ExistingApps_Quota_Period_midmorning = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_afternoon
  # ----------------------------------------------------------------- #
  b_Quota_Period_afternoon = 0,
  
  b_Gender_Quota_Period_afternoon = 0,
  b_Age_Quota_Period_afternoon = 0,
  b_Education_Quota_Period_afternoon = 0,
  b_ExistingApps_Quota_Period_afternoon = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_evening
  # ----------------------------------------------------------------- #
  b_Quota_Period_evening = 0,
  
  b_Gender_Quota_Period_evening = 0,
  b_Age_Quota_Period_evening = 0,
  b_Education_Quota_Period_evening = 0,
  b_ExistingApps_Quota_Period_evening = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_night
  # ----------------------------------------------------------------- #
  b_Quota_Period_night = 0,
  
  b_Gender_Quota_Period_night = 0,
  b_Age_Quota_Period_night = 0,
  b_Education_Quota_Period_night = 0,
  b_ExistingApps_Quota_Period_night = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Compensation
  # ----------------------------------------------------------------- #
  b_Quota_Compensation = -0.1,
  
  b_Gender_Quota_Compensation = 0,
  b_Age_Quota_Compensation = 0,
  b_Education_Quota_Compensation = 0,
  b_Income_Quota_Compensation = 0,
  
  b_CheapTalk_Quota_Compensation = 0
)


### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c("b_asc_3",
                 "b_Quota_Period_morning",
                 "mu_CBC")

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities.functionality = "estimate"

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality=apollo_probabilities.functionality) {
  
  ### Function initialisation: do not change the following three commands
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P = list()
  
  # ----------------------------------------------------------------- #
  #---- Scaling parameters
  # ----------------------------------------------------------------- #

  
  # ----------------------------------------------------------------- #
  #---- Alternative specific constants (ASC)
  # ----------------------------------------------------------------- #
  b_asc_1_value = b_asc_1
  b_asc_2_value = b_asc_2
  b_asc_3_value = b_asc_3

  
  # ----------------------------------------------------------------- #
  #---- DR_asc_0
  # ----------------------------------------------------------------- #
  
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_1
  # ----------------------------------------------------------------- #
  b_DR_asc_1_value =
    ### Sociodemographics
    b_Gender_DR_asc_1 * COV_Gender +
    b_Age_DR_asc_1 * COV_Age_centered +
    b_Education_DR_asc_1 * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_DR_asc_1 * COV_IncomeClasses_centered +
    
    ### Current conditions
    b_ApartmentOwner_DR_asc_1 * COV_ApartmentOwner +
    b_ExistingApps_DR_asc_1 * COV_ExistingApps +
    b_CurrentMix_DR_asc_1 * COV_CurrentMix
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Count
  # ----------------------------------------------------------------- #
  b_Quota_Count_value = b_Quota_Count +
    b_Gender_Quota_Count * COV_Gender +
    b_Age_Quota_Count * COV_Age_centered +
    b_Education_Quota_Count * COV_Education_centered +
    b_ExistingApps_Quota_Count * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Max_Duration
  # ----------------------------------------------------------------- #
  b_Quota_Max_Duration_value = b_Quota_Max_Duration +
    b_Gender_Quota_Duration * COV_Gender +
    b_Age_Quota_Duration * COV_Age_centered +
    b_Education_Quota_Duration * COV_Education_centered +
    b_ExistingApps_Quota_Duration * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_morning
  # ----------------------------------------------------------------- #
  b_Quota_Period_morning_value = b_Quota_Period_morning
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_midmorning
  # ----------------------------------------------------------------- #
  b_Quota_Period_midmorning_value = b_Quota_Period_midmorning +
    b_Gender_Quota_Period_midmorning * COV_Gender +
    b_Age_Quota_Period_midmorning * COV_Age_centered +
    b_Education_Quota_Period_midmorning * COV_Education_centered +
    b_ExistingApps_Quota_Period_midmorning * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_afternoon
  # ----------------------------------------------------------------- #
  b_Quota_Period_afternoon_value = b_Quota_Period_afternoon +
    b_Gender_Quota_Period_afternoon * COV_Gender +
    b_Age_Quota_Period_afternoon * COV_Age_centered +
    b_Education_Quota_Period_afternoon * COV_Education_centered +
    b_ExistingApps_Quota_Period_afternoon * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_evening
  # ----------------------------------------------------------------- #
  b_Quota_Period_evening_value = b_Quota_Period_evening +
    b_Gender_Quota_Period_evening * COV_Gender +
    b_Age_Quota_Period_evening * COV_Age_centered +
    b_Education_Quota_Period_evening * COV_Education_centered +
    b_ExistingApps_Quota_Period_evening * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_night
  # ----------------------------------------------------------------- #
  b_Quota_Period_night_value = b_Quota_Period_night +
    b_Gender_Quota_Period_night * COV_Gender +
    b_Age_Quota_Period_night * COV_Age_centered +
    b_Education_Quota_Period_night * COV_Education_centered +
    b_ExistingApps_Quota_Period_night * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Compensation
  # ----------------------------------------------------------------- #
  b_Quota_Compensation_value = b_Quota_Compensation +
    
    ### Sociodemographics
    b_Gender_Quota_Compensation * COV_Gender +
    b_Age_Quota_Compensation * COV_Age_centered +
    b_Education_Quota_Compensation * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_Quota_Compensation * COV_IncomeClasses_centered +
    b_CheapTalk_Quota_Compensation * COV_CheapTalk
  
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  
  ### CBC (Forced Choices)
  V[['alt1']] = b_asc_1_value +
    b_Quota_Count_value * Quota_Count.1 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.1 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
    b_Quota_Compensation_value * Quota_Compensation.1
  
  V[['alt2']] = b_asc_2_value +
    b_Quota_Count_value * Quota_Count.2 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.2 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.2 + b_Quota_Period_afternoon_value * Quota_Period3.2 + b_Quota_Period_evening_value * Quota_Period4.2 + b_Quota_Period_night_value * Quota_Period5.2 +
    b_Quota_Compensation_value * Quota_Compensation.2
  
  V[['alt3']] = b_asc_3_value +
    b_Quota_Count_value * Quota_Count.3 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.3 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.3 + b_Quota_Period_afternoon_value * Quota_Period3.3 + b_Quota_Period_evening_value * Quota_Period4.3 + b_Quota_Period_night_value * Quota_Period5.3 +
    b_Quota_Compensation_value * Quota_Compensation.3
  
  # Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = Choice,
    V             = lapply(V, "*", (1+mu_CBC)),
    rows          = (CBC_type=="R"),
    componentName = "CBC"
  )
  
  # Compute probabilities using MNL model
  P[['CBC']] = apollo_mnl(mnl_settings, functionality)
  
  
  ### DR (Free Choices)
  ol_settings1 = list(outcomeOrdered=Choice,
                      V = (1+mu_DR) * (b_Quota_Count_value * Quota_Count.1  +
                                         b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
                                         b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
                                         b_Quota_Compensation_value * Quota_Compensation.1 +
                                         b_DR_asc_1_value),
                      tau=list(tau_DR_asc_1_1, tau_DR_asc_1_2, tau_DR_asc_1_3, tau_DR_asc_1_4, tau_DR_asc_1_5, tau_DR_asc_1_6),
                      rows=(CBC_type=="RD"),
                      componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_ol(ol_settings1, functionality)
  
  
  ### Combined 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)
}

# ################################################################# #
#### MODEL ESTIMATION                                            ####
# ################################################################# #
model = apollo_estimate(apollo_beta,
                        apollo_fixed,
                        apollo_probabilities,
                        apollo_inputs)

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #

saveOutput_settings = list(
  printClassical = TRUE,
  printPVal = TRUE,
  printDiagnostics = TRUE,
  printCovar = TRUE,
  printCorr = TRUE,
  printOutliers = 50,
  printChange = TRUE,
  saveEst = TRUE,
  saveCov = TRUE,
  saveCorr = TRUE,
  saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)
Linear:

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)
library(alr4) # for deltaMethod

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName = "flexQgrid_MNL_OL1",
  modelDescr = "MNL model in utility space (classical estimation)",
  indivID ="sys_RespNum",
  nCores = 1, # parallel :: detectCores()
  panelData = TRUE,
  mixing = FALSE
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #
input.file <- "./input/CBC_results_merged.wide.rds"
database <- readRDS(input.file)

### Dual response
# database[which(database$CBC_type == "RD" & database$Choice <6), "Choice"] <- 0
# database[which(database$CBC_type == "RD" & database$Choice >=6), "Choice"] <- 1
database$Choice[database$CBC_type=="RD"] <- scale(database$Choice[database$CBC_type=="RD"], center=TRUE, scale=FALSE)

### Dummy for first row of a respondent
database[, "sys_RespNum.first"] <- rep(0, times=nrow(database))
database[match( unique(database$sys_RespNum), database$sys_RespNum), "sys_RespNum.first"] <- 1

### Mean centering of indicators
NonAttendance_colIndices <- grep("NonAttendance", colnames(database))
database[,NonAttendance_colIndices] <- apply(database[,NonAttendance_colIndices], 2, scale, center=TRUE, scale=FALSE)

### Rescale data
# old <- c(1,2,3,4,5)
# new <- c(30,60,90,120,150)
# 
# database$Quota_Compensation.1[database$Quota_Compensation.1 %in% old] <- new[match(database$Quota_Compensation.1, old, nomatch = 0)]
# database$Quota_Compensation.2[database$Quota_Compensation.2 %in% old] <- new[match(database$Quota_Compensation.2, old, nomatch = 0)]
# database$Quota_Compensation.3[database$Quota_Compensation.3 %in% old] <- new[match(database$Quota_Compensation.3, old, nomatch = 0)]

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation

## Dummy Coding + Constants
apollo_beta = c(
  # ----------------------------------------------------------------- #
  #---- Scaling parameter
  # ----------------------------------------------------------------- #
  mu_CBC = 0,
  mu_DR = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Alternative specific constants (ASC)
  # ----------------------------------------------------------------- #
  b_asc_1 = 0,
  b_asc_2 = 0,
  b_asc_3 = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_0
  # ----------------------------------------------------------------- #
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_1
  # ----------------------------------------------------------------- #
  b_DR_asc_1 = 0,
  
  ### Continuous measurement model
  sigma_DR_asc_1 = 1,
  
  ### Sociodemographics
  b_Gender_DR_asc_1 = 0,
  b_Age_DR_asc_1 = 0,
  b_Education_DR_asc_1 = 0,
  b_Income_DR_asc_1 = 0,
  
  ### Current conditions
  b_ApartmentOwner_DR_asc_1 = 0,
  b_ExistingApps_DR_asc_1 = 0,
  b_CurrentMix_DR_asc_1 = 0,
  # b_EnergyDecisionMaker = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Count
  # ----------------------------------------------------------------- #
  b_Quota_Count = 0,
  
  b_Gender_Quota_Count = 0,
  b_Age_Quota_Count = 0,
  b_Education_Quota_Count = 0,
  b_ExistingApps_Quota_Count = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Max_Duration
  # ----------------------------------------------------------------- #
  b_Quota_Max_Duration = 0,
  
  b_Gender_Quota_Duration = 0,
  b_Age_Quota_Duration = 0,
  b_Education_Quota_Duration = 0,
  b_ExistingApps_Quota_Duration = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_morning
  # ----------------------------------------------------------------- #
  b_Quota_Period_morning = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_midmorning
  # ----------------------------------------------------------------- #
  b_Quota_Period_midmorning = 0,
  
  b_Gender_Quota_Period_midmorning = 0,
  b_Age_Quota_Period_midmorning = 0,
  b_Education_Quota_Period_midmorning = 0,
  b_ExistingApps_Quota_Period_midmorning = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_afternoon
  # ----------------------------------------------------------------- #
  b_Quota_Period_afternoon = 0,
  
  b_Gender_Quota_Period_afternoon = 0,
  b_Age_Quota_Period_afternoon = 0,
  b_Education_Quota_Period_afternoon = 0,
  b_ExistingApps_Quota_Period_afternoon = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_evening
  # ----------------------------------------------------------------- #
  b_Quota_Period_evening = 0,
  
  b_Gender_Quota_Period_evening = 0,
  b_Age_Quota_Period_evening = 0,
  b_Education_Quota_Period_evening = 0,
  b_ExistingApps_Quota_Period_evening = 0,
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_night
  # ----------------------------------------------------------------- #
  b_Quota_Period_night = 0,
  
  b_Gender_Quota_Period_night = 0,
  b_Age_Quota_Period_night = 0,
  b_Education_Quota_Period_night = 0,
  b_ExistingApps_Quota_Period_night = 0,
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Compensation
  # ----------------------------------------------------------------- #
  b_Quota_Compensation = -0.1,
  
  b_Gender_Quota_Compensation = 0,
  b_Age_Quota_Compensation = 0,
  b_Education_Quota_Compensation = 0,
  b_Income_Quota_Compensation = 0,
  
  b_CheapTalk_Quota_Compensation = 0
)


### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c("b_asc_3",
                 "b_Quota_Period_morning",
                 "mu_CBC")

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities.functionality = "estimate"

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality=apollo_probabilities.functionality) {
  
  ### Function initialisation: do not change the following three commands
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P = list()
  
  # ----------------------------------------------------------------- #
  #---- Scaling parameters
  # ----------------------------------------------------------------- #

  
  # ----------------------------------------------------------------- #
  #---- Alternative specific constants (ASC)
  # ----------------------------------------------------------------- #
  b_asc_1_value = b_asc_1
  b_asc_2_value = b_asc_2
  b_asc_3_value = b_asc_3

  
  # ----------------------------------------------------------------- #
  #---- DR_asc_0
  # ----------------------------------------------------------------- #
  
  
  # ----------------------------------------------------------------- #
  #---- DR_asc_1
  # ----------------------------------------------------------------- #
  b_DR_asc_1_value = b_DR_asc_1 +
    ### Sociodemographics
    b_Gender_DR_asc_1 * COV_Gender +
    b_Age_DR_asc_1 * COV_Age_centered +
    b_Education_DR_asc_1 * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_DR_asc_1 * COV_IncomeClasses_centered +
    
    ### Current conditions
    b_ApartmentOwner_DR_asc_1 * COV_ApartmentOwner +
    b_ExistingApps_DR_asc_1 * COV_ExistingApps +
    b_CurrentMix_DR_asc_1 * COV_CurrentMix
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Count
  # ----------------------------------------------------------------- #
  b_Quota_Count_value = b_Quota_Count +
    b_Gender_Quota_Count * COV_Gender +
    b_Age_Quota_Count * COV_Age_centered +
    b_Education_Quota_Count * COV_Education_centered +
    b_ExistingApps_Quota_Count * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Max_Duration
  # ----------------------------------------------------------------- #
  b_Quota_Max_Duration_value = b_Quota_Max_Duration +
    b_Gender_Quota_Duration * COV_Gender +
    b_Age_Quota_Duration * COV_Age_centered +
    b_Education_Quota_Duration * COV_Education_centered +
    b_ExistingApps_Quota_Duration * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_morning
  # ----------------------------------------------------------------- #
  b_Quota_Period_morning_value = b_Quota_Period_morning
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_midmorning
  # ----------------------------------------------------------------- #
  b_Quota_Period_midmorning_value = b_Quota_Period_midmorning +
    b_Gender_Quota_Period_midmorning * COV_Gender +
    b_Age_Quota_Period_midmorning * COV_Age_centered +
    b_Education_Quota_Period_midmorning * COV_Education_centered +
    b_ExistingApps_Quota_Period_midmorning * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_afternoon
  # ----------------------------------------------------------------- #
  b_Quota_Period_afternoon_value = b_Quota_Period_afternoon +
    b_Gender_Quota_Period_afternoon * COV_Gender +
    b_Age_Quota_Period_afternoon * COV_Age_centered +
    b_Education_Quota_Period_afternoon * COV_Education_centered +
    b_ExistingApps_Quota_Period_afternoon * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_evening
  # ----------------------------------------------------------------- #
  b_Quota_Period_evening_value = b_Quota_Period_evening +
    b_Gender_Quota_Period_evening * COV_Gender +
    b_Age_Quota_Period_evening * COV_Age_centered +
    b_Education_Quota_Period_evening * COV_Education_centered +
    b_ExistingApps_Quota_Period_evening * COV_ExistingApps
  
  # ----------------------------------------------------------------- #
  #---- Quota_Period_night
  # ----------------------------------------------------------------- #
  b_Quota_Period_night_value = b_Quota_Period_night +
    b_Gender_Quota_Period_night * COV_Gender +
    b_Age_Quota_Period_night * COV_Age_centered +
    b_Education_Quota_Period_night * COV_Education_centered +
    b_ExistingApps_Quota_Period_night * COV_ExistingApps
  
  
  # ----------------------------------------------------------------- #
  #---- Quota_Compensation
  # ----------------------------------------------------------------- #
  b_Quota_Compensation_value = b_Quota_Compensation +
    
    ### Sociodemographics
    b_Gender_Quota_Compensation * COV_Gender +
    b_Age_Quota_Compensation * COV_Age_centered +
    b_Education_Quota_Compensation * COV_Education_centered +
    (1-IncomeQuan_miss) * b_Income_Quota_Compensation * COV_IncomeClasses_centered +
    b_CheapTalk_Quota_Compensation * COV_CheapTalk
  
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  
  ### CBC (Forced Choices)
  V[['alt1']] = b_asc_1_value +
    b_Quota_Count_value * Quota_Count.1 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.1 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
    b_Quota_Compensation_value * Quota_Compensation.1
  
  V[['alt2']] = b_asc_2_value +
    b_Quota_Count_value * Quota_Count.2 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.2 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.2 + b_Quota_Period_afternoon_value * Quota_Period3.2 + b_Quota_Period_evening_value * Quota_Period4.2 + b_Quota_Period_night_value * Quota_Period5.2 +
    b_Quota_Compensation_value * Quota_Compensation.2
  
  V[['alt3']] = b_asc_3_value +
    b_Quota_Count_value * Quota_Count.3 +
    b_Quota_Max_Duration_value * Quota_Max_Duration.3 +
    b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.3 + b_Quota_Period_afternoon_value * Quota_Period3.3 + b_Quota_Period_evening_value * Quota_Period4.3 + b_Quota_Period_night_value * Quota_Period5.3 +
    b_Quota_Compensation_value * Quota_Compensation.3
  
  # Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = list(alt1=1, alt2=1, alt3=1),
    choiceVar     = Choice,
    V             = lapply(V, "*", (1+mu_CBC)),
    rows          = (CBC_type=="R"),
    componentName = "CBC"
  )
  
  # Compute probabilities using MNL model
  P[['CBC']] = apollo_mnl(mnl_settings, functionality)
  
  
  ### DR (Free Choices)
  normalDensity_settings1 = list(outcomeNormal = Choice, 
                                 xNormal       = (1+mu_DR) * (b_Quota_Count_value * Quota_Count.1  +
                                                                b_Quota_Max_Duration_value * Quota_Max_Duration.1  +
                                                                b_Quota_Period_morning_value * 0 + b_Quota_Period_midmorning_value * Quota_Period2.1 + b_Quota_Period_afternoon_value * Quota_Period3.1 + b_Quota_Period_evening_value * Quota_Period4.1 + b_Quota_Period_night_value * Quota_Period5.1 +
                                                                b_Quota_Compensation_value * Quota_Compensation.1 +
                                                                b_DR_asc_1_value),
                                 mu            = 0, 
                                 sigma         = sigma_DR_asc_1,
                                 rows=(CBC_type=="RD"),
                                 componentName = "OL_DualResponse")
  
  P[["OL_DualResponse"]] = apollo_normalDensity(normalDensity_settings1, functionality)
  
  
  ### Combined 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)
}

# ################################################################# #
#### MODEL ESTIMATION                                            ####
# ################################################################# #
model = apollo_estimate(apollo_beta,
                        apollo_fixed,
                        apollo_probabilities,
                        apollo_inputs)

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #

saveOutput_settings = list(
  printClassical = TRUE,
  printPVal = TRUE,
  printDiagnostics = TRUE,
  printCovar = TRUE,
  printCorr = TRUE,
  printOutliers = 50,
  printChange = TRUE,
  saveEst = TRUE,
  saveCov = TRUE,
  saveCorr = TRUE,
  saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Scaling in joint model estimation with non-discrete dual response

Post by stephanehess »

Nico

so what you're doing is sharing the part of the utility between stages that does not include the ASC, and you use a separate ASC in the two stages? If that's the case, there is no need to include the ASC in the multiplication by the scale parameter. In your linear model, if you have mean centred the response, then you should not be estimating b_DR_asc_1, I believe

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Re: Scaling in joint model estimation with non-discrete dual response

Post by cybey »

Yes, exactely. In the forced-choices (3 alternatives per choice set) I check for left-right effects (b_asc_1, b_asc_2, b_asc_3), whereas in the free choices, I would like to explain the responses to the Likert scale questions.

About the scale parameter: Does it matter if I multiply only the coefficients shared by two model components or all coefficients by the scale parameter? I checked the results on both approaches and they are the same (log likelihood value). This is also consistent with Bradley and Daly's (1997)* statement that the scale is arbitrary and the scale parameter only needs to be considered when interpreting the results. The reason I multiply all coefficient in the linear and ordinal model, respectively, is consistency, as I use the lapply() function in the MNL model component and multiply the whole utility function as well.

And thank you for the advice regarding the intercept. For ease of interpretation, I decided to a use a model with an intercept.

*Bradley, Mark; Daly, Andrew (1997): Estimation of Logit Choice Models Using Mixed Stated-Preference and Revealed-Preference Information. In Peter R. Stopher, Martin Lee-Gosselin (Eds.): Understanding travel behaviour in an era of change. 1st ed. Oxford: Pergamon, pp. 209–231.
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Scaling in joint model estimation with non-discrete dual response

Post by stephanehess »

Nico

yes, it's fine to include the parameter in the scale multiplication, just not necessary per se

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
Post Reply