The design of experiments for discrete mixed effect models is challenging due to the unavailability of a closed-form expression for the Fisher information matrix (FIM), on which most optimality criteria depend. Existing approaches for the computation of the FIM for those models are all based on approximations of the likelihood. A new method is presented which is based on derivatives of the exact conditional likelihood and which uses Monte Carlo (MC) simulations as well as adaptive Gaussian quadrature (AGQ) to integrate those derivatives over the data and random effects. The method is implemented in R and evaluated with respect to the influence of the tuning parameter, the accuracy of the FIM approximation, and computational complexity. The accuracy evaluation is performed by comparing the expected relative standard errors (RSE) from the MC/AGQ FIM with RSE obtained in a simulation study with four different discrete data models (two binary, one count and one repeated time-to-event model) and three different estimation algorithms. Additionally, the results from the MC/AGQ FIM are compared with expected RSE from a marginal quasi-likelihood (MQL) approximated FIM. The comparison resulted in close agreement between the MC/AGQ-based RSE and empirical RSE for all models investigated, and better performance of MC/AGQ than the MQL approximated FIM for variance parameters. The MC/AGQ method also proved to be well suited to calculate the expected power to detect a group effect for a model with binary outcomes.