two new examples added

Projects/Tofu/parameters_tofu.R

+
+
+designpath<- "Projects/Tofu/Designs/"
+
+resps =360  # number of respondents
+nosim=2 # number of simulations to run (about 500 is minimum)
+
+#betacoefficients should not include "-"
+basc = -2.5
+bcost = -0.05
+bmanufEU = -0.4
+bmanufnonEU = -1
+bcultEU = -0.2
+bcultnonEU = -0.4
+bCONV = -0.4
+
+
+
+manipulations = list(alt1.professional=     expr(alt1.initiator==1), 
+                     alt2.professional=     expr(alt2.initiator==1),
+                     alt1.expert      =     expr(alt1.initiator==2), 
+                     alt2.expert      =     expr(alt2.initiator==2),
+                     alt1.domestic    =     expr(alt1.funding==1), 
+                     alt2.domestic    =     expr(alt2.funding==1),
+                     alt1.foreign     =     expr(alt1.funding==2), 
+                     alt2.foreign     =     expr(alt2.funding==2))
+
+
+
+
+
+
+#place your utility functions here
+u<- list(u1= 
+  list(
+  v1 =V.1 ~  bprof*alt1.professional+ bexp * alt1.expert + bdomestic * alt1.domestic + bforeign * alt1.foreign + bdamage*alt1.damage + bprice * alt1.compensation,
+  v2 =V.2 ~  bprof*alt2.professional + bexp * alt2.expert + bdomestic * alt2.domestic + bforeign * alt2.foreign + bdamage*alt2.damage + bprice * alt2.compensation,
+  v3 =V.3 ~ basc)
+)

Projects/Tofu/simulation_SWE_after_pilot.Rmd

