From bc0b38850e5d528b7c38a08014d8c1e736421161 Mon Sep 17 00:00:00 2001
From: Fabian Marder <fabian.marder@idiv.de>
Date: Fri, 24 Mar 2023 15:51:57 +0000
Subject: [PATCH] Upload New File
---
R_Scripts/mnl_CMI.R | 160 ++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 160 insertions(+)
create mode 100644 R_Scripts/mnl_CMI.R
diff --git a/R_Scripts/mnl_CMI.R b/R_Scripts/mnl_CMI.R
new file mode 100644
index 0000000..17dc950
--- /dev/null
+++ b/R_Scripts/mnl_CMI.R
@@ -0,0 +1,160 @@
+### Initialise code
+apollo_initialise()
+
+### Set core controls
+apollo_control = list(
+ modelName = "MNL_SP",
+ modelDescr = "Simple MNL model on mode choice SP data",
+ indivID = "id",
+ outputDirectory = "Output",
+ panelData = TRUE
+)
+
+
+
+# ################################################################# #
+#### DEFINE MODEL PARAMETERS ####
+# ################################################################# #
+
+### Vector of parameters, including any that are kept fixed in estimation
+apollo_beta=c(asc_a = 0,
+ asc_b = 0,
+ b_profit =1,
+ b_hold = 0,
+ b_loc = 0,
+ b_admin1 = 0,
+ b_admin2 = 0,
+ b_hold_CMI = 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()
+
+# ################################################################# #
+#### 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()
+
+ ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
+ V = list()
+ V[["alt1"]] = asc_a +
+ b_profit * profit1 -
+ b_hold * hold1 -
+ b_hold_CMI*hold1*CMI+
+ b_loc * loc1 +
+ b_admin1 *admin_11 +
+ b_admin2 *admin_21
+
+
+
+ V[["alt2"]] =asc_b +
+ b_profit * profit2 -
+ b_hold * hold2 -
+ b_hold_CMI*hold2*CMI+
+ b_loc * loc2 +
+ b_admin1 *admin_12 +
+ b_admin2 *admin_22
+
+
+
+
+ V[["alt3"]] = b_profit * profit3
+
+ ### Define settings for MNL model component
+ mnl_settings = list(
+ alternatives = c(alt1=1, alt2=2, alt3=3),
+ avail = 1,
+ choiceVar = choosen,
+ utilities = V
+ )
+
+ ### 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 ####
+# ################################################################# #
+
+model_CMI = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
+
+# ################################################################# #
+#### MODEL OUTPUTS ####
+# ################################################################# #
+
+# ----------------------------------------------------------------- #
+#---- FORMATTED OUTPUT (TO SCREEN) ----
+# ----------------------------------------------------------------- #
+
+apollo_modelOutput(model_CMI)
+
+# ----------------------------------------------------------------- #
+#---- FORMATTED OUTPUT (TO FILE, using model name) ----
+# ----------------------------------------------------------------- #
+summary(model_CMI)
+
+deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold/b_profit)*100")))
+apollo_deltaMethod(model_CMI, deltaMethod_settings)
+deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold_CMI/b_profit)*100")))
+apollo_deltaMethod(model_CMI, deltaMethod_settings)
+
+min<-min(database$CMI)
+max<-max(database$CMI)
+x_list=seq(min,max,0.1)
+
+df <- data.frame(Mean=double(),
+ Index=double(),
+ UB=double(),
+ LB=double())
+i=1
+for(j in x_list){
+ deltaMethod_settings=list(expression=c(paste0(mean_rate="100*(b_hold/b_profit+b_hold_CMI/b_profit*",j,")")))
+ a<-apollo_deltaMethod(model_CMI, deltaMethod_settings)
+ df[i,"Mean"]<-a[["Value"]]
+ df[i,"Index"]<-j
+ df[i,"UB"]<-a[["Value"]]+1.96*a[["Robust s.e."]]
+ df[i,"LB"]<-a[["Value"]]-1.96*a[["Robust s.e."]]
+ i=i+1
+}
+
+
+fig_CMI<-ggplot(df, aes(x=Index, y=Mean)) +
+ geom_ribbon(aes(ymin=LB, ymax=UB),alpha=.2, fill = "#F8766D") +
+ geom_line(color = "#F8766D") +
+ scale_color_manual(values="blue")+
+ labs(x = "CMI", y= "Money discount rate in %")
+
+ggsave("Graphs/CMI.png")
+ggsave("Graphs/CMI.tiff")
+
+
+
+ggarrange(fig_PMI, fig_CMI,
+ labels = c("a", "b"),
+ ncol = 2, nrow = 1)
+ggsave("Graphs/PMI_CMI.tiff")
+
--
GitLab