How do I extract standard errors or variation from predicted ordinal logistic regression analyses?

user_11468006 0
Viewed 42

I am undertaking a ordinal logistic regression using R package MASS.

For example:

library(MASS)
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
summary(house.plr, digits = 3)

I am using the s3 method predict() to obtain the predicted values

test_dat <- data.frame(Infl = factor(rep("Low",4)),
                       Cont = factor(rep("Low",4)),
                       Type = unique(housing$Type))
predict(house.plr, test_dat, type = "p")

        Low    Medium      High
1 0.3784493 0.2876752 0.3338755
2 0.5190445 0.2605077 0.2204478
3 0.4675584 0.2745383 0.2579033
4 0.6444840 0.2114256 0.1440905

The result is a table of predicted means for each level of Sat given the variables defined in the test_dat.

How might I extract the variation around each of these means in the form of a standard error or standard deviation?

predict regression r
1 Answers

First, your predicted values are the predicted probability of each outcome for each observation. It is not the predicted mean on the response scale.

Second, you can use the marginaleffects package to get the standard errors for the predicted probabilities and then calculate the confidence intervals yourself. Alternatively, you may implement the non-parametric bootstrap. I implement both in the below. Note that I shifted the order of the columns around in the test data to match the training data.

# Packages
library(MASS)
library(marginaleffects)
library(dplyr)

# Create a test set
N <- 4
test_dat <- data.frame(
  Infl = factor(rep("Low", N)),
  Type = unique(housing$Type),
  Cont = factor(rep("Low", N))
)

# Fit ordered logistic regression model
house.plr <- polr(Sat ~ Infl + Type + Cont, 
                  weights = Freq, 
                  data = housing,
                  Hess = TRUE)

# Demonstrate that predict() doesn't provide any measure of variability
# for the predicted class probabilities, as shown in question
predict(house.plr, test_dat, type = "probs")

# Use the marginaleffects package to get delta method standard errors for
# each predicted probability
probs <- marginaleffects::predictions(house.plr, 
                                      newdata = test_dat,
                                      type = "probs")

# Compute CIs from the standard error using normal approximation
probs$predicted - 1.96*probs$std.error
probs$predicted + 1.96*probs$std.error

# Alternatively, use non-parametric bootstrapped confidence intervals.
# note that this does not adjust the weights to a constant sum for
# each bootstrap, although it is easy to implement. You're free to
# determine how to handle the weights, including resampling based 
# on the weights.

# Generate bootstrapped data.frames
set.seed(123)
sims <- 5
samples <- vector(mode = "list", length = sims)
samples <- lapply(samples, function(x){ slice_sample(housing, n = nrow(housing), replace = TRUE)})

# Fit model on each bootstrapped data.frame
models <- lapply(samples, function(x){polr(Sat ~ Infl + Type + Cont, 
                                           weights = Freq, 
                                           data = x,
                                           Hess = TRUE)})

# Get test predictions into a data.frame
probs_boot <- lapply(models, function(x) {
  marginaleffects::predictions(x,
                               newdata = test_dat,
                               type = "probs")
})

probs_boot_df <- bind_rows(probs_boot)

# Compute CIs
probs_boot_df %>% 
  group_by(group, Type.x, Infl, Type.y, Cont) %>%
  summarise(ci_low = quantile(predicted, probs = 0.025),
            ci_high = quantile(predicted, probs = 0.975))
user_36146480
Related