Code: Select all
library(apollo)
apollo_initialise()
apollo_control = list(
modelName = "Apollo_Jaeren_MXL" ,
modelDescr = "MNL model of repeated data" ,
indivID = "RespID",
mixing=TRUE ,
nCores = 3 #,
# calculateLLC=FALSE
)
apollo_beta=c(beta.TC = -4.399 ,
beta.Dist = 0.663 ,
beta.asc1 = 0 ,
beta.asc2 = 0.205198 ,
beta.asc3 = 0.349196 ,
beta.asc4 = -1.792595,
beta.asc5 = -1.557427,
beta.asc6 = 0.003143 ,
beta.asc7 = -1.303693,
beta.asc8 = -0.176013,
beta.asc9 = -2.028402,
beta.asc10 = -1.228613,
beta.asc11 = -1.492752,
beta.asc12 = 0.61742 ,
beta.asc13 = 1.087001 ,
beta.asc14 = -3.774711,
beta.asc15 = 0.245711 ,
beta.asc16 = -1.825746 ,
beta.asc17 = -0.291979,
beta.asc18 = -1.343846,
beta.asc19 = -2.573377 ,
beta.asc20 = -1.837075,
beta.asc21 = -0.725732,
beta.asc22 = -2.518054,
beta.asc23 = -6.16141,
beta.asc24 = -4.09383,
beta.asc25 = -2.845921,
beta.asc26 = -2.935369,
beta.asc27 = -6.743 ,
beta.asc28 = -2.495 ,
beta.asc29 = -4.805 ,
sigma.TC = 1.323 ,
sigma.Dist = 0.517 ,
mu_RP = 1.0000 ,
mu_SP = 0.674 ,
beta.TC.SP = 0 , # 0.01104 ,
beta.SP.asc1 = 0 ,
beta.SP.asc2 = 0 ,
beta.SP.asc3 = 0 ,
beta.SP.asc4 = 0 ,
beta.SP.asc5 = 0 ,
beta.SP.asc6 = 0 ,
beta.SP.asc7 = 0 ,
beta.SP.asc8 = 0 ,
beta.SP.asc9 = 0 ,
beta.SP.asc10 = 0 ,
beta.SP.asc11 = 0 ,
beta.SP.asc12 = 0 ,
beta.SP.asc13 = 0 ,
beta.SP.asc14 = 0 ,
beta.SP.asc15 = 0 ,
beta.SP.asc16 = 0 ,
beta.SP.asc17 = 0 ,
beta.SP.asc18 = 0 ,
beta.SP.asc19 = 0 ,
beta.SP.asc20 = 0 ,
beta.SP.asc21 = 0 ,
beta.SP.asc22 = 0 ,
beta.SP.asc23 = 0 ,
beta.SP.asc24 = 0 ,
beta.SP.asc25 = 0 ,
beta.SP.asc26 = 0 ,
beta.SP.asc27 = 0 ,
beta.SP.asc28 = 0 ,
beta.SP.asc29 = 0 ,
sigma.Dist.TC = 0
)
apollo_fixed = c("beta.asc1", "mu_RP" , #"beta.TC.SP" ,
"beta.asc2",
"beta.asc3",
"beta.asc4",
"beta.asc5",
"beta.asc6",
"beta.asc7",
"beta.asc8",
"beta.asc9",
"beta.asc10",
"beta.asc11",
"beta.asc12",
"beta.asc13",
"beta.asc14",
"beta.asc15",
"beta.asc16",
"beta.asc17",
"beta.asc18",
"beta.asc19",
"beta.asc20",
"beta.asc21",
"beta.asc22",
"beta.asc23",
"beta.asc24",
"beta.asc25",
"beta.asc26" ,
"beta.SP.asc1" ,
"beta.SP.asc2" ,
"beta.SP.asc3" ,
"beta.SP.asc4" ,
"beta.SP.asc5" ,
"beta.SP.asc6" ,
"beta.SP.asc7" ,
"beta.SP.asc8" ,
"beta.SP.asc9" ,
"beta.SP.asc10" ,
"beta.SP.asc11" ,
"beta.SP.asc12" ,
"beta.SP.asc13" ,
"beta.SP.asc14" ,
"beta.SP.asc15" ,
"beta.SP.asc16" ,
"beta.SP.asc17" ,
"beta.SP.asc18" ,
"beta.SP.asc19" ,
"beta.SP.asc20" ,
"beta.SP.asc21" ,
"beta.SP.asc22" ,
"beta.SP.asc23" ,
"beta.SP.asc24" ,
"beta.SP.asc25" ,
"beta.SP.asc26" ,
"beta.SP.asc27" ,
"beta.SP.asc28" ,
"beta.SP.asc29"
)
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interUnifDraws = c(), #"draws2a", "draws2b"
interNormDraws = c("draws_tc_inter",
"draws_Diste_inter"
),
intraDrawsType = "halton",
intraNDraws = 10,
intraUnifDraws = c(),
intraNormDraws = c()
)
## Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["TC"]] = -exp(beta.TC + sigma.TC * draws_tc_inter) # beta.TC + sigma.TC*(draws2a + draws2b) #beta.TC + sigma.TC * draws_tc_inter #
randcoeff[["Dist"]] = beta.Dist + sigma.Dist * draws_Diste_inter + sigma.Dist.TC * draws_tc_inter
return(randcoeff)
}
#apollo_beta=apollo_readBeta(apollo_beta,apollo_fixed,"Repeated_RP_27_R.csv",overwriteFixed=FALSE)
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: these must use the same names as in mnl_settings, order is irrelevant
V = list()
V[['B1']] = TC*TC1 + beta.asc1 + Dist*DIST1 +beta.TC.SP*TC1 *SP + beta.SP.asc1 * SP
V[['B2']] = TC*TC2 + beta.asc2 + Dist*DIST2 +beta.TC.SP*TC2 *SP + beta.SP.asc2 * SP
V[['B3']] = TC*TC3 + beta.asc3 + Dist*DIST3 +beta.TC.SP*TC3 *SP + beta.SP.asc3 * SP
V[['B4']] = TC*TC4 + beta.asc4 + Dist*DIST4 +beta.TC.SP*TC4 *SP + beta.SP.asc4 * SP
V[['B5']] = TC*TC5 + beta.asc5 + Dist*DIST5 +beta.TC.SP*TC5 *SP + beta.SP.asc5 * SP
V[['B6']] = TC*TC6 + beta.asc6 + Dist*DIST6 +beta.TC.SP*TC6 *SP + beta.SP.asc6 * SP
V[['B7']] = TC*TC7 + beta.asc7 + Dist*DIST7 +beta.TC.SP*TC7 *SP + beta.SP.asc7 * SP
V[['B8']] = TC*TC8 + beta.asc8 + Dist*DIST8 +beta.TC.SP*TC8 *SP + beta.SP.asc8 * SP
V[['B9']] = TC*TC9 + beta.asc9 + Dist*DIST9 +beta.TC.SP*TC9 *SP + beta.SP.asc9 * SP
V[['B10']] = TC*TC10 + beta.asc10 + Dist*DIST10 +beta.TC.SP*TC10*SP + beta.SP.asc10 * SP
V[['B11']] = TC*TC11 + beta.asc11 + Dist*DIST11 +beta.TC.SP*TC11*SP + beta.SP.asc11 * SP
V[['B12']] = TC*TC12 + beta.asc12 + Dist*DIST12 +beta.TC.SP*TC12*SP + beta.SP.asc12 * SP
V[['B13']] = TC*TC13 + beta.asc13 + Dist*DIST13 +beta.TC.SP*TC13*SP + beta.SP.asc13 * SP
V[['B14']] = TC*TC14 + beta.asc14 + Dist*DIST14 +beta.TC.SP*TC14*SP + beta.SP.asc14 * SP
V[['B15']] = TC*TC15 + beta.asc15 + Dist*DIST15 +beta.TC.SP*TC15*SP + beta.SP.asc15 * SP
V[['B16']] = TC*TC16 + beta.asc16 + Dist*DIST16 +beta.TC.SP*TC16*SP + beta.SP.asc16 * SP
V[['B17']] = TC*TC17 + beta.asc17 + Dist*DIST17 +beta.TC.SP*TC17*SP + beta.SP.asc17 * SP
V[['B18']] = TC*TC18 + beta.asc18 + Dist*DIST18 +beta.TC.SP*TC18*SP + beta.SP.asc18 * SP
V[['B19']] = TC*TC19 + beta.asc19 + Dist*DIST19 +beta.TC.SP*TC19*SP + beta.SP.asc19 * SP
V[['B20']] = TC*TC20 + beta.asc20 + Dist*DIST20 +beta.TC.SP*TC20*SP + beta.SP.asc20 * SP
V[['B21']] = TC*TC21 + beta.asc21 + Dist*DIST21 +beta.TC.SP*TC21*SP + beta.SP.asc21 * SP
V[['B22']] = TC*TC22 + beta.asc22 + Dist*DIST22 +beta.TC.SP*TC22*SP + beta.SP.asc22 * SP
V[['B23']] = TC*TC23 + beta.asc23 + Dist*DIST23 +beta.TC.SP*TC23*SP + beta.SP.asc23 * SP
V[['B24']] = TC*TC24 + beta.asc24 + Dist*DIST24 +beta.TC.SP*TC24*SP + beta.SP.asc24 * SP
V[['B25']] = TC*TC25 + beta.asc25 + Dist*DIST25 +beta.TC.SP*TC25*SP + beta.SP.asc25 * SP
V[['B26']] = TC*TC26 + beta.asc26 + Dist*DIST26 +beta.TC.SP*TC26*SP + beta.SP.asc26 * SP
V[['B27']] = TC*TC27 + beta.asc27 + Dist*DIST27 +beta.TC.SP*TC27*SP + beta.SP.asc27 * SP
V[['B28']] = TC*TC28 + beta.asc28 + Dist*DIST28 +beta.TC.SP*TC28*SP + beta.SP.asc28 * SP
V[['B29']] = TC*TC29 + beta.asc29 + Dist*DIST29 +beta.TC.SP*TC29*SP + beta.SP.asc29 * SP
### Define settings for MNL model component
mnl_settings = list(
alternatives = c('B1' = 1,
'B2' = 2,
'B3' = 3,
'B4' = 4,
'B5' = 5,
'B6' = 6,
'B7' = 7,
'B8' = 8,
'B9' = 9,
'B10' = 10,
'B11' = 11,
'B12' = 12,
'B13' = 13,
'B14' = 14,
'B15' = 15,
'B16' = 16,
'B17' = 17,
'B18' = 18,
'B19' = 19,
'B20' = 20,
'B21' = 21,
'B22' = 22,
'B23' = 23,
'B24' = 24,
'B25' = 25,
'B26' = 26,
'B27' = 27,
'B28' = 28,
'B29' = 29
),
avail = list('B1' = CS1,
'B2' = CS2,
'B3' = CS3,
'B4' = CS4,
'B5' = CS5,
'B6' = CS6,
'B7' = CS7,
'B8' = CS8,
'B9' = CS9,
'B10' = CS10,
'B11' = CS11,
'B12' = CS12,
'B13' = CS13,
'B14' = CS14,
'B15' = CS15,
'B16' = CS16,
'B17' = CS17,
'B18' = CS18,
'B19' = CS19,
'B20' = CS20,
'B21' = CS21,
'B22' = CS22,
'B23' = CS23,
'B24' = CS24,
'B25' = CS25,
'B26' = CS26,
'B27' = CS27,
'B28' = CS28,
'B29' = CS29
) ,
choiceVar = choice,
V = lapply(V, "*", mu_RP) ,
rows = (SP==0)
)
P[["RP"]] = apollo_mnl(mnl_settings,functionality)
mnl_settings$V = lapply(V, "*", mu_SP)
mnl_settings$rows = (SP==1)
P[["SP"]] = apollo_mnl(mnl_settings,functionality)
### Combined model
P = apollo_combineModels(P, apollo_inputs, functionality)
### Average across intra-individual draws
#P = apollo_avgIntraDraws(P, apollo_inputs, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
### Average across inter-individual draws
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
### Optional speedTest
#speedTest_settings=list(
# nDrawsTry = c(50, 75, 100),
# nCoresTry = 1:3,
# nRep = 10
#)
#apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
#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(hessianRoutine="maxLik"))
apollo_modelOutput(model, list(printPVal=TRUE, printDiagnostics=TRUE))
before P[["SP"]], the code runs well.
It is still strange that this setting needs to be specified for a different number of Halton draws.