Discrete choice simulations

multchoices <- function(nCoefs = 5, nObs = 1000, nChoices = 2, grpSize = 5, error = "varying") {
	require(mclogit)
	require(dplyr)
	require(broom)
	require(fExtremes)
	require(mitools)
	require(survival)
### Generate coefs
    coefs <- runif(n = nCoefs, min = -10, max = 10)
    nObs <- nObs
    grpSize <- grpSize
    if (nObs %% grpSize != 0) {
        stop("Observations must be perfectly divisible by group size")
    }
    ### Set up the RHS
    X <- rnorm(n = nObs * nCoefs)
    X <- matrix(X, nrow = nObs, ncol = nCoefs)
### Set up the groups
    group <- rep(1:(nObs / grpSize), each = grpSize)
    ystar <- tcrossprod(coefs, X)
    ystar <- as.vector(ystar)
    dat <- data.frame(ystar = ystar, group = group, chosen = 0)

### If the error is fixed we can add it on now
    if (error == "fixed") {
        e <- rgev(n = nObs,
            xi = 0, ## Gumbel
            mu = 0, beta = 1)
        dat$mu <- dat$ystar + e
        dat <- dat %>%
            group_by(group) %>%
            mutate(chosen = as.numeric(mu >= sort(mu, decreasing = TRUE)[nChoices]))
    } else {
        for (i in 1:nChoices) {
            e <- rgev(n = nObs,
                      xi = 0, ## Gumbel
                      mu = 0, beta = 1)
            dat$mu <- dat$ystar + e
### Get the chosen entry from the currently non-chosen
            dat <- dat %>%
                group_by(group) %>%
                mutate(chosen = pmax(chosen, as.numeric(mu == max(mu[chosen == 0]))))
        }
    }
### We have the data, now estimate the models
    dat$ystar <- NULL
    dat$mu <- NULL
    X <- as.data.frame(X)
    class(dat) <- "data.frame"
    dat <- cbind(dat, X)
    
    mclogit_mod <- tidy(mclogit(cbind(chosen, group) ~ . - group, data = dat))

    mclogit_mod$method <- "mclogit"

    true <- data.frame(term = paste0("V", 1:nCoefs),
                       estimate = coefs,
                       std.error = 0,
                       statistic = NA,
                       p.value = NA,
                       method = "known")
    retval <- merge(mclogit_mod, true, all = TRUE)
    retval$nChoices <- nChoices
    retval$error <- error
    retval$nCoefs <- nCoefs
    return(retval)
}

set.seed(1982)
a <- replicate(100, multchoices(nChoices = 1),simplify = FALSE)
b <- replicate(100, multchoices(nChoices = 2),simplify = FALSE)

a <- do.call("rbind", a)
b <- do.call("rbind", b)
summary(abs(a$estimate) / abs(b$estimate))