Massively categorical variables like this require special handling simply due to the memory (RAM) required for inversion of a large cross-products matrix (assuming that's what is needed here). Depending on the installation and module, R software can be quite limited in how many levels it can invert. There are workarounds, depending on the hardware platform and software you have available.
Here's a link to the best paper I'm aware of about a Bayesian multilevel model for crunching massively categorical variables. Regrettably, I'm not aware of any R modules that do this but they probably exist somewhere. It would be readily programmable in R, Stan or Winbugs -- depending on which software you're most comfortable with:
Bayes may not be your preferred approach. Here are a couple of other things to try:
1) 1,000 levels really isn't that many, not when you think of a truly massive variable like zip codes which can contain 40,000 levels or more (if business zips are included). Other software is able to crunch this many levels (and more) as standard practice on today's laptops. Given that, there must be an R module which can do classic matrix inversion for a few thousand levels. I don't know what that module is.
2) Gain access to a massively parallel platform and use a "divide and conquer" approach to slicing your data up into many smaller datasets, distributing them across the core. The goal is to size them appropriately for the memory capacities of each CPU. Then, on the backend, roll them up to summarize the test statistics that are needed. Here's a link to a paper about this "D&C" methodology: