Skip to content
Snippets Groups Projects
Commit 51c31bb8 authored by dj44vuri's avatar dj44vuri
Browse files

simulate choices and read_design now outside simchoice

parent e8c61dc3
No related branches found
No related tags found
No related merge requests found
rm(list=ls())
source("functions.R")
designpath<- "Projects/Rbook/Designs/"
#notes <- "This design consists of different heuristics. One group did not attend the methan attribute, another group only decided based on the payment"
......@@ -23,13 +25,6 @@ bheight3=0.50
destype <- "spdesign"
#place your utility functions here
u<- list(u1= list(
v1 =V.1 ~ bsq * alt1.sq,
......@@ -40,7 +35,7 @@ u<- list(u1= list(
## logBonus
sim_all()
......
......@@ -78,7 +78,113 @@ all_designs[["time"]]=time
return(all_designs)
}
##test
readdesign <- function(design = designfile, designtype=destype) {
design <- switch(designtype,
"ngene" =readr::read_delim(design,delim = "\t",
escape_double = FALSE,
trim_ws = TRUE ,
col_select = c(-Design, -starts_with("...")) ,
name_repair = "universal" , show_col_types = FALSE) %>%
filter(!is.na(Choice.situation)),
"spdesign" =as.data.frame(readRDS(design)) %>%
mutate(Choice.situation=1:n()) %>%
rename_with(~ stringr::str_replace(.,pattern = "_","\\."), everything()),
stop("Invalid value for design. Please provide either 'ngene' or 'spdesign'.")
)
}
simulate_choices <- function(data=datadet, utility =utils, setspp) { #the part in dataset that needs to be repeated in each run
by_formula <- function(equation){ #used to take formulas as inputs in simulation utility function
# //! cur_data_all may get deprecated in favor of pick
# pick(everything()) %>%
cur_data_all() %>%
transmute(!!lhs(equation) := !!rhs(equation) )
}
cat(" \n does sou_gis exist: ", exists("sou_gis"), "\n")
if (exists("sou_gis") && is.function(sou_gis)) {
sou_gis()
cat("\n source of gis has been done \n")
}
# source("Projects/ValuGaps/code/GIS_data.R")
if(!exists("manipulations")) manipulations=list() ## If no user input on further data manipulations
n=seq_along(1:length(utility)) # number of utility functions
cat("\n dataset final_set exists: ",exists("final_set"), "\n")
if(exists("final_set")) data = left_join(data,final_set, by = "ID")
cat("\n decisiongroups exists: " ,exists("decisiongroups"))
if(exists("decisiongroups")) { ### create a new variable to classify decision groups.
data = mutate(data,group = as.numeric(cut(row_number(),
breaks = decisiongroups * n(),
labels = seq_along(decisiongroups[-length(decisiongroups)]),
include.lowest = TRUE)))
print(table(data$group))
} else {
data$group=1
}
data<- data %>%
group_by(ID) %>%
mutate(!!! manipulations)
subsets<- split(data,data$group)
subsets <- map2(.x = seq_along(u),.y = subsets,
~ mutate(.y,map_dfc(u[[.x]],by_formula)))
data <-bind_rows(subsets)
data<- data %>%
rename_with(~ stringr::str_replace(.,pattern = "\\.","_"), everything()) %>%
mutate(across(.cols=n,.fns = ~ rgumbel(setspp,loc=0, scale=1), .names = "{'e'}_{n}" ),
across(starts_with("V_"), .names = "{'U'}_{n}") + across(starts_with("e_")) ) %>% ungroup() %>%
mutate(CHOICE=max.col(.[,grep("U_",names(.))])
) %>%
as.data.frame()
cat("\n data has been made \n")
cat("\n First few observations \n ")
print(head(data))
cat("\n \n ")
return(data)
}
sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1]] ,destype) {
......@@ -101,34 +207,13 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1]
require("rlang")
readdesign <- function(design = designfile, designtype=destype) {
design <- switch(designtype,
"ngene" =readr::read_delim(design,delim = "\t",
escape_double = FALSE,
trim_ws = TRUE ,
col_select = c(-Design, -starts_with("...")) ,
name_repair = "universal" , show_col_types = FALSE) %>%
filter(!is.na(Choice.situation)),
"spdesign" =as.data.frame(readRDS(design)) %>%
mutate(Choice.situation=1:n()) %>%
rename_with(~ stringr::str_replace(.,pattern = "_","\\."), everything()),
stop("Invalid value for design. Please provide either 'ngene' or 'spdesign'.")
)
}
mnl_U <-paste(map_chr(utils,as.character,keep.source.attr = TRUE),collapse = "",";") %>%
str_replace_all( c( "priors\\[\"" = "" , "\"\\]" = "" , "~" = "=", "\\." = "_" , " b" = " @b" , "V_"="U_", " alt"="$alt"))
cat("mixl \n")
print(mnl_U)
cat(mnl_U)
cat("\n Simulation \n")
......@@ -136,93 +221,19 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1]
by_formula <- function(equation){ #used to take formulas as inputs in simulation utility function
# //! cur_data_all may get deprecated in favor of pick
# pick(everything()) %>%
cur_data_all() %>%
transmute(!!lhs(equation) := !!rhs(equation) )
}
simulate_choices <- function(data=datadet) { #the part in dataset that needs to be repeated in each run
cat(" \n does sou_gis exist: ", exists("sou_gis"), "\n")
if (exists("sou_gis") && is.function(sou_gis)) {
sou_gis()
cat("\n source of gis has been done \n")
}
# source("Projects/ValuGaps/code/GIS_data.R")
if(!exists("manipulations")) manipulations=list() ## If no user input on further data manipulations
n=seq_along(1:length(utils)) # number of utility functions
#browser()
cat("\n dataset final_set exists: ",exists("final_set"), "\n")
if(exists("final_set")) data = left_join(data,final_set, by = "ID")
cat("\n decisiongroups exists: " ,exists("decisiongroups"))
if(exists("decisiongroups")) { ### create a new variable to classify decision groups.
data = mutate(data,group = as.numeric(cut(row_number(),
breaks = decisiongroups * n(),
labels = seq_along(decisiongroups[-length(decisiongroups)]),
include.lowest = TRUE)))
print(table(data$group))
} else {
data$group=1
}
data<- data %>%
group_by(ID) %>%
mutate(!!! manipulations)
subsets<- split(data,data$group)
subsets <- map2(.x = seq_along(u),.y = subsets,
~ mutate(.y,map_dfc(u[[.x]],by_formula)))
data <-bind_rows(subsets)
data<- data %>%
rename_with(~ stringr::str_replace(.,pattern = "\\.","_"), everything()) %>%
mutate(across(.cols=n,.fns = ~ rgumbel(setpp,loc=0, scale=1), .names = "{'e'}_{n}" ),
across(starts_with("V_"), .names = "{'U'}_{n}") + across(starts_with("e_")) ) %>% ungroup() %>%
mutate(CHOICE=max.col(.[,grep("U_",names(.))])
) %>%
as.data.frame()
print("\n data has been made \n")
cat("\n First few observations \n ")
print(head(data))
cat("\n \n ")
return(data)
}
estimate_sim <- function(run=1) { #start loop
cat("This is Run number ", run)
database <- simulate_choices(datadet)
database <- simulate_choices(datadet, utility = utils, setspp=setpp )
cat("This is the utility functions \n" , mnl_U)
......@@ -236,15 +247,15 @@ print(head(data))
designs_all <- list()
design<- readdesign()
design<- readdesign(design = designfile)
if (!exists("design$Block")) design$Block=1
nsets<-nrow(design)
nblocks<-max(design$Block)
setpp <- nsets/nblocks # Choice Sets per respondent; in this 'no blocks' design everyone sees all 24 sets
#respondents <- replications*nblocks
replications <- respondents/nblocks
#browser()
......@@ -259,24 +270,23 @@ as.data.frame()
database <- simulate_choices()
database <- simulate_choices(data=datadet, utility = utils, setspp = setpp)
model_spec <- mixl::specify_model(mnl_U, database, disable_multicore=F)
model_spec <- mixl::specify_model(mnl_U, database, disable_multicore=F)
est=setNames(rep(0,length(model_spec$beta_names)), model_spec$beta_names)
est=setNames(rep(0,length(model_spec$beta_names)), model_spec$beta_names)
availabilities <- mixl::generate_default_availabilities(
availabilities <- mixl::generate_default_availabilities(
database, model_spec$num_utility_functions)
plan(multisession, workers = 8)
output<- 1:no_sim %>% map(estimate_sim)
output<- 1:no_sim %>% map(estimate_sim)
coefs<-map(1:length(output),~summary(output[[.]])[["coefTable"]][c(1,8)] %>%
coefs<-map(1:length(output),~summary(output[[.]])[["coefTable"]][c(1,8)] %>%
tibble::rownames_to_column() %>%
pivot_wider(names_from = rowname, values_from = c(est, rob_pval0)) ) %>%
bind_rows(.id = "run")
......
......@@ -3,9 +3,9 @@ rm(list=ls())
source("functions.R")
file <- "Projects/feedadditives/parameters_feedadd.R"
file <- "Projects/Rbook/parameters_Rbook.R"
# file <- "Projects/CSA/parameters_csa.R"
make_md()
\ No newline at end of file
make_md()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment