diff --git a/Projects/Tofu/Designs/tofu_main_swe_tofu.xlsx b/Projects/Tofu/Designs/tofu_main_swe_tofu.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..0a765f9df0df22fd2797425d8f1934dbcf03bd21
Binary files /dev/null and b/Projects/Tofu/Designs/tofu_main_swe_tofu.xlsx differ
diff --git a/Projects/Tofu/parameters_tofu.R b/Projects/Tofu/parameters_tofu.R
new file mode 100644
index 0000000000000000000000000000000000000000..e5a0bbf84aed9f9a852f061673185830e5e4f927
--- /dev/null
+++ b/Projects/Tofu/parameters_tofu.R
@@ -0,0 +1,39 @@
+
+
+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)
+)
diff --git a/Projects/Tofu/simulation_SWE_after_pilot.Rmd b/Projects/Tofu/simulation_SWE_after_pilot.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..bc2f390eb9cd748c0ec48759c407ffa5c62f2852
--- /dev/null
+++ b/Projects/Tofu/simulation_SWE_after_pilot.Rmd
@@ -0,0 +1,332 @@
+---
+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)
+
+
+
+```
+
+
+
+
+
diff --git a/Projects/feedadditives/Designs/Design 10112023.ngd b/Projects/feedadditives/Designs/Design 10112023.ngd
new file mode 100644
index 0000000000000000000000000000000000000000..d60e4d6a9352c77c04a02ff6b280fa7c45fdc585
--- /dev/null
+++ b/Projects/feedadditives/Designs/Design 10112023.ngd
@@ -0,0 +1,44 @@
+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
diff --git a/Projects/feedadditives/parameters_feedadd.R b/Projects/feedadditives/parameters_feedadd.R
new file mode 100644
index 0000000000000000000000000000000000000000..33b0ae5e4f29d440ebb4e87a0de5078818811276
--- /dev/null
+++ b/Projects/feedadditives/parameters_feedadd.R
@@ -0,0 +1,25 @@
+
+
+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)
+)
diff --git a/generatemd.R b/generatemd.R
index 5808e8bc81acb2d410945c73d35b034f173dc641..a03acb82feb97b3e762ebe66f65540cc39eac4ea 100644
--- a/generatemd.R
+++ b/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)
)