Consider the following xtmelogit.

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear

xtmelogit honors i.female read || cid:

Refining starting values: 

Iteration 0:   log likelihood = -82.156921  
Iteration 1:   log likelihood = -80.132249  
Iteration 2:   log likelihood = -79.765159  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -79.765159  
Iteration 1:   log likelihood = -79.721813  
Iteration 2:   log likelihood =  -79.72173  
Iteration 3:   log likelihood =  -79.72173  

Mixed-effects logistic regression               Number of obs      =       200
Group variable: cid                             Number of groups   =        20

                                                Obs per group: min =         7
                                                               avg =      10.0
                                                               max =        12

Integration points =   7                        Wald chi2(2)       =      6.39
Log likelihood =  -79.72173                     Prob > chi2        =    0.0410

------------------------------------------------------------------------------
      honors |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.female |   1.203685    .498816     2.41   0.016     .2260232    2.181346
        read |   .0426387   .0579243     0.74   0.462    -.0708908    .1561682
       _cons |  -4.827711   2.919208    -1.65   0.098    -10.54925    .8938321
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Identity                |
                   sd(_cons) |   2.097709   .9227794      .8857332     4.96807
------------------------------------------------------------------------------
LR test vs. logistic regression: chibar2(01) =    11.44 Prob>=chibar2 = 0.0004

Say that we wanted to know the predicted probabilities for males and females at five
different values of read, 30, 40, 50, 60 and 70. We might try using the margins
command.

margins female, at(read=(30(10)70)) vsquish

default prediction is a function of possibly stochastic quantities other than e(b)
r(498);

That didn’t work because margins can’t compute predicted values for models that
have both fixed and random components.

We can, however, get the predicted probabilities from just the fixed part of the model,
like this.

margins female, at(read=(30(10)70)) predict(mu fixedonly) vsquish

Adjusted predictions                              Number of obs   =        200

Expression   : Predicted mean, fixed portion only, predict(mu fixedonly)
1._at        : read            =          30
2._at        : read            =          40
3._at        : read            =          50
4._at        : read            =          60
5._at        : read            =          70

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  _at#female |
        1 0  |    .027962   .0353643     0.79   0.429    -.0413509    .0972748
        1 1  |   .0874748   .1006056     0.87   0.385    -.1097085    .2846582
        2 0  |   .0422023   .0352079     1.20   0.231    -.0268039    .1112085
        2 1  |   .1280315   .0900429     1.42   0.155    -.0484494    .3045123
        3 0  |   .0632231   .0416911     1.52   0.129      -.01849    .1449363
        3 1  |   .1836082   .0928378     1.98   0.048     .0016494     .365567
        4 0  |   .0936902   .0807071     1.16   0.246    -.0644929    .2518733
        4 1  |   .2562211   .1691374     1.51   0.130     -.075282    .5877243
        5 0  |   .1366968   .1661537     0.82   0.411    -.1889584     .462352
        5 1  |   .3454012   .3085949     1.12   0.263    -.2594336    .9502361
------------------------------------------------------------------------------

For many purposes these probabilities from the fixed effects only will be all that we
will need and these probabilities could be graphed using marginsplot. If we are specifically interested in the estimated of probabilities that
include both fixed and random effects we can make use of the predict command.

First, we will estimate the predicted probabilities from the fixed and random parts of the
model directly.

predict mu1, mu

tabstat  mu1, by(female)

Summary for variables: mu1
     by categories of: female 

female |      mean
-------+----------
  male |  .1909927
female |  .3182445
-------+----------
 Total |  .2603449
------------------

The values produced by tabstat give us the predicted probabilities separately for males
and females while read is held at its observed values. We can estimate these
same values in two steps by estimating the linear predictor for the random and fixed
effects separately.

predict re*, reffects  //  linear predictor for the random effects

predict xb, xb         // linear predictor for the fixed effects

gen mu2 = 1 /(1+exp(-1*(xb + re1)))  // compute probabilities using both fixed and random components

tabstat mu2, by(female) 

Summary for variables: mu2
     by categories of: female 

female |      mean
-------+----------
  male |  .1909927
female |  .3182445
-------+----------
 Total |  .2603449
------------------

