I want to plot a similar plot as this one in the buttom of the page: ordered logit
They use a variabel on the x-axis that is categorical (0-10) and therefore they use seq(0, 10, 1) and hold all other variables constant at their means when they determine the predicted probabilities for the plot.
However in my data the variabel I want to have on the x-axis is a dummy (0-1) and therefore seq() does not work. What can I do instead to do a similar plot? Do you have any exapmples.
It's difficult to walk you through your own code in the absence of a reproducible example, but let's create one with a dummy variable on the x axis.
Here, we have a binary outcome variable, a dummy predictor variable, and a numeric predictor variable that we want to hold at its mean value for the purposes of the plot:
set.seed(1)
df <- data.frame(dummy = rep(0:1, each = 20),
other_var = runif(40),
outcome = rbinom(40, 1, rep(c(0.3, 0.7), each = 20)))
We can create a logistic regression like this:
model <- glm(outcome ~ dummy + other_var, family = binomial, data = df)
summary(model)
#>
#> Call:
#> glm(formula = outcome ~ dummy + other_var, family = binomial,
#> data = df)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -1.4942 -0.9675 -0.5557 0.9465 2.0245
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.006877 0.831579 -0.008 0.9934
#> dummy 1.447983 0.713932 2.028 0.0425 *
#> other_var -2.120011 1.339460 -1.583 0.1135
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 54.548 on 39 degrees of freedom
#> Residual deviance: 46.726 on 37 degrees of freedom
#> AIC: 52.726
#>
#> Number of Fisher Scoring iterations: 4
To predict the outcome for each value of dummy at the mean of other_var, we create a little prediction data frame with one column holding each unique value of dummy and another column filled with the mean value of other_var:
pred_df <- data.frame(dummy = 0:1, other_var = mean(df$other_var))
Now we feed this to predict
preds <- predict(model, pred_df, se.fit = TRUE)
Since we are dealing with log odds here, we need a little function to convert log odds to probabilities:
log_odds_to_probs <- function(x) exp(x) / (1 + exp(x))
Now we can get our predicted values as well as the 95% confidence intervals for the predictions as probabilities:
pred_df$fit <- log_odds_to_probs(preds$fit)
pred_df$lower <- log_odds_to_probs(preds$fit - 1.96 * preds$se.fit)
pred_df$upper <- log_odds_to_probs(preds$fit + 1.96 * preds$se.fit)
Finally, we plot the result using whatever style you like:
library(ggplot2)
ggplot(pred_df, aes(factor(dummy), fit)) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, size = 1.5,
color = 'deepskyblue4') +
geom_point(size = 4) +
labs(x = 'dummy', y = 'probability') +
ylim(c(0, 1)) +
theme_minimal(base_size = 16)
Created on 2022-12-06 with reprex v2.0.2
I want to offer a much easier solution. For comparability, I also start with the model proposed by Allan.
model <- glm(outcome ~ dummy + other_var, family = binomial, data = df)
Now to calculate the predicted probabilities, use the glm.predict package by Benjamin Schlegel. This package does all the hard stuff for you with just one command. You can even set the method (boostrap or simulation) and it calculates discrete changes if you like. Also, the package fits within a tidyverse pipeline.
You set you dummy to 0 and 1 and your other variable to the mean. From that you can directly calculate the same plot as Allan has. Below you see the code for that.
library(glm.predict)
predicts(model, "0,1;mean") %>%
ggplot(aes(x=as.factor(dummy), y=mean, ymin=lower, ymax=upper))+
geom_pointrange()
