Please specify whether your issue is about:

• a possible bug
• a question about package functionality
• a suggested code or documentation change, improvement to the code, or feature request

Hi @leeper. Thanks again for your work on this. In general, I think expecting "coefficients" to be in the output of a model object should be deprecated in favor of using s3 methods like "coef". You have done this in some places in your package (since the last CRAN submission).

However, in a crucial location you still overrwite model[["coefficients"]] and then expect coef(model) to extract that object. For example, in get_effect_variances.R, see the following lines:

FUN <- gradient_factory(data = data,
                                model = model,
                                variables = variables,
                                type = type,
                                weights = weights,
                                eps = eps,
                                varslist = varslist,
                                ...)
# get jacobian
jacobian <- jacobian(FUN, coef(model)[names(coef(model)) %in% c("(Intercept)", colnames(vcov))], weights = weights, eps = eps)

Now, in gradient_factory.default, you do the following:

model[["coefficients"]][names(coefs)] <- coefs

So if while in the gradient factory you think you are clobbering model[["coefficients"]], the next line in get_effect_variances uses coef(model) and won't get those new coefficients if coef.model_class != model$coefficients.