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