Skip to content
Snippets Groups Projects
Commit 43f6772f authored by dj44vuri's avatar dj44vuri
Browse files

small change in code structure

parent 51c31bb8
Branches
No related tags found
No related merge requests found
...@@ -123,16 +123,12 @@ simulate_choices <- function(data=datadet, utility =utils, setspp) { #the part ...@@ -123,16 +123,12 @@ simulate_choices <- function(data=datadet, utility =utils, setspp) { #the part
cat("\n source of gis has been done \n") 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 if(!exists("manipulations")) manipulations=list() ## If no user input on further data manipulations
n=seq_along(1:length(utility)) # number of utility functions n=seq_along(1:length(utility)) # number of utility functions
cat("\n dataset final_set exists: ",exists("final_set"), "\n") cat("\n dataset final_set exists: ",exists("final_set"), "\n")
if(exists("final_set")) data = left_join(data,final_set, by = "ID") if(exists("final_set")) data = left_join(data,final_set, by = "ID")
...@@ -207,28 +203,6 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1] ...@@ -207,28 +203,6 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1]
require("rlang") require("rlang")
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")
cat(mnl_U)
cat("\n Simulation \n")
print(u)
estimate_sim <- function(run=1) { #start loop estimate_sim <- function(run=1) { #start loop
cat("This is Run number ", run) cat("This is Run number ", run)
...@@ -236,40 +210,45 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1] ...@@ -236,40 +210,45 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u[[1]
database <- simulate_choices(datadet, utility = utils, setspp=setpp ) database <- simulate_choices(datadet, utility = utils, setspp=setpp )
cat("This is the utility functions \n" , mnl_U) cat("This is the utility functions \n" , mnl_U)
model<-mixl::estimate(model_spec,start_values = est, availabilities = availabilities, data= database,) model<-mixl::estimate(model_spec,start_values = est, availabilities = availabilities, data= database,)
return(model) return(model)
} }
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")
cat(mnl_U)
cat("\n Simulation \n")
print(u)
designs_all <- list() designs_all <- list()
design<- readdesign(design = designfile) design<- readdesign(design = designfile)
if (!exists("design$Block")) design$Block=1 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
replications <- respondents/nblocks nsets<-nrow(design)
nblocks<-max(design$Block)
setpp <- nsets/nblocks # Choice Sets per respondent; in this 'no blocks' design everyone sees all 24 sets
replications <- respondents/nblocks
datadet<- design %>%
datadet<- design %>%
arrange(Block,Choice.situation) %>% arrange(Block,Choice.situation) %>%
slice(rep(row_number(), replications)) %>% ## replicate design according to number of replications slice(rep(row_number(), replications)) %>% ## replicate design according to number of replications
mutate(ID = rep(1:respondents, each=setpp)) %>% # create Respondent ID. mutate(ID = rep(1:respondents, each=setpp)) %>% # create Respondent ID.
relocate(ID,`Choice.situation`) %>% relocate(ID,`Choice.situation`) %>%
as.data.frame() as.data.frame()
database <- simulate_choices(data=datadet, utility = utils, setspp = setpp) 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)
...@@ -291,22 +270,22 @@ coefs<-map(1:length(output),~summary(output[[.]])[["coefTable"]][c(1,8)] %>% ...@@ -291,22 +270,22 @@ coefs<-map(1:length(output),~summary(output[[.]])[["coefTable"]][c(1,8)] %>%
pivot_wider(names_from = rowname, values_from = c(est, rob_pval0)) ) %>% pivot_wider(names_from = rowname, values_from = c(est, rob_pval0)) ) %>%
bind_rows(.id = "run") bind_rows(.id = "run")
output[["summary"]] <-psych::describe(coefs[,-1], fast = TRUE) output[["summary"]] <-psych::describe(coefs[,-1], fast = TRUE)
output[["coefs"]] <-coefs output[["coefs"]] <-coefs
pvals <- output[["coefs"]] %>% dplyr::select(starts_with("rob_pval0")) pvals <- output[["coefs"]] %>% dplyr::select(starts_with("rob_pval0"))
output[["power"]] <- 100*table(apply(pvals,1, function(x) all(x<0.05)))/nrow(pvals) output[["power"]] <- 100*table(apply(pvals,1, function(x) all(x<0.05)))/nrow(pvals)
output[["metainfo"]] <- c(Path = designfile, NoSim = no_sim, NoResp =respondents) output[["metainfo"]] <- c(Path = designfile, NoSim = no_sim, NoResp =respondents)
print(kable(output[["summary"]],digits = 2, format = "rst")) print(kable(output[["summary"]],digits = 2, format = "rst"))
print(output[["power"]]) print(output[["power"]])
return(output) return(output)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment