We present mexhaz, an R package for fitting flexible hazard-based regression models with the possibility to add time-dependent effects of covariates and to account for a two-level hierarchical structure in the data through the inclusion of a normally distributed random intercept (i.e., a log-normally distributed shared frailty). Moreover, mexhaz-based models can be fitted within the excess hazard setting by allowing the specification of an expected hazard in the model. These models are of common use in the context of the analysis of population-based cancer registry data. Follow-up time can be entered in the right-censored or counting process input style, the latter allowing models with delayed entries. The logarithm of the baseline hazard can be flexibly modeled with B-splines or restricted cubic splines of time. Parameters estimation is based on likelihood maximization: in deriving the contribution of each observation to the cluster-specific conditional likelihood, Gauss-Legendre quadrature is used to calculate the cumulative hazard; the cluster-specific marginal likelihoods are then obtained by in-tegrating over the random effects distribution, using adaptive Gauss-Hermite quadrature. Functions to compute and plot the predicted (excess) hazard and (net) survival (possibly with cluster-specific predictions in the case of random effect models) are provided. We illustrate the use of the different options of the mexhaz package and compare the results obtained with those of other available R packages.