From 3eeb6820cd7bf13a18930671438679ab74a26e69 Mon Sep 17 00:00:00 2001 From: Julian Sagebiel <julian.sagebiel@idiv.de> Date: Fri, 21 Jun 2024 23:49:29 +0200 Subject: [PATCH] added tests for decisiongroups and latent class model --- inst/extdata/spdesigns/designs/twoattr.RDS | Bin 0 -> 631 bytes inst/extdata/spdesigns/twoattr.R | 26 ++ tests/manual-tests/test_decisiongroups.R | 273 +++++++++++++++++++++ 3 files changed, 299 insertions(+) create mode 100644 inst/extdata/spdesigns/designs/twoattr.RDS create mode 100644 inst/extdata/spdesigns/twoattr.R create mode 100644 tests/manual-tests/test_decisiongroups.R diff --git a/inst/extdata/spdesigns/designs/twoattr.RDS b/inst/extdata/spdesigns/designs/twoattr.RDS new file mode 100644 index 0000000000000000000000000000000000000000..e43085bd732b5cf4297e08b3593992b6dfb325bf GIT binary patch literal 631 zcmV--0*L(|iwFP!000001C>?Hi_<_9pC(hAT|X$a;K570h`ZX_1O@SuiYVwo1r-+7 z64E4H0{K`b(`v<oR|O9u2nrqq{SWluK}GPW;?c#6H$5oq(bLA6WF}3<joN|fo8SBR z&CHv;*Of&8fD9n40l7|=pnm1@na**dYZ1@@K%!D2vq@$RK(IjMk)e5ATkUnVfvPU6 zWT~AvHc$_#M6DC`(E-)!tM<X(VO3iisN03B*`VRsDhtoGY()nEP%4DiLo0A&#wjr+ zw0x`;);f_}v*XKqD8b2%$2Z3BZ}&t(+qYtouyOIirL!wGW^hLeamG4hyZLkf*Nja{ z<{+Q}H8Efjdx2}l*g`x;Ni}ZhWKpFg%0cK8(MnZ7q9>A!mk;T9WNc)zVk1cNUs$tT zrB4pYq+IR30@BfMGudQPDNAEhus6$bxg7|&El8#8m-B$K<njERK+ua)9P|3ktKY_+ z<+GWAK2!dYKQSMozYr(-<vbPdBEEyj!_PVIy{qh1`S@R4ozGIW$-KPq&^MhCFAjou zRXG;6x+6q-hEVnkD7x0=y<|6;(k=3nB*_n<ZXmB{+GA#=!=PS(dfBSywUyQ<G(H5c z9z65!C%@R;UO)YK;9g^2f9cz~&i8kBD);TU2t|lPy9}K&j2SHM`{WdSyLo=^o9WKG z=I7<-OujI_oW8oXHvA+cMs~b_r(`xuE(&qssVn|N{H?|O%P#iC2EWl5d%lOK+?_T> z-noOwarsN{T1Q@l%!%cX$+_5P;kdE47BbXwM<dVn+|ZtyHuA8GJd3kV>}Vo<Tlr<< Rbw)!D`5P*~mU%J<001#RL<|4` literal 0 HcmV?d00001 diff --git a/inst/extdata/spdesigns/twoattr.R b/inst/extdata/spdesigns/twoattr.R new file mode 100644 index 0000000..ea90f0f --- /dev/null +++ b/inst/extdata/spdesigns/twoattr.R @@ -0,0 +1,26 @@ + + +# Taken from Erlend Sandorf +# Example file for creating a simple MNL design +# +rm(list = ls(all = TRUE)) +library(spdesign) + +# Define the list of utility functions ---- +#' Specifying a utility function with 3 attributes and a constant for the +#' SQ alternative. The design has 20 rows. +utility <- list( + alt1 = "b_x1[0.1] * x1[1:5] + b_x2[0.4] * x2[c(0, 1)] ", + alt2 = "b_x1 * x1 + b_x2 * x2 " +) + +# Generate designs ---- +design <- generate_design(utility, rows = 20, + model = "mnl", efficiency_criteria = "d-error", + algorithm = "rsc", draws = "scrambled-sobol", + control = list( + max_iter = 2000, + max_no_improve = 5000 + )) + +saveRDS(design,"inst/extdata/spdesigns/designs/twoattr.RDS") diff --git a/tests/manual-tests/test_decisiongroups.R b/tests/manual-tests/test_decisiongroups.R new file mode 100644 index 0000000..76c3511 --- /dev/null +++ b/tests/manual-tests/test_decisiongroups.R @@ -0,0 +1,273 @@ + +rm(list=ls()) +devtools::load_all() + +library(dplyr) +designpath<- system.file("extdata","spdesigns","designs", package = "simulateDCE") + +destype = 'spdesign' +resps =4000 # number of respondents +nosim=2 # number of simulations to run (about 500 is minimum) + +#betacoefficients should not include "-" +bcoeff <- list( + bx1A = 0.1, ## very high asc + bx1B = -0.1, + bx2A = 0.2, + bx2B = -0.2 + ) + + +desisiongroups =c(0,0.3,1) + +ul<- list(uA= + list( + v1 = V.1 ~ bx1A*alt1.x1+ bx2A * alt1.x2, + v2 = V.2 ~ bx1A*alt2.x1+ bx2A * alt2.x2), + uB = list( + v1 = V.1 ~ bx1B*alt1.x1+ bx2B * alt1.x2, + v2 = V.2 ~ bx1B*alt2.x1+ bx2B * alt2.x2) + +) + +simplesim <- sim_all(nosim = nosim, resps = resps,designpath = designpath, bcoeff = bcoeff, u = ul,destype = "spdesign" , decisiongroups = desisiongroups) + +testdata = simplesim[["twoattr"]][[1]][["data"]] %>% + distinct(group, Choice_situation ,.keep_all = T) %>% + arrange(Choice_situation) + +## estimate models + +library(apollo) + +apollo_initialise() + + +apollo_control = list( + modelName = "Clogit_simpledesign", + modelDescr = "Simple conditional logit model for group 1 ", + indivID = "ID", + nCores = 1, + outputDirectory = "/~" +) + + +## group 1 + +database <- simplesim[["twoattr"]][[1]][["data"]] +apollo_beta = c( + + b_x1 = 0, + b_x2 = 0 +) + +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() + + ### Define settings for MNL model component that are generic across classes + mnl_settings = list( + alternatives = c(alt1=1, alt2=2), + avail = list(alt1=1, alt2=1), + choiceVar = CHOICE + ) + + + + + ### Compute class-specific utilities + V=list() + V[["alt1"]] = b_x1*alt1_x1 + b_x2*alt1_x2 + V[["alt2"]] = b_x1*alt2_x1 + b_x2*alt2_x2 + + mnl_settings$utilities = V + #mnl_settings$componentName = paste0("Class_",s) + + + + ### 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 #### +# ################################################################# # + +### Optional starting values search +# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs) + +model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs ,estimate_settings =list(estimationRoutine="BFGS" , hessianRoutine= "numDeriv")) + + + +apollo_modelOutput(model) + + + + +### latent class model + +apollo_control = list( + modelName = "Latent Class logit_simpledesign", + modelDescr = "Latent Class logit model", + indivID = "ID", + nCores = 1, + outputDirectory = "/~" +) + + +## group 1 + +database <- simplesim[["twoattr"]][[1]][["data"]] + +apollo_beta = c( + b_x1_a = 0, + b_x2_a = 0, + b_x1_b = 0, + b_x2_b = 0, + delta_a = 0, + delta_b = 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("delta_b") + +# ################################################################# # +#### DEFINE LATENT CLASS COMPONENTS #### +# ################################################################# # + +apollo_lcPars=function(apollo_beta, apollo_inputs){ + lcpars = list() + + lcpars[["b_x1"]] = list(b_x1_a, b_x1_b) + lcpars[["b_x2"]] = list(b_x2_a, b_x2_b) + + V=list() + V[["class_a"]] = delta_a + V[["class_b"]] = delta_b + + + classAlloc_settings = list( + classes = c(class_a = 1, class_b = 2), + utilities = V + ) + + lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings) + + return(lcpars) +} + + + +# ################################################################# # +#### GROUP AND VALIDATE INPUTS #### +# ################################################################# # + +apollo_inputs = apollo_validateInputs() + +# ################################################################# # +#### DEFINE MODEL AND LIKELIHOOD FUNCTION #### +# ################################################################# # + + +# ################################################################# # +#### 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() + + ### Define settings for MNL model component that are generic across classes + mnl_settings = list( + alternatives = c(alt1=1, alt2=2), + avail = list(alt1=1, alt2=1), + choiceVar = CHOICE + ) + + ### Loop over classes + for(s in 1:2){ + + ### Compute class-specific utilities + V=list() + V=list() + V[["alt1"]] = b_x1[[s]]*alt1_x1 + b_x2[[s]]*alt1_x2 + V[["alt2"]] = b_x1[[s]]*alt2_x1 + b_x2[[s]]*alt2_x2 + + + + mnl_settings$utilities = V + #mnl_settings$componentName = paste0("Class_",s) + + ### Compute within-class choice probabilities using MNL model + P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality) + + ### Take product across observation for same individual + P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs ,functionality) + } + + ### Compute latent class model probabilities + lc_settings = list(inClassProb = P, classProb=pi_values) + P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality) + + ### Prepare and return outputs of function + P = apollo_prepareProb(P, apollo_inputs, functionality) + return(P) +} + +# ################################################################# # +#### MODEL ESTIMATION #### +# ################################################################# # + +### Optional starting values search +# apollo_beta=apollo_searchStart(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs) + +model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs ,estimate_settings =list(estimationRoutine="bhhh" , hessianRoutine= "maxLik")) + +# ################################################################# # +#### MODEL OUTPUTS #### +# ################################################################# # + +# ----------------------------------------------------------------- # +#---- FORMATTED OUTPUT (TO SCREEN) ---- +# ----------------------------------------------------------------- # + +apollo_modelOutput(model) + + -- GitLab