As you can see the predicted probabilities computed this way are the same as when
we predicted mu1 directly.
The term (xb + re1) combines the fixed and random linear predictors while
1 /(1+exp(-1*(xb + re1))) converts the predictions to the probability metric.

We will now use the same approach to fix read at its mean value while letting
female vary as observed.

sum read
replace read = r(mean)
predict xb3, xb
gen mu3 = 1 /(1+exp(-1*(xb3 + re1)))

tabstat mu3, by(female)

Summary for variables: mu3
     by categories of: female 

female |      mean
-------+----------
  male |  .1588983
female |  .3037404
-------+----------
 Total |  .2378372
------------------

We can use this same method to compute predicted probabilities for each gender at the
five values of read of interest by computing a separate xb for each combination of
values. We can continue to use the the linear random predictor, re1, because the observations
remain nested in the same cid. We will do the computations in a forvalues loop.

forvalues i=30(10)70 {
   quietly replace read = `i'
   quietly replace female = 0  // begin with the males
   predict xbm`i' , xb
   gen mum`i' = 1 /(1+exp(-1*(xbm`i' + re1))) 
   quietly replace female = 1  // now for the females
   predict xbf`i' , xb
   gen muf`i' = 1 /(1+exp(-1*(xbf`i' + re1))) 
}
 
sum mum* muf*

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       mum30 |       200    .0878644    .1277633   .0062536    .452949
       mum40 |       200    .1181626     .161502   .0095469   .5591283
       mum50 |       200    .1545762    .1959134   .0145493    .660161
       mum60 |       200    .1969638    .2284821   .0221143   .7484569
       mum70 |       200    .2450655    .2570279   .0334791   .8200648
-------------+--------------------------------------------------------
       muf30 |       200     .189036    .2229538   .0205396   .7339823
       muf40 |       200    .2361479    .2523496   .0311209   .8086573
       muf50 |       200    .2887174    .2762624   .0468924   .8661915
       muf60 |       200    .3463831    .2931488   .0700785   .9083859
       muf70 |       200    .4087083    .3014779   .1034842   .9382238

The variables that begin mum give the predicted values for the males for each of
the five different values of read. The muf do the same thing for the
females. We will plot these values by collapsing the data to a single observation and
then reshaping the data long.

collapse (mean) mum* muf*
gen tid = 1
reshape long mum muf, i(tid) j(read)

twoway (line mum read)(line muf read), legend(order(1 "male" 2 "female")) ///
  ytitle(predicted probability) ///
  title("Predicted probabilities with both" "fixed and random effects") ///
  name(mu, replace)

With this approach you can replace the observed values with the mean values or any value that you
wish. So, this approach worked okay for random intercepts, how would it work if we include a random
coefficient (slope).

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear

xtmelogit honors i.female read || cid: female, cov(unstr)

Refining starting values: 

Iteration 0:   log likelihood = -81.308446  
Iteration 1:   log likelihood = -79.509044  
Iteration 2:   log likelihood = -79.233888  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -79.233888  
Iteration 1:   log likelihood = -79.203192  
Iteration 2:   log likelihood =  -79.20305  
Iteration 3:   log likelihood =  -79.20305  

Mixed-effects logistic regression               Number of obs      =       200
Group variable: cid                             Number of groups   =        20

                                                Obs per group: min =         7
                                                               avg =      10.0
                                                               max =        12

Integration points =   7                        Wald chi2(2)       =      2.07
Log likelihood =  -79.20305                     Prob > chi2        =    0.3544

------------------------------------------------------------------------------
      honors |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.female |   1.107571   .9144827     1.21   0.226    -.6847819    2.899924
        read |   .0548958   .0614466     0.89   0.372    -.0655374    .1753289
       _cons |  -5.459728   3.237309    -1.69   0.092    -11.80474    .8852801
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Unstructured            |
                  sd(female) |   1.270624    .877583      .3281837    4.919457
                   sd(_cons) |   2.104286   1.168272      .7088079    6.247133
          corr(female,_cons) |  -.1370656   1.016493     -.9741795    .9555909
------------------------------------------------------------------------------
LR test vs. logistic regression:     chi2(3) =    12.48   Prob > chi2 = 0.0059

