One possible approach was suggested by Nakagawa and Schielzeth (2013) where $R^2$ is decomposed into variance explained by fixed effects ($R^2_m$) and random effects ($R^2_c$).
Let's say your model looks like the one below, where $\alpha$ and $\gamma$ are random effects:
$$y_{ijk} = \beta_0 + \sum^M_{m=1} \beta_m X_{mijk} + \alpha_j + \gamma_k + \epsilon_{ijk}$$
if so, you can define $R^2$ for fixed effects as:
$$R^2_m = \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\gamma + \sigma^2_\epsilon}$$
where:
$$\sigma^2_f = var\left( \sum^M_{m=1} \beta_m X_{mijk} \right)$$
And $R^2$ for random effects as:
$$R^2_c = \frac{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\gamma}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\gamma + \sigma^2_\epsilon}$$
And this approach could be extended to GLMM's by decomposing $\sigma^2_\epsilon$ into (Nakagawa and Schielzeth, 2013, p. 137):
(i) multiplicative dispersion ($\omega$), (ii) additive dispersion
($\sigma^2_e$) and (iii) distribution-specific variance ($\sigma^2_d$)
Johnson (2014) suggested how this approach could be used for random slopes models.
Functions for estimating $R^2_m$ and $R^2_c$ can be found in MuMIn package using r.squaredGLMM function.
