I am estimating a series of models with varying combinations of outcomes and treatments and I am seeking to estimate the ATT via average marginal effects using the {marginaleffects} package. Below, I will provide syntax for the estimation of the marginal effects for exposed units:

ame.robust <- avg_comparisons(model, variables = treatment,
                              newdata = subset(final, treatment == 1),
                              vcov = ~id)

However, when I run this code, I get the following error:

Error: Unable to compute predicted values with this model. This error can arise when `insight::get_data()` is unable to extract
  the dataset from the model object, or when the data frame was modified since fitting the model. You can try to supply a
  different dataset to the `newdata` argument.

I have referred to this query and removed NAs from the final data set. But I still get the same error. Further, all of the data in the loop (which I will provide) are all stored in the final data set, so I don’t see how this is an issue. Importantly, the entire loop ran just fine until I supplied the newdata() argument.

Here is the syntax for the loop:

outcomes <- c("y1", "y2", "y3", "y4")

treatments <- c("x1", "x2", "x3", "x4")

adjustment <- c("z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z9 + z10")

model.list <- list()
ame.list <- list()
results.list <- list()
or.list <- list()
  
  for (outcome in outcomes) {
    
    model.list[[outcome]] <- list()
    
    for (treatment in treatments) {
      
      # Fit the Model
      model <- polr(as.formula(paste(outcome, "~", treatment, "+", adjustment)), data = final)
      
      # Compute Clustered Robust Standard Errors 
      model.robust <- coeftest(model, vcov = vcovCL, type = "HC1", df = 2400, cluster = ~id)
      
      # Store Results
      model.name <- paste(outcome, treatment, sep = ".")
      model.list[[model.name]] <- model.robust
      
      # Calculate/Store Odds Ratios for Sensitivity Analysis
      treatment.coef <- coef(model)[[treatment]]
      
      odds.ratio <- exp(treatment.coef)
      
      # Calculate Average Marginal Effects
      ame.robust <- avg_comparisons(model,
                                    variables = treatment,
                                    newdata = subset(final, treatment == 1),
                                    vcov = ~id)
      
      # Store Results
      ame.name <- paste(outcome, treatment, sep = ".")
      ame.list[[ame.name]] <- ame.robust
      
      # Store Odds Ratio
      or.list[[ame.name]] <- odds.ratio
      
      # Manually Extract Estimates and Confidence Intervals
      estimate <- ame.robust$estimate
      conf.low <- ame.robust$conf.low
      conf.high <- ame.robust$conf.high
      
      new.row <- data.frame(outcome = outcome,
                            treatment = treatment,
                            estimate = estimate,
                            conf.low = conf.low,
                            conf.high = conf.high)
      
      results.list <- rbind(results.list, new.row)
      
    }
    
  }