A popular method for reducing the mean and median bias of the maximum likelihood estimator in regular parametric models is through the additive adjustment of the score equation (Firth, 1993; Kenne Pagui et al., 2017). The current work focuses on mean and median bias-reducing adjusted score equations in models with latent variables. First, we give estimating equations based on a mean bias-reducing adjustment of the score function for mean bias reduction in linear mixed models. Second, we propose an extension of the adjusted score equation approach (Firth, 1993) to obtain bias-reduced estimates for models with either computationally infeasible adjusted score equations and/or intractable likelihood. The proposed bias-reduced estimator is obtained by solving an approximate adjusted score equation, which uses an approximation of the log-likelihood to obtain tractable derivatives, and Monte Carlo approximation of the bias function to get feasible expressions. Under certain general conditions, we prove that the feasible and tractable bias-reduced estimator is consistent and asymptotically normally distributed. The “iterated bootstrap with likelihood adjustment” algorithm is presented that can compute the solution of the new bias-reducing adjusted score equation. The effectiveness of the proposed method is demonstrated via simulation studies and real data examples in the case of generalised linear models and generalised linear mixed models. Finally, we derive the median bias-reducing adjusted scores for linear mixed models and random-effects meta-analysis and meta-regression models
