# plot probabilities from mixed effect logistic regression with random intercept and fixed continuous and categorical covariate

I have created a mixed-effect logistic regression model with a random intercept, a fixed continuous covariate and a fixed categorical covariate. I want to plot the probability of success (y = 1) vs the continuous covariate for each level of the categorical covariate.

When creating the plot should I just use the fixed intercept; or should I use the mean value of the fixed + continuous intercept for each observation in each level of the categorical covariate; or should I use a different approach?

Below I have provided functional and reproducible R code that generates a fake data set for such a model and creates plots using both of the above approaches. The plots are similar but are not identical. Sorry the code is so long, but the model is a little complex.

library(arm)

# Generate data
set.seed (1234)
n.groups <- 10                                          # number of groups
n        <- 100                                         # number of observations
group    <- round(runif(n, 1, n.groups))                # group
mu.a     <- 5                                           # mean used for random group intercept
sigma.a  <- 2                                           # SD used for random group intercept
B0       <- rnorm(n.groups, mu.a, sigma.a)              # random group intercepts
B1       <- -1.5                                        # slope on continuous covariate
x1       <- round(runif(n, -6 , 6), digits = 3)         # continuous covariate x1
B2.2     <- 0.5                                         # slope on level 2 of categorical covariate x2
B2.3     <- 2.5                                         # slope on level 3 of categorical covariate x2
x2       <- round(runif(n, 1, 3))                       # categorical covariate x2
x2.2     <- ifelse(x2 == 2, 1, 0)                       # indicator for level 2 of categorical covariate x2
x2.3     <- ifelse(x2 == 3, 1, 0)                       # indicator for level 3 of categorical covariate x2
lin.pred <- (B0[group] + B1*x1 + B2.2*x2.2 + B2.3*x2.3) # linear predictor
e        <- exp(1)                                      # e
logit.p  <- e^(lin.pred) / (1 + e^(lin.pred))           # linear predictor with logit link
y        <- rbinom(n, size=1, prob=logit.p)             # simulate random outcome of model

# Analyze data
group    <- as.factor(group)
x2       <- as.factor(x2)
my.model <- glmer (y ~ x1 + x2 + (1 | group), family=binomial(link="logit"))
summary(my.model)

# Extract estimated coefficients
fixed.effects      <- fixef(my.model)
random.effects     <- ranef(my.model)
my.coef            <- coef(my.model)
fitted.values      <- fitted(my.model)
fixed.int          <- as.numeric(fixed.effects[1])
x1.effect          <- as.numeric(fixed.effects[2])
x22.effect         <- as.numeric(fixed.effects[3])
x23.effect         <- as.numeric(fixed.effects[4])

