rm(list = ls()) devtools::load_all() library(dplyr) designpath <- system.file("extdata", "spdesigns", "designs", package = "simulateDCE") destype <- "spdesign" resps <- 4000 # number of respondents nosim <- 2 # number of simulations to run (about 500 is minimum) # betacoefficients should not include "-" bcoeff <- list( bx1A = 0.1, ## very high asc bx1B = -0.1, bx2A = 0.2, bx2B = -0.2 ) desisiongroups <- c(0, 0.3, 1) ul <- list( uA = list( v1 = V.1 ~ bx1A * alt1.x1 + bx2A * alt1.x2, v2 = V.2 ~ bx1A * alt2.x1 + bx2A * alt2.x2 ), uB = list( v1 = V.1 ~ bx1B * alt1.x1 + bx2B * alt1.x2, v2 = V.2 ~ bx1B * alt2.x1 + bx2B * alt2.x2 ) ) simplesim <- sim_all(nosim = nosim, resps = resps, designpath = designpath, bcoeff = bcoeff, u = ul, designtype = "spdesign", decisiongroups = desisiongroups) testdata <- simplesim[["twoattr"]][[1]][["data"]] %>% distinct(group, Choice_situation, .keep_all = T) %>% arrange(Choice_situation) ## estimate models library(apollo) apollo_initialise() apollo_control <- list( modelName = "Clogit_simpledesign", modelDescr = "Simple conditional logit model for group 1 ", indivID = "ID", nCores = 1, outputDirectory = "/~" ) ## group 1 database <- simplesim[["twoattr"]][[1]][["data"]] apollo_beta <- c( b_x1 = 0, b_x2 = 0 ) apollo_fixed <- c() # ################################################################# # #### GROUP AND VALIDATE INPUTS #### # ################################################################# # apollo_inputs <- apollo_validateInputs() # ################################################################# # #### DEFINE MODEL AND LIKELIHOOD FUNCTION #### # ################################################################# # 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 of probabilities P P <- list() ### Define settings for MNL model component that are generic across classes mnl_settings <- list( alternatives = c(alt1 = 1, alt2 = 2), avail = list(alt1 = 1, alt2 = 1), choiceVar = CHOICE ) ### Compute class-specific utilities V <- list() V[["alt1"]] <- b_x1 * alt1_x1 + b_x2 * alt1_x2 V[["alt2"]] <- b_x1 * alt2_x1 + b_x2 * alt2_x2 mnl_settings$utilities <- V # mnl_settings$componentName = paste0("Class_",s) ### Compute probabilities using MNL model P[["model"]] <- apollo_mnl(mnl_settings, functionality) ### Take product across observation for same individual P <- apollo_panelProd(P, apollo_inputs, functionality) ### Prepare and return outputs of function P <- apollo_prepareProb(P, apollo_inputs, functionality) return(P) } # ################################################################# # #### MODEL ESTIMATION #### # ################################################################# # ### Optional starting values search # apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs) model <- apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings = list(estimationRoutine = "BFGS", hessianRoutine = "numDeriv")) apollo_modelOutput(model) ### latent class model apollo_control <- list( modelName = "Latent Class logit_simpledesign", modelDescr = "Latent Class logit model", indivID = "ID", nCores = 1, outputDirectory = "/~" ) ## group 1 database <- simplesim[["twoattr"]][[1]][["data"]] apollo_beta <- c( b_x1_a = 0, b_x2_a = 0, b_x1_b = 0, b_x2_b = 0, delta_a = 0, delta_b = 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("delta_b") # ################################################################# # #### DEFINE LATENT CLASS COMPONENTS #### # ################################################################# # apollo_lcPars <- function(apollo_beta, apollo_inputs) { lcpars <- list() lcpars[["b_x1"]] <- list(b_x1_a, b_x1_b) lcpars[["b_x2"]] <- list(b_x2_a, b_x2_b) V <- list() V[["class_a"]] <- delta_a V[["class_b"]] <- delta_b classAlloc_settings <- list( classes = c(class_a = 1, class_b = 2), utilities = V ) lcpars[["pi_values"]] <- apollo_classAlloc(classAlloc_settings) return(lcpars) } # ################################################################# # #### GROUP AND VALIDATE INPUTS #### # ################################################################# # apollo_inputs <- apollo_validateInputs() # ################################################################# # #### DEFINE MODEL AND LIKELIHOOD FUNCTION #### # ################################################################# # # ################################################################# # #### DEFINE MODEL AND LIKELIHOOD FUNCTION #### # ################################################################# # 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 of probabilities P P <- list() ### Define settings for MNL model component that are generic across classes mnl_settings <- list( alternatives = c(alt1 = 1, alt2 = 2), avail = list(alt1 = 1, alt2 = 1), choiceVar = CHOICE ) ### Loop over classes for (s in 1:2) { ### Compute class-specific utilities V <- list() V <- list() V[["alt1"]] <- b_x1[[s]] * alt1_x1 + b_x2[[s]] * alt1_x2 V[["alt2"]] <- b_x1[[s]] * alt2_x1 + b_x2[[s]] * alt2_x2 mnl_settings$utilities <- V # mnl_settings$componentName = paste0("Class_",s) ### Compute within-class choice probabilities using MNL model P[[paste0("Class_", s)]] <- apollo_mnl(mnl_settings, functionality) ### Take product across observation for same individual P[[paste0("Class_", s)]] <- apollo_panelProd(P[[paste0("Class_", s)]], apollo_inputs, functionality) } ### Compute latent class model probabilities lc_settings <- list(inClassProb = P, classProb = pi_values) P[["model"]] <- apollo_lc(lc_settings, apollo_inputs, functionality) ### Prepare and return outputs of function P <- apollo_prepareProb(P, apollo_inputs, functionality) return(P) } # ################################################################# # #### MODEL ESTIMATION #### # ################################################################# # ### Optional starting values search apollo_beta <- apollo_searchStart(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs) model <- apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs) # ,estimate_settings =list(estimationRoutine="bhhh" , hessianRoutine= "maxLik") # ################################################################# # #### MODEL OUTPUTS #### # ################################################################# # # ----------------------------------------------------------------- # #---- FORMATTED OUTPUT (TO SCREEN) ---- # ----------------------------------------------------------------- # apollo_modelOutput(model)