Note: LR test is conservative and provided only for reference.

predict mu1, mu

tabstat mu1, by(female)

Summary for variables: mu1
     by categories of: female 

female |      mean
-------+----------
  male |  .1884487
female |  .3188527
-------+----------
 Total |  .2595188
------------------

To reproduce the predicted probabilities that include both fixed and random effects, we
need to obtain the linear predictor of the random effects as follows:

predict re*, reffects // obtain the random effects

des re*

              storage  display     value
variable name   type   format      label      variable label
--------------------------------------------------------------------------------------------
read            byte   %9.0g                  reading score
re1             float  %9.0g                  random effects for cid: female
re2             float  %9.0g                  random effects for cid: _cons

As you can see, re1 is the linear predictor for the random coefficient while re2
linear the linear predictor for the random intercept.

Now we need to compute the interaction of female with re1, obtain the
linear predictor for the fixed effect and compute the probability.

gen fre1=female*re1

predict xb2, xb

gen mu2 = 1 /(1+exp(-1*(xb2 + fre1 + re2)))

tabstat mu2, by(female)

Summary for variables: mu2
     by categories of: female 

female |      mean
-------+----------
  male |  .1884487
female |  .3188527
-------+----------
 Total |  .2595188
------------------

Once again mu1 and mu2 are the same.

Next, with this knowledge, we are going to compute the predicted probabilities for
read for the values
30, 40 50 60 and 70 for both males and females, just like we did in the first section.
This time we need to compute fre1 for both males and female. For male the
value of female*re1 is zero. While for female the value of female*re1 is
equal to re1.

forvalues i=30(10)70 {
  quietly replace read = `i'
  quietly replace female = 0
  quietly replace fre1 = 0  // this value is 0 for males
  predict xbm`i' , xb
  gen mum`i' = 1 /(1+exp(-1*(xbm`i' +fre1 + re2))) 
  quietly replace female = 1
  quietly replace fre1 = re1  // this value is re1 for females
  predict xbf`i' , xb
  gen muf`i' = 1 /(1+exp(-1*(xbf`i' + fre1 + re2))) 
}

sum mum* muf*

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       mum30 |       200     .067521    .1073255   .0062084   .4335063
       mum40 |       200    .1006205    .1453719    .010701   .5698889
       mum50 |       200    .1441323    .1850342   .0183843   .6964301
       mum60 |       200    .1988144    .2226292    .031409   .7988807
       mum70 |       200    .2647876    .2545738   .0531617   .8730579
-------------+--------------------------------------------------------
       muf30 |       200    .1651364    .2082666   .0140494   .5906218
       muf40 |       200    .2210162     .252848   .0240784   .7141231
       muf50 |       200    .2841479    .2869894    .040969   .8122126
       muf60 |       200    .3544923    .3071029   .0688718   .8821978
       muf70 |       200     .432452    .3113492   .1135287   .9283999

As before we will collapse the data then reshape them prior to plotting.

collapse (mean) mum* muf*
gen tid = 1
reshape long mum muf, i(tid) j(read)

twoway (line mum read)(line muf read), legend(order(1 "male" 2 "female")) ///
  ytitle(predicted probability) ///
  title("Predicted probabilities with both" "fixed and random effects") ///
  name(mu2, replace)

The two graphs that we generated for this page look very similar but inspection of the
tables for the predicted probabilities show that adding the random coefficient to the model
changes the predicted probabilities. Let’s check the two models using a likelihood-ratio test.

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear

quietly xtmelogit honors i.female read || cid: female, cov(unstr)
estimates store m1
quietly xtmelogit honors i.female read || cid:
estimates store m2

lrtest m1 m2

Likelihood-ratio test                                 LR chi2(2)  =      1.04
(Assumption: m1 nested in m2)                         Prob > chi2 =    0.5953

Note: The reported degrees of freedom assumes the null hypothesis is not on the boundary
      of the parameter space.  If this is not true, then the reported test is 
      conservative.

So, in this instance, the model with the random coefficient is not significantly better than the model with just
the random intercept.

Acknowledgements: Thanks to Jurjen Iedema of the Netherlands Institute for Social Research / SCP
for suggesting this approach.