2
$\begingroup$

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.

enter image description here

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()

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Robert Long Jul 1 at 17:40
  • $\begingroup$ 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. $\endgroup$ – Mark Miller Jul 1 at 18:27
  • 1
    $\begingroup$ 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 ? $\endgroup$ – Robert Long Jul 1 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.