title: How to run a multinomial logit regression with both individual and time fixed effects in R
body: <p>Long story short:</p>
  
  <p>I need to run a multinomial logit regression with both individual and time fixed effects in R. 
  I thought I could use the packages <code>mlogit</code> and <code>survival</code> to this purpose, but I am cannot find a way to include fixed effects.</p>
  
  <p>Now the long story:</p>
  
  <p>I have found many questions on this topic on various stack-related websites, none of them were able to provide an answer. Also, I have noticed a lot of confusion regarding what a multinomial logit regression with fixed effects is (people use different names) and about the R packages implementing this function.
  So I think it would be beneficial to provide some background before getting to the point.</p>
  
  <p>Consider the following.
  In a multiple choice question, each respondent take one choice.
  Respondents are asked the same question every year. There is no apriori on the extent to which choice at time t is affected by the choice at t-1.
  Now imagine to have a panel data recording these choices. The data, would look like this:</p>
  
  <pre><code>set.seed(123)
  # number of observations
  n &lt;- 100
  # number of possible choice
  possible_choice &lt;- letters[1:4]
  # number of years
  years &lt;- 3
  # individual characteristics
  x1 &lt;- runif(n * 3, 5.0, 70.5)
  x2 &lt;- sample(1:n^2, n * 3, replace = F)
  # actual choice at time 1
  actual_choice_year_1 &lt;- possible_choice[sample(1:4, n, replace = T, prob = rep(1/4, 4))]
  actual_choice_year_2 &lt;- possible_choice[sample(1:4, n, replace = T, prob = c(0.4, 0.3, 0.2, 0.1))]
  actual_choice_year_3 &lt;- possible_choice[sample(1:4, n, replace = T, prob = c(0.2, 0.5, 0.2, 0.1))]
  # create long dataset
  df &lt;- data.frame(choice = c(actual_choice_year_1, actual_choice_year_2, actual_choice_year_3),
             x1 = x1, x2 = x2, 
             individual_fixed_effect = as.character(rep(1:n, years)),
             time_fixed_effect = as.character(rep(1:years, each = n)),
             stringsAsFactors = F)
  </code></pre>
  
  <p>I am new to this kind of analysis. But if I understand correctly, if I want to estimate the effects of respondents' characteristics on their choice, I may use a multinomial logit regression.</p>
  
  <p>In order to take advantage of the longitudinal structure of the data, I want to include in my specification individual and time fixed effects.</p>
  
  <p>To the best of my knowledge, the multinomial logit regression with fixed effects was first proposed by Chamberlain (1980, Review of Economic Studies 47: 225–238). Recently, Stata users have been provided with the routines to implement this model (femlogit).</p>
  
  <p>In the vignette of the femlogit package, the author refers to the R function <code>clogit</code>, in the <code>survival</code> package. </p>
  
  <p>According to the help page, <code>clogit</code> requires data to be rearranged in a different format:</p>
  
  <pre><code>library(mlogit)
  # create wide dataset
  data_mlogit &lt;- mlogit.data(df, id.var = "individual_fixed_effect", 
              group.var = "time_fixed_effect", 
              choice = "choice", 
              shape = "wide")
  </code></pre>
  
  <p>Now, if I understand correctly how <code>clogit</code> works, fixed effects can be passed through the function <code>strata</code> (see for additional details this <a href="https://data.princeton.edu/wws509/r/fixedrandom3" rel="noreferrer">tutorial</a>). However, I am afraid that it is not clear to me how to use this function, as no coefficient values are returned for the individual characteristic variables (i.e. I get only <code>NAs</code>).</p>
  
  <pre><code>library(survival)
  fit &lt;- clogit(formula("choice ~ alt + x1 + x2 + strata(individual_fixed_effect, time_fixed_effect)"), as.data.frame(data_mlogit))
  summary(fit)
  </code></pre>
  
  <p>Since I was not able to find a reason for this (there must be something that I am missing on the way these functions are estimated), I have looked for a solution using other packages in R: e.g., <code>glmnet</code>, <code>VGAM</code>, <code>nnet</code>, <code>globaltest</code>, and <code>mlogit</code>.</p>
  
  <p>Only the latter seems to be able to explicitly deal with panel structures using appropriate estimation strategy. For this reason, I have decided to give it a try. However, I was only able to run a multinomial logit regression without fixed effects.</p>
  
  <pre><code># state formula
  formula_mlogit &lt;- formula("choice ~ 1| x1 + x2")
  
  # run multinomial regression
  fit &lt;- mlogit(formula_mlogit, data_mlogit)
  summary(fit)
  </code></pre>
  
  <p>If I understand correctly how <code>mlogit</code> works, here's what I have done.</p>
  
  <p>By using the function <code>mlogit.data</code>, I have created a dataset compatible with the function <code>mlogit</code>. Here, I have also specified the id of each individual (<code>id.var = individual_fixed_effect</code>) and the group to which individuals belongs to (<code>group.var = "time_fixed_effect"</code>). In my case, the group represents the observations registered in the same year.</p>
  
  <p>My formula specifies that there are no variables correlated with a specific choice, and which are randomly distributed among individuals (i.e., the variables before the <code>|</code>). By contrast, choices are only motivated by individual characteristics (i.e., <code>x1</code> and <code>x2</code>).</p>
  
  <p>In the help of the function <code>mlogit</code>, it is specified that one can use the argument <code>panel</code> to use panel techniques. To set <code>panel = TRUE</code> is what I am after here.</p>
  
  <p>The problem is that <code>panel</code> can be set to <code>TRUE</code> only if another argument of <code>mlogit</code>, i.e. <code>rpar</code>, is not <code>NULL</code>.</p>
  
  <p>The argument <code>rpar</code> is used to specify the distribution of the random variables: i.e. the variables before the <code>|</code>.
  The problem is that, since these variables does not exist in my case, I can't use the argument <code>rpar</code> and then set <code>panel = TRUE</code>.</p>
  
  <p>An interesting question related to this is <a href="https://stats.stackexchange.com/questions/51148/unable-to-provide-random-parameter-with-mlogit">here</a>. A few suggestions were given, and one seems to go in my direction. Unfortunately, no examples that I can replicate are provided, and I do not understand how to follow this strategy to solve my problem.</p>
  
  <p>Moreover, I am not particularly interested in using <code>mlogit</code>, any efficient way to perform this task would be fine for me (e.g., I am ok with <code>survival</code> or other packages).</p>
  
  <p>Do you know any solution to this problem?</p>
  
  <p>Two caveats for those interested in answering:</p>
  
  <ol>
  <li>I am interested in fixed effects, not in random effects. However, if you believe there is no other way to take advantage of the longitudinal structure of my data in R (there is indeed in Stata but I don't want to use it), please feel free to share your code.</li>
  <li>I am not interested in going Bayesian. So if possible, please do not suggest this approach.</li>
  </ol>