Is it possible to plot R glmer model predictions using Python?

Question

I have an glmer model in R which I want to plot predictions for. I found the plot_model function from the sjPlot library and it works fine.

Here is a MWE:

library(lattice)

cbpp$response <- sample(c(0,1), replace=TRUE, size=nrow(cbpp))
gm1 <- glmer(response ~ size + incidence + (1 | herd),
              data = cbpp, family = binomial)

For example, calling plot_model(gm1, type = "pred", show.data = TRUE) yields the following figure:

However, I am not familiar with R and I am having a hard time trying to control the plot aesthetics and plotting multiple models into the same figure (already asked a question regarding that issue here). I am familiar with Python and matplotlib and getting these figures to work on a Python environment would be much simpler for me.

I'm guessing one way to accomplish this would be taking the y values (predicted probabilities of fire) from R and exporting them so I could read them in Python in order to plot them against each covariate (evi prev) in this example. However, I am not sure how to do this. Furthermore, I tried to read sjPlot source code to figure out how it plots the predictions but could not figure it out either.

Answer 1

ggpredict() actually returns more values (and along the x-axis, i.e. for the term in question - size in your example - these are even-spaced), but only prints fewer values.

library(lme4)
#> Loading required package: Matrix
library(ggeffects)

cbpp$response <- sample(c(0,1), replace=TRUE, size=nrow(cbpp))
gm1 <- glmer(response ~ size + incidence + (1 | herd), data = cbpp, family = binomial)

pr1 <- ggpredict(gm1, term = "size")

pr1
#> 
#> # Predicted probabilities of response
#> # x = size
#> 
#>   x predicted std.error conf.low conf.high
#>   2     0.632     0.717    0.297     0.875
#>   6     0.610     0.550    0.347     0.821
#>  10     0.587     0.407    0.390     0.759
#>  14     0.563     0.321    0.407     0.708
#>  18     0.539     0.339    0.376     0.695
#>  22     0.515     0.448    0.306     0.719
#>  26     0.491     0.601    0.229     0.758
#>  34     0.444     0.951    0.110     0.837
#> 
#> Adjusted for:
#> * incidence = 1.77
#> *      herd = 0 (population-level)
#> Standard errors are on link-scale (untransformed).

as.data.frame(pr1)
#>     x predicted std.error  conf.low conf.high group
#> 1   2 0.6323758 0.7168742 0.2967912 0.8751705     1
#> 2   4 0.6211339 0.6316777 0.3221952 0.8497229     1
#> 3   6 0.6097603 0.5501862 0.3470481 0.8212222     1
#> 4   8 0.5982662 0.4743133 0.3701925 0.7904902     1
#> 5  10 0.5866630 0.4072118 0.3898523 0.7592017     1
#> 6  12 0.5749627 0.3539066 0.4033525 0.7302266     1
#> 7  14 0.5631779 0.3213384 0.4071542 0.7076259     1
#> 8  16 0.5513213 0.3159857 0.3981187 0.6953669     1
#> 9  18 0.5394060 0.3391396 0.3759558 0.6947993     1
#> 10 20 0.5274456 0.3857000 0.3438768 0.7038817     1
#> 11 22 0.5154536 0.4484344 0.3063836 0.7192510     1
#> 12 24 0.5034437 0.5215385 0.2672889 0.7380720     1
#> 13 26 0.4914299 0.6012416 0.2292244 0.7584368     1
#> 14 28 0.4794260 0.6852450 0.1938167 0.7791488     1
#> 15 30 0.4674458 0.7721464 0.1619513 0.7994688     1
#> 16 32 0.4555030 0.8610687 0.1339908 0.8189431     1
#> 17 34 0.4436111 0.9514457 0.1099435 0.8373008     1

^{Created on 2019-05-06 by the reprex package (v0.2.1)}

There are some vignettes that show the different features of the package, this one here demonstrates how to compute marginal effects at specific values / levels of focal terms.

The recipe posted by Ben that shows how to calculate the confidence intervals (conditioned or not conditioned on random effects) is implemented in ggpredict(), a short vignette explaining the differences is here.

Answer 2

The easiest way to do this is probably with ggeffects::ggpredict().

Something like

library(ggeffects)
pred_frame <- ggpredict(myModel, term="evi_prev")

should produce a data frame with predictions, lower and upper confidence levels. I'm not sure whether it will make the predictions for evenly spaced values along the x-axis (which would be nice), or how to trick it into doing so. (If you provide a reproducible example I might give it a shot.)

Playing around with the MWE you posted does suggest that it's hard to get predictions for evenly spaced values (or more generally, for values that aren't in the original data); I tried things like terms="size [1:35]", but this restricts the range of the predicted values rather than filling them in.

More basically, the built-in predict() method for merMod objects can be used (possibly with newdata to specify e.g. evenly spaced values) to get predictions [use type="response" to get predictions on the probability rather than the log-odds scale]; confidence intervals are harder but can be generated with the recipe shown here