+---
+title: "Pre-Registration and Simulation after pretest"
+author: "Julian"
+date: "11 Feb 2021"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+t<- knitr::opts_chunk$set(echo = TRUE)
+source("functions/help_functions.R")
+```
+
+# Introduction
+
+This document has two purposes. First, it simulates the data for the preregistration at OFS with ID XXX, and second, it describes and exemplifies the analysis of the data as described in the preregistration. The document contains the r code for the simulation and the analysis, and provides details on each step. The user is free to change some parameters and values to inspect the limits of the design. Here, we simulate the data with a rather small sample size and realistic parameter values. While this html file is static, all required files to replicate and adapt this simulation are stored on github under the link
+
+
+# Preparation and base dataset
+
+In the following R chunk, we read in all packages needed (including noting down their versions). It also reads in functions written by the authors to automate the simulation and read in the data. The functions are stored under "functions/help_functions.R". Further, it generates objects which will be used to store runs from each run in the simulation. 
+
+```{r prep}
+
+rm(list = ls())
+
+
+library("dplyr")          
+library("evd")           
+library("boot")
+library("apollo")
+library("sjlabelled")
+library("readxl")
+library("httr")
+library("ggplot2")
+library("tidyr")
+library("magrittr")
+library("psych")
+library("stringr")
+dce_data <- read_excel("simulated_data/simulation_SWE_after_pilot/tofu_main_swe_sim02_tofu.xlsx")
+dce_dictionary <- read_excel("simulated_data/simulation_SWE_after_pilot/tofu_main_swe_sim02_tofu.xlsx", sheet = "dictionary")
+
+
+source("functions/help_functions.R") 
+
+
+models=list()
+results <- data.frame()
+
+sessionInfo()
+```
+
+
+
+Next we define the utility weights and the utility function. In this case, we simulate data for a conditional logit model, with no observed or unobserved heterogeneity. All attributes but cost are dummy coded and parameter values for the simulation are derived from a previous study. We also define the sample size (i.e. how many replications for the simulation) and create further objects containing the number of blocks, choice sets etc. 
+
+```{r}
+results <- data.frame()
+
+# Setting the beta for utility function 
+basc = -2.5
+bcost = -0.05
+bmanufEU = -0.4
+bmanufnonEU = -1
+bcultEU = -0.2
+bcultnonEU = -0.4
+bCONV = -0.4
+
+
+
+# Simulate reveaved preferences
+respondents <-length(unique(dce_data$RID))
+no_sim <- 100
+
+nsets <- length(unique(dce_data$SCENARIO))
+nblocks <- 1
+#max(dce_data$Block)
+
+setpp <- nsets/nblocks
+
+``` 
+
+We will now create the base dataset. It contains all factors as well as the deterministic utility values of each alternative. This dataset can be used as the basis for the simulation. It does not contain choices, as choices are subject to a random variable. For each run of the simulation, new random variables will be generated.
+
+``` {r}
+
+  database <-  dce_data%>% 
+    #slice(rep(row_number(), 2)) %>%  
+ mutate(RID = rep(1:respondents, each=nsets)) %>% 
+    #group_by(RID) %>% 
+    mutate(across(matches("._x3$") , ~ .x - 1 ) ) %>% 
+    mutate(a1_manuf_EU = a1_x1==2, a1_manuf_nonEU = a1_x1==3 , a2_manuf_EU = a2_x1==2, a2_manuf_nonEU = a2_x1==3 ,
+           a1_cult_EU = a1_x2==2, a1_cult_nonEU = a1_x2==3 , a2_cult_EU = a2_x2==2, a2_cult_nonEU = a2_x2==3 ,
+           across(matches("._x4$") ,  ~ recode(  .x ,`1` = 15,`2` = 20  , `3` = 25 , `4` = 30 , `5` = 35 , `6` = 40 , `7` = 50 , `8` = 60 , .default =99999)),
+                    V.1 =  bcost*a1_x4 + bmanufEU*a1_manuf_EU + bmanufnonEU*a1_manuf_nonEU + bcultEU*a1_cult_EU + bcultnonEU*a1_cult_nonEU + bCONV*a1_x3 , 
+           V.2 = bcost*a2_x4 + bmanufEU*a2_manuf_EU + bmanufnonEU*a2_manuf_nonEU + bcultEU*a2_cult_EU + bcultnonEU*a2_cult_nonEU + bCONV*a2_x3  ,
+           V.3 = basc ,  
+    )   %>% 
+    as.data.frame()
+  
+
+  
+ rm(list = grep(x= ls() , pattern = "^b" , value = TRUE))
+   
+``` 
+
+
+# Inspect the results
+
+
+```{r desc}
+
+
+table(database$DESIGN_ROW)
+  
+  ## For all respondents
+  table(database$a1_x1)
+  table(database$a1_x2)
+  table(database$a1_x3)
+  table(database$a1_x4)
+  table(database$a1_x5)
+  
+  #Only for the first respondents
+  table(database$a1_x1[unique(database$DESIGN_ROW)], useNA = "always")
+  table(database$a1_x2[unique(database$DESIGN_ROW)], useNA = "always")
+  table(database$a1_x3[unique(database$DESIGN_ROW)], useNA = "always")
+  table(database$a1_x4[unique(database$DESIGN_ROW)], useNA = "always")
+  table(database$a1_x5[unique(database$DESIGN_ROW)], useNA = "always")
+  
+
+  ## Over all observations
+  cor(database[, 6:13])
+  
+  ## Only design
+  cor(database[(unique(database$DESIGN_ROW)), 6:13])
+  
+```
+  
+# Simulation
+The next steps are a loop, which replicates the simulation the desired number of times. Replicating the simulation allows us to investigate if our estimates are unbiased. 
+Within the loop we are generating the choices and estimate a conditional logit model with the package apollo. 
+
+
+## Descriptives
+
+Here, we plot some relevant desciptive statistics. The first table shows the frequency of choices for each alternative. Then we create a new dataframe in long format, so that we can plot conditional choice probabilites, i.e. the probability to choose an alternative, given a certain attribute.
+
+```{r}
+
+database<- simulate_choices()
+
+table(database$pref1)
+
+  barplot(prop.table(table(database$pref1)), ylim = c(0,1))
+
+
+long= database %>% rename_with( ~ paste0(.,"_",str_extract(.,"._" )), starts_with("a") )%>% 
+  rename_with( ~ sub("^..." , "", .), starts_with("a")) %>% 
+  pivot_longer(matches("^x"), names_to = c(".value", "set"), names_sep = "_") %>% mutate(set=as.numeric(set), choice=pref1==set, across(starts_with("x"),  ~if_else(set==3,0,.) )  )
+
+
+
+```
+
+```{r}
+for (att in c("x1","x2","x3","x4")) {
+  
+
+
+plot(aggregate(as.numeric(choice)==1 ~ get(att) , data=long , mean), type="o", main=att, ylim=0:1)
+
+}
+
+
+```
+
+
+
+```{r}
+
+for (i in 1:no_sim) {
+  # database <-  database%>% 
+  #   group_by(RID) %>% 
+  #       mutate(
+  #          e.1 = rgumbel(setpp,loc=0, scale=1) ,
+  #          e.2 = rgumbel(setpp,loc=0, scale=1) ,
+  #          e.3 = rgumbel(setpp,loc=0, scale=1) ,
+  #          U.1 = V.1 + e.1 ,
+  #          U.2 = V.2 + e.2 ,
+  #          U.3 = V.3 + e.3 
+  #   )   %>% 
+  #   as.data.frame()
+  # 
+  # # Preferences derived from utility
+  # database$pref1 <- max.col(database[,c("U.1" , "U.2" , "U.3" )])
+
+  database<-simulate_choices()
+  
+apollo_initialise()
+  
+  
+  modelOutput_settings = list(printPVal=T)
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  ="Simulated Data Tofu",
+    modelDescr ="Simple MNL model",
+    indivID    ="RID"
+  )
+  
+  
+  apollo_beta=c(basc = 1.2,
+bcost = 0.2,
+bmanufEU = -0.4,
+bmanufnonEU = -1,
+bcultEU = -0.2,
+bcultnonEU = -0.4,
+bCONV = -0.4)
+  
+  
+  ### keine Parameter fix halten
+  apollo_fixed = c()
+  
+  ### validieren
+  apollo_inputs = apollo_validateInputs()
+  
+  
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### 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()
+    
+    ### List of utilities (later integrated in mnl_settings below)
+    V = list()
+    V[['alt1']] =  bcost*a1_x4 + bmanufEU*a1_manuf_EU + bmanufnonEU*a1_manuf_nonEU + bcultEU*a1_cult_EU + bcultnonEU*a1_cult_nonEU + bCONV*a1_x3 
+    V[['alt2']] = bcost*a2_x4 + bmanufEU*a2_manuf_EU + bmanufnonEU*a2_manuf_nonEU + bcultEU*a2_cult_EU + bcultnonEU*a2_cult_nonEU + bCONV*a2_x3    
+    V[['alt3']] = basc
+
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3) ,
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = pref1,
+      V             = V  # tell function to use list vector defined above
+      
+    )
+    
+    ### 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)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    ### P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  model = apollo_estimate(apollo_beta, apollo_fixed,
+                          apollo_probabilities, apollo_inputs, 
+                          estimate_settings=list(hessianRoutine="maxLik"))
+  
+  apollo_modelOutput(model)
+  
+  models[[i]]=model
+  
+paras <- length(model$estimate)*2
+  
+  results[i,1:paras]<-apollo_modelOutput(model)[,1:2]
+  
+  if (i==1) {
+    names(results) <- c(rownames(apollo_modelOutput(model)) ,paste0("se_",rownames(apollo_modelOutput(model))))
+  }  
+  
+} # end of loop
+```
+
+  
+  
+# Simulation results and Post estimation
+
+The results are stored in the object `results`. Looking at the summary table allows us to inspect the accuracy and biasedness. 
+  
+
+```{r}
+
+#inspect any abritrary model
+modelOutput_settings = list(printPVal=T)
+
+apollo_modelOutput(models[[3]] , modelOutput_settings)
+
+
+describe(results, fast = TRUE)
+
+
+# basc = 1.2
+# bcost = 0.05
+# bmanufEU = -0.4
+# bmanufnonEU = -1
+# bcultEU = -0.2
+# bcultnonEU = -0.4
+# bCONV = -0.4
+```
+  
+
+
+Next, we can estimate WTP values. These are not relevant to inspect the model, but help to interpret the parameters better. We do that only for the last model
+
+
+
+```{r}
+
+wtp("bcost", names(model$estimate), model=model)
+
+
+
+```
+
+
+
+
+

Projects/feedadditives/Designs/Design 10112023.ngd

+Design	Choice situation	alt1.cow	alt1.adv	alt1.vet	alt1.far	alt1.met	alt1.bon	alt2.cow	alt2.adv	alt2.vet	alt2.far	alt2.met	alt2.bon	Block	
+1	1	0	1	1	1	3	6	1	0	1	0	1	0	2	
+1	2	1	0	0	1	0	7	0	1	0	0	1	2	2	
+1	3	1	0	0	1	1	1	1	1	1	0	3	5	2	
+1	4	0	0	1	1	3	0	0	1	0	0	2	3	1	
+1	5	0	0	1	0	2	4	0	1	0	1	1	7	2	
+1	6	1	1	0	1	2	0	0	0	1	0	1	5	2	
+1	7	1	0	1	0	2	3	0	1	0	1	0	4	1	
+1	8	0	0	0	1	0	5	0	1	1	0	2	7	1	
+1	9	1	1	0	0	0	2	0	0	0	1	2	0	1	
+1	10	1	1	1	0	1	4	0	0	0	1	0	2	2	
+1	11	0	1	1	0	0	1	1	0	0	1	3	4	2	
+1	12	1	1	0	0	1	6	1	0	1	1	3	1	1	
+1	13	0	0	1	1	1	2	1	1	0	0	3	1	1	
+1	14	1	1	0	1	2	5	1	0	1	0	0	6	1	
+1	15	0	0	0	0	3	7	1	1	1	1	0	3	2	
+1	16	0	1	1	0	3	3	1	0	1	1	2	6	1	
+||||||||||
+design  
+ ;alts = alt1*, alt2*, alt3  
+ ;eff = (mnl, d)  
+ ;alg = swap  
+ ;rows = 16  
+ ;block = 2  
+ ;model:  
+ U(alt1) = b1[0.3] * COW[0,1]    
+         + b2[0.3] * ADV[0,1]    
+         + b3[0.3] * VET[0,1]    
+         + b4[0.3] * FAR[0,1]    
+         + b5.dummy[0.3|0.3|0.3] * MET[1,2,3,0]      
+         + b6.dummy[0.6|0.3|0.3|0.3|0.3|0.3|0.3] * BON[1,2,3,4,5,6,7,0]           
+         + i1[0] * COW.dummy[0] * VET.dummy[1]  
+ /  
+ U(alt2) = b1 * COW  
+         + b2 * ADV  
+         + b3 * VET  
+         + b4 * FAR  
+         + b5 * MET      
+         + b6 * BON     
+         + i1 * COW.dummy[0] * VET.dummy[1]  
+/  
+ U(alt3) = asc3[0]  
+;
+$
\ No newline at end of file

Projects/feedadditives/parameters_feedadd.R

+
+
+designpath<- "Projects/feedadditives/Designs/"
+
+resps =360  # number of respondents
+nosim=500 # number of simulations to run (about 500 is minimum)
+
+#betacoefficients should not include "-"
+basc = 0.2
+bcow = 0.3
+badv = 0.3
+bvet = 0.3
+bfar = 0.3
+bmet = 0.3
+bbon = 0.3
+
+
+
+#place your utility functions here
+u<- list(u1= 
+  list(
+  v1 =V.1 ~  bcow*alt1.cow + badv * alt1.adv + bvet * alt1.vet + bfar * alt1.far + bmet*alt1.met + bbon * alt1.bon,
+  v2 =V.2 ~  bcow*alt2.cow + badv * alt2.adv + bvet * alt2.vet + bfar * alt2.far + bmet*alt2.met + bbon * alt2.bon,
+  v3 =V.3 ~ basc)
+)

generatemd.R

@@ -2,13 +2,13 @@
 rm(list=ls())
 
 #file <- "Projects/ValuGaps/parameters_valugaps.R"
-  file <- "Projects/SE_AGRI/parameters_SE Design-Agri.R"
+  file <- "Projects/feedadditives/parameters_feedadd.R"
 
 
 rmarkdown::render("simulation_output.rmd",
                   output_file = paste0(
                     stringr::str_remove_all(
-                      file,"parameters_|.R$"),".html"),
+                      file,"parameters_|.R$"),"_noheur.html"),
                   params = list(file=file)
                   )