Simulating Multinomial Logit Data with R

Question

I'm looking to generate fake data to fit a multinomial logit in R? Any code/suggestions on material to look at would be very much appreciated ...

Answer 1

It is really simple to generate multinomial logit regression data. All you need to keep in mind are the normalizing assumptions.

# covariate matrix
mX = matrix(rnorm(1000), 200, 5)

# coefficients for each choice
vCoef1 = rep(0, 5)
vCoef2 = rnorm(5)
vCoef3 = rnorm(5)

# vector of probabilities
vProb = cbind(exp(mX%*%vCoef1), exp(mX%*%vCoef2), exp(mX%*%vCoef3))

# multinomial draws
mChoices = t(apply(vProb, 1, rmultinom, n = 1, size = 1))
dfM = cbind.data.frame(y = apply(mChoices, 1, function(x) which(x==1)), mX)

Here mChoices and dfM$y encode the same information differently.

Answer 2

In this model with $k$ variables and $d$ categories, let $x = (1,x_1, x_2,\ldots, x_k)$ be one data point including a constant $1$ for the intercept. There are $d-1$ column vectors $\beta_2, \beta_3, \ldots, \beta_d,$ each of length $k+1,$ that give the relative chances of each category in the ratios

$$p_1:p_2:\cdots:p_d = 1:e^{x\beta_2}:e^{x\beta_3}:\cdots:e^{x\beta_d}.$$

The random response for this data point places a specified number $s$ balls into $d$ bins (one for each category) according to these probabilities.

Thus, to simulate such data you need to

1. Specify the variables--their number $k,$ the number of data points $n,$ and all their values.

2. Specify the number of categories and the $\beta_j$ parameters.

3. Specify the value of $s$ for each data point.

4. Compute the probabilities determined by (1) and (2) according to this model.

5. Use those sizes (3) and probabilities (4) to generate random multinomial outcomes.

This leads to a straightforward R implementation. It randomly generates and stores the variables (1) in an array X and, given a randomly-generated array of coefficients (2) and randomly-generated size array (3), computes the probabilities (4) and applies rmultinom to each data point (5) to obtain the matrix of responses (one column per category) and store it in the array y.

n <- 5e3            # Number of observations
k <- 3              # Number of variables
d <- 4              # Number of categories
size <- 15          # Expected size of each outcome (number of balls selected)

set.seed(17)
xnames <-  paste0("X", seq_len(k))
beta <- matrix(rnorm((k+1) * d), ncol = d, dimnames=list(c("Intercept", xnames), seq_len(d))) 
beta = beta - beta[, 1] # Standardize: category 1 is the reference category
X <- matrix(runif(n * k), n, dimnames = list(NULL, xnames))

p <- (function(h) h / rowSums(h))(exp(cbind(1, X) %*% beta))
s <- 1 + rpois(n, size-1)
y <- t(sapply(seq_len(n), function(i) rmultinom(1, s[i], p[i, ])))

With the $n=5000$ observations specified here, the estimates $\hat\beta$ had better be close to the stipulated value of $\beta$! Using nnet::multinom I estimated the coefficients.

library(nnet)
fit <- multinom(y ~ X)

Here is a plot comparing all $12$ coefficients to their specified values. The gray lines are 95% confidence intervals around each estimate. The agreement is good, indicating this approach agrees with the model assumed by multinom.

The total number of results in each category depend on all the inputs. The expected totals can be computed by multiplying the probabilities $(p_1, \ldots, p_d)$ for each data point times its value of $s$ and adding these up by category:

round(rbind(Expected = s %*% p, Observed = colSums(y)))

           1     2     3     4
         880 28799 13837 31623
Observed 902 28941 13852 31444

The counts observed in this simulation, by category, are close to the expected counts.

This indicates how you can adjust your simulation to achieve a desired count, set of counts, or (as asked in the question) proportion of counts: you can alter the data points, the coefficients, and the sizes as you will. How you do this will depend on what aspects of the situation you are willing to vary. There are so many possibilities that it would take us too far afield to discuss them all here.

Answer 3

#Genarating 500 random numbers with zero mean
x = rnorm(500,0)
#Assigning the values of beta1 and beta2
Beta1 = 2
Beta2 = .5
#Calculation of denominator for probability calculation
Denominator= 1+exp(Beta1*x)+exp(Beta2*x)
#Calculating the matrix of probabilities for three choices
vProb = cbind(1/Denominator, exp(x*Beta1)/Denominator, exp(x*Beta2)/Denominator )
# Assigning the value one to maximum probability and zero for rest to get the appropriate choices for value of x
mChoices = t(apply(vProb, 1, rmultinom, n = 1, size = 1))
# Value of Y and X together
dfM = cbind.data.frame(y = apply(mChoices, 1, function(x) which(x==1)), x)
#Adding library for multinomial logit regression
library("nnet")
#We want zero intercept hence x+0 hence the foumula of regression as below
fit<-(multinom(y ~ x + 0, dfM))
#This function uses first y as base class 
#hence upper probability calculation is changed
summary(fit)
#In case we do not keep intercept as zero
fit2<-multinom(y ~ x, dfM)
summary(fit2)
#This also result intercept very close to zero and non significant
#and value of beta as modeled earlier and significant
#running from mlogit package
library(mlogit)
DM<-mlogit.data(dfM, shape="wide",sep="",choice="y",alt.levels=1:3)
#Do not know why -1 is used at two places. I will appreciate if some one can explain
fit3<-mlogit(y~-1|-1+x,data=DM)
summary(fit3)