### R code from vignette source 'RSGHB_HowTo.Rnw' ################################################### ### code chunk number 1: RSGHB_HowTo.Rnw:44-72 (eval = FALSE) ################################################### ## ## library(RSGHB) ## ## # Load example data ## data(choicedata) ## ## # This data set is organized as one row per choice observation. This isn't necessary ## # but it does make writing the likelihood function straightforward. ## head(choicedata) ## ## # Any variables or transformations of variables from the choicedata data.frame that are ## # needed for evaluating the likelihood function can be extracted as a series of vectors ## # for convenience (e.g., TT1, TT2, TOLL2). Alternatively, the likelihood function could ## # refer to the choicedata data.frame directly. In this example, each alternative is ## # defined by travel times and toll costs. ## TT1 <- choicedata$tt1 ## TT2 <- choicedata$tt2 ## TOLL2 <- choicedata$toll2 ## ## # Similarly for the vectors of choices. Note that in this example, there are only two ## # alternatives. Dummying coding the choice vector is not necessary but again makes ## # writing the likelihood function straightforward. ## choice1 <- (choicedata$Choice==1) ## choice2 <- (choicedata$Choice==2) ## ## # Frequency of choices ## table(choicedata$Choice) ## ################################################### ### code chunk number 2: RSGHB_HowTo.Rnw:83-102 (eval = FALSE) ################################################### ## ## # The likelihood function ## likelihood <- function(fc, b) { ## ## # Assign fixed parameters to named variables for convenience ## cc <- 1 ## ASC1 <- fc[cc]; cc <- cc + 1 ## Btime <- fc[cc]; cc <- cc + 1 ## Btoll <- fc[cc]; cc <- cc + 1 ## ## # Utility functions ## v1 <- ASC1 + Btime * TT1 ## v2 <- Btime * TT2 + Btoll * TOLL2 ## ## # MNL probability statement ## p <- (exp(v1) * choice1 + exp(v2) * choice2) / (exp(v1) + exp(v2)) ## ## return(p) ## } ################################################### ### code chunk number 3: RSGHB_HowTo.Rnw:109-146 (eval = FALSE) ################################################### ## ## ### Setting control list for estimation (see ?doHB for more estimation options) ## ## # The model name/description ## modelname <- "MNL" ## ## # gVarNamesFixed contains the names for the fixed (non-random) variables in the model ## # These will be used in the output and also when displaying iteration detail to ## # the screen ## gVarNamesFixed <- c("ASC1", "BTime", "BCost") ## ## # FC contains the starting values for the fixed coefficients ## FC <- c(0, 0, 0) ## ## # gNCREP contains the number of iterations to use prior to convergence ## gNCREP <- 2500 ## # gNEREP contains the number of iterations to keep for averaging after convergence ## # has been reached ## gNEREP <- 2500 ## # gNSKIP contains the number of iterations between retaining draws for averaging ## gNSKIP <- 1 ## # gINFOSKIP controls how frequently to print info about the iteration process ## gINFOSKIP <- 10 ## # gSeed ensures reproducible results ## gSeed <- 1987 ## ## # To simplify the doHB function call, all of the control arguments are placed in ## # a single list that can be passed directly to doHB ## control <- list(modelname = modelname, ## gVarNamesFixed = gVarNamesFixed, ## FC = FC, ## gNCREP = gNCREP, ## gNEREP = gNEREP, ## gNSKIP = gNSKIP, ## gINFOSKIP = gINFOSKIP, ## gSeed = gSeed) ## ################################################### ### code chunk number 4: RSGHB_HowTo.Rnw:157-163 (eval = FALSE) ################################################### ## ## # Estimate model ## control$verbose <- FALSE ## control$nodiagnostics <- TRUE ## model <- doHB(likelihood, choicedata, control) ## ################################################### ### code chunk number 5: RSGHB_HowTo.Rnw:185-189 (eval = FALSE) ################################################### ## ## # Model iteration details ## head(model[["iter.detail"]]) ## ################################################### ### code chunk number 6: RSGHB_HowTo.Rnw:194-198 (eval = FALSE) ################################################### ## ## # Plot model statistics ## plot(model) ## ################################################### ### code chunk number 7: RSGHB_HowTo.Rnw:203-207 (eval = FALSE) ################################################### ## ## # Plot parameter estimates (see ?plot.RSGHB for more uses) ## plot(model, type = "F") ## ################################################### ### code chunk number 8: RSGHB_HowTo.Rnw:222-243 (eval = FALSE) ################################################### ## ## likelihood <- function(fc, b) { ## ## # Note that the fc argument is still supplied, but is unused ## ## # Using b instead of fc is the only change ## cc <- 1 ## ASC1 <- b[, cc]; cc <- cc + 1 ## Btime <- b[, cc]; cc <- cc + 1 ## Btoll <- b[, cc]; cc <- cc + 1 ## ## # Utility functions ## v1 <- ASC1 + Btime * TT1 ## v2 <- Btime * TT2 + Btoll * TOLL2 ## ## # MNL probability statement ## p <- (exp(v1) * choice1 + exp(v2) * choice2) / (exp(v1) + exp(v2)) ## ## return(p) ## } ## ################################################### ### code chunk number 9: RSGHB_HowTo.Rnw:250-301 (eval = FALSE) ################################################### ## ### Setting control list for estimation (see ?doHB for more estimation options) ## ## # The model name/description ## modelname <- "MMNL" ## ## # gVarNamesNormal provides names for the random parameters in the same way ## # gVarNamesFixed does for the fixed parameters ## gVarNamesNormal <- c("ASC1","BTime","BCost") ## ## # svN contains the starting values for the means of the normal distributions for each ## # of the random parameters ## svN <- c(0, 0, 0) ## ## # gDIST specifies the type of continuous distribution to use for the random parameters ## # gDIST must have an entry for each value in gVarNamesNormal ## # The options are: ## # 1. normal ## # 2. log-nomal ## # 3. negative log-normal ## # 4. normal with all values below zero massed at zero ## # 5. normal with all values greater than zero massed at zero ## # 6. Johnson SB with a specified min and max ## ## # In this example, normal distributions are used for all 3 parameters ## gDIST <- c(1, 1, 1) ## ## # gNCREP contains the number of iterations to use prior to convergence ## gNCREP <- 2500 ## # gNEREP contains the number of iterations to keep for averaging after convergence ## # has been reached ## gNEREP <- 2500 ## # gNSKIP contains the number of iterations between retaining draws for averaging ## gNSKIP <- 1 ## # gINFOSKIP controls how frequently to print info about the iteration process ## gINFOSKIP <- 10 ## # gSeed ensures reproducible results ## gSeed <- 1987 ## ## # To simplify the doHB function call, all of the control parameters are placed in ## # a single list that can be passed directly to doHB ## control <- list(modelname = modelname, ## gVarNamesNormal = gVarNamesNormal, ## gDIST = gDIST, ## svN = svN, ## gNCREP = gNCREP, ## gNEREP = gNEREP, ## gNSKIP = gNSKIP, ## gINFOSKIP = gINFOSKIP, ## gSeed = gSeed) ## ## ################################################### ### code chunk number 10: RSGHB_HowTo.Rnw:312-318 (eval = FALSE) ################################################### ## ## # Estimate model ## control$verbose <- FALSE ## control$nodiagnostics <- TRUE ## model <- doHB(likelihood, choicedata, control) ## ################################################### ### code chunk number 11: RSGHB_HowTo.Rnw:336-343 (eval = FALSE) ################################################### ## ## # Model iteration details ## head(model[["iter.detail"]]) ## ## # Plot model statistics ## plot(model) ## ################################################### ### code chunk number 12: RSGHB_HowTo.Rnw:348-355 (eval = FALSE) ################################################### ## ## # Sample-level means ## head(model[["A"]]) ## ## # Plot sample-level means ## plot(model, type = "A") ## ################################################### ### code chunk number 13: RSGHB_HowTo.Rnw:360-367 (eval = FALSE) ################################################### ## ## # Average individual-level draws ## head(model[["B"]]) ## ## # Standard deviations of individual-level draws ## head(model[["Bsd"]]) ## ################################################### ### code chunk number 14: RSGHB_HowTo.Rnw:374-381 (eval = FALSE) ################################################### ## ## # Average individual-level draws (transformed; if applicable) ## head(model[["C"]]) ## ## # Standard deviations of individual-level draws (transformed; if applicable) ## head(model[["Csd"]]) ##