# Create data set with observed values of y, group, x1, x2 and estimated coefficients
x2                 <- as.numeric(as.character(x2))
my.coef2           <- data.frame(group = 1:n.groups, my.coef$group) colnames(my.coef2) <- c('group', 'intercept', 'B.x1', 'B.x22', 'B.x23') my.data <- data.frame(y, group, x1, x2) my.data2 <- merge(my.data, my.coef2, by = 'group', all = TRUE) # Create plots png('cross_validated_logistic.png', width = 1080, height = 720) par(mfrow=c(1,2)) # Do I create a plot using the fixed intercept? plot(x1, y, ylab="Probability of Success", xlab="x1", cex.lab=1.5, pch = 19, bty = "l", cex.main = 1.5, cex = 1.5, lwd = 2) curve (invlogit(fixed.int + x1.effect*x ), min(x1), max(x1), lwd=2, lty=1, add=TRUE) curve (invlogit(fixed.int + x1.effect*x + x22.effect), min(x1), max(x1), lwd=2, lty=2, add=TRUE) curve (invlogit(fixed.int + x1.effect*x + x23.effect), min(x1), max(x1), lwd=2, lty=9, add=TRUE) title('fixed intercepts') # Do I create a plot using the means of the (fixed+random) intercepts for all observations within # each level of x2? plot(x1, y, ylab="Probability of Success", xlab="x1", cex.lab=1.5, pch = 19, bty = "l", cex.main = 1.5, cex = 1.5, lwd = 2) curve (invlogit(mean(my.data2$$intercept[my.data2$$x2 == 1]) + x1.effect*x ), min(x1), max(x1), lwd=2, lty=1, add=TRUE) curve (invlogit(mean(my.data2$$intercept[my.data2$$x2 == 2]) + x1.effect*x + x22.effect), min(x1), max(x1), lwd=2, lty=2, add=TRUE) curve (invlogit(mean(my.data2$$intercept[my.data2$$x2 == 3]) + x1.effect*x + x23.effect), min(x1), max(x1), lwd=2, lty=9, add=TRUE) title('mean fixed+random intercepts') dev.off()  EDIT: July 1, 2020 The predict function uses the fixed intercept + the random intercept for each observation in the data set: as.numeric(predict(my.model)) [1] -0.1355467 -2.1389760 -5.2377738 0.6662584 10.2387598 17.5963610 9.5228168 3.7000504 12.1100313 12.1377452 3.1985063 4.0260289 13.3135166 [14] 18.3857911 -3.4260311 10.7098732 -3.2675657 1.2716749 0.3979737 -2.6588965 7.1504649 -8.6896208 7.1710257 0.5748609 10.5022989 1.7939458 [27] 4.3076396 11.6691996 -4.2395046 1.1648367 8.6207565 10.7581858 17.9582347 -1.5288265 7.1199373 -7.4741397 1.1646037 -12.5147846 10.2446393 [40] 3.6527023 12.3294160 6.6001777 8.9556858 19.1039239 14.4035017 0.4605888 13.8060396 14.5525344 -7.4695883 9.5608927 -0.4696090 8.1002000 [53] 6.6665300 12.3915173 6.6784981 0.8138920 12.0524027 3.5049440 13.6469675 18.4275355 13.9091297 4.0825143 7.6761979 14.0326569 -7.9343093 [66] 24.3118685 -1.5608503 4.0688574 4.9754884 4.2819425 2.5712139 0.7093886 12.6921251 -8.2600107 -1.8272722 2.2454477 -1.3788292 -3.3219547 [79] -5.1927063 8.4567610 9.2343998 9.3713392 10.8306885 0.5890424 6.1509242 7.3949161 8.5158585 12.7278185 13.7366016 3.9290338 7.4804431 [92] 19.0078616 14.2334121 -0.4910232 0.9437631 -7.1157713 0.1758769 -5.0635100 7.7543466 -0.8881751 > as.numeric(unlist(my.coef$group[1]))[as.numeric(group)] + x1.effect*x1 + x22.effect*x2.2 + x23.effect*x2.3
[1]  -0.1355467  -2.1389760  -5.2377738   0.6662584  10.2387598  17.5963610   9.5228168   3.7000504  12.1100313  12.1377452   3.1985063   4.0260289  13.3135166
[14]  18.3857911  -3.4260311  10.7098732  -3.2675657   1.2716749   0.3979737  -2.6588965   7.1504649  -8.6896208   7.1710257   0.5748609  10.5022989   1.7939458
[27]   4.3076396  11.6691996  -4.2395046   1.1648367   8.6207565  10.7581858  17.9582347  -1.5288265   7.1199373  -7.4741397   1.1646037 -12.5147846  10.2446393
[40]   3.6527023  12.3294160   6.6001777   8.9556858  19.1039239  14.4035017   0.4605888  13.8060396  14.5525344  -7.4695883   9.5608927  -0.4696090   8.1002000
[53]   6.6665300  12.3915173   6.6784981   0.8138920  12.0524027   3.5049440  13.6469675  18.4275355  13.9091297   4.0825143   7.6761979  14.0326569  -7.9343093
[66]  24.3118685  -1.5608503   4.0688574   4.9754884   4.2819425   2.5712139   0.7093886  12.6921251  -8.2600107  -1.8272722   2.2454477  -1.3788292  -3.3219547
[79]  -5.1927063   8.4567610   9.2343998   9.3713392  10.8306885   0.5890424   6.1509242   7.3949161   8.5158585  12.7278185  13.7366016   3.9290338   7.4804431
[92]  19.0078616  14.2334121  -0.4910232   0.9437631  -7.1157713   0.1758769  -5.0635100   7.7543466  -0.8881751


I might now try to express my question using mathematical symbols in addition to the R code above.

Here I plot the predicted values of probability of success for each of the 100 observations in the data set across a range of x1. I color code the observations by their value of x2. The pattern generally appears to be what I would expect: red lines occur everywhere; green lines are mostly to the right of blue lines. But this plot seems too messy, which is why I was hoping to have just one line for each level of x2. I just am not sure of the best way to create such a plot.

x1.for.expected <- seq(min(x1), max(x1), length = length(unique(x1)))

expected.matrix <- matrix(0, nrow=nrow(my.data2), ncol=length(unique(x1)))

for(i in 1:length(x1.for.expected)) {

expected.matrix[i,] <- invlogit(my.data2$$intercept[i] + x1.effect*x1.for.expected + x22.effect * my.data2$$x2.2[i] + x23.effect * my.data2\$x2.3[i])

}

png('cross_validated_logistic_July1_2020.png', width = 1080, height = 720)

plot(x1, y, ylab="Probability of Success", xlab="x1",
cex.lab=1.5, pch = 19, bty = "l", cex.main = 1.5, cex = 1.5, lwd = 2)

for(i in 1:n) {

if(x2[i] == 1) lines(x1.for.expected, expected.matrix[i,], col = "green", lwd = 2)
if(x2[i] == 2) lines(x1.for.expected, expected.matrix[i,], col = "red",   lwd = 2)
if(x2[i] == 3) lines(x1.for.expected, expected.matrix[i,], col = "blue",  lwd = 2)

}

title('expected values for each observation')

dev.off()


• I'm a little confused. Do you want to plot predicted probabilities ? It might be useful if you write down the model (maths) and then explain what your two approaches mean in terms of the maths. – Robert Long Jul 1 at 17:40
• I my mind I was trying to plot the expected values, but instead of plotting them for each observation I wanted to plot them for each level of x2. A third approach I did not show would be to obtain the expected values for each observation and then average each point of every observation that fell within the same level of x2. Maybe I will write the R code for that third approach and add it to my question. – Mark Miller Jul 1 at 18:27
• Please also write out the equation(s) for the model, it makes these kinds of questions much easier to answer. Is there any reason you didn't use the predict function ? Also, while your simulation code is very nice, there is probably no need to actually deal with the inverse logit etc to plot probabilities in order to answer your question - you would have the same question if it was just a linear model that you wanted to plot expeced values for, right ? – Robert Long Jul 1 at 18:32