Firth's regression with many fixed effects I've got a sample (n=7000) with a binary dependent variable with many 0s and only a few 1s (150 total) and therefore want to run a Firth's biased reduced regression to estimate the probability of the 1s. So far so good, but my explanatory variables are fixed effects for the individuals (=their individual names). Each individual has between 2 and 50 cases in the data, and all in all there are about 120 individuals. An individual's covariate values can change from trial to trial. As I want to include the individuals as fixed effects I convert them as Dummies. Running Firth's regression with the R logistf package produces several errors now, because of the big number of 'IndividualDummies'. Is there a way to include a big amount of dummies in another way like (1|Individual) in glmer package? Would be great to make it more handy for Firth's regression.

The big problem here is the small number of events per predictor, as you want to include the individuals as fixed effects. It's not clear that the Firth penalization is the best solution to that problem.
To avoid overfitting you typically need about 10-20 cases in the minority class (events) per predictor  in the model. With only 150 events and 120 individuals treated as fixed effects, plus other covariates, you are approaching just 1 event per predictor. Some type of penalization is called for, but it's not clear that Firth's is the best choice.
First, the original Firth method penalizes both the regression coefficients and the intercept toward values of 0. As it reduces small-sample bias in predictor coefficients it thus also biases the intercept toward 0 so that probability predictions are biased toward 0.5. The logistf package now provides modifications that help avoid that problem.
Second, penalization is tricky with dummy variables.* As Frank Harrell explains in Section 9.10 of the second edition of Regression Modeling Strategies, with penalization "one will get different predictions depending on which cell is chosen as the reference cell when constructing dummy variables." For a categorical predictor having multiple levels, he recommends penalizing the sum of squared differences of the individual regression coefficients from their mean. The lrm() and pentrace() functions in his rms package can implement that approach in the broader context of penalized maximum likelihood.
Another choice would be a simple mixed-effect model, with the individuals treated as random effects and fixed effects limited to the covariates. In that case, the much higher event/predictor ratio might get around the small-sample coefficient bias for covariates adequately for your purposes. Random intercepts for the individuals would represent their differences from the log-odds at baseline covariate values and thus still provide some basis for examining differences among them.
If you are intent on treating the individuals as fixed effects, I suppose that Björn's suggestion to use a Bayesian approach (+1) could be extended to a fixed-effects model.
Finally, I'd suggest that you read this related thread. The answer from StasK shows how to treat clusters (individuals in your case) as fixed effects with a Firth regression.* If you want to use a Bayesian approach in a mixed model, the answer from Ben Bolker illustrates how to use the blme package, a Bayesian-modeling wrapper for lme4, to use a weak prior on the fixed effects to accomplish something similar to the Jeffreys prior that corresponds to Firth penalization.

*Did you set up the dummies yourself, or did you set up a single Individual categorical predictor including all individuals and let R define the corresponding dummies? The latter is much less error-prone. Errors in defining dummy variables can lead to error messages when you try to fit a model.
