Combining estimates from multivariate multiple regression using MICE in R? A similar question was posted last year, but it didn't get any love, so!:
I'm trying to calculate a pooled estimate after fitting a multivariate regression in multiply imputed data (having used MICE to do the imputation). Specifically, I fit a weighted, multivariate linear regression model with 4 outcomes. Unfortunately, pool() doesn't seem able to combine estimates from multivariate regression (Error: Glance does not support linear models with multiple responses.), and I can't find (A) any other packages that do or (B) any resources for how to combine estimates from multivariate regression using Rubin's rules so that I could do it manually! (Is it the same as if I was just combining the output across imputations for a typical linear regression? I fear not...)
Does anyone have any ideas/resources they could point me to? I really only care about getting pooled estimates for the relationship between my main exposure variable (let's call it X1) and the 4 outcomes + the standard errors, p-values, CIs for those relationships; plus, would like to be able to run a joint test of whether all of the X1-related betas are equal to 0 (beta_x1_outcome1 = beta_x1_outcome2 = beta_x1_outcome3 = beta_x1_outcome4 = 0). I have combined estimates using Rubin's rules before, but only for a model with one outcome.
Any help would be so appreciated!
PS - I'm using R 3.4; here's a generic version of the model:
fit <- with(data=imp, lm(cbind(y1, y2, y3, y4) ~ x1 + x2 + x3 + x4, weight=weight, subset=(analysis_indcator==1)))
summary(pool(fit))
 A: There are two separate issues here. First is the problem of getting coefficient covariance matrices from a weighted multivariate (in the sense of multiple outcomes) regression. The second is how to apply Rubin's rules to put together the results of multivariate regressions on a group of imputed data sets.
The first problem is the most vexing. As noted in a comment on the question, although the R stats package will fit a weighted multivariate regression, it does not support tools for multivariate inference on the resulting "mlm" objects that are available for unweighted regressions.
I suspect that's because it's not completely clear what a weight in a multivariate analysis corresponds to. The help page for lm() says:

Non-NULL weights can be used to indicate that different observations have different variances (with the values in weights being inversely proportional to the variances)...

That's easy to understand when there's a single outcome. What does that mean when there are multiple outcomes? Do all of the outcomes have the same relative variances from case to case? Does that even make sense when you're trying to model the covariances among the error terms? At the least, it seems that you are putting some strong restrictions on those covariances.
Another use of weights is for unit-weighting of observations: for example, when outcome values have been summarized over different numbers of observations for each set of covariates. But in that case, the help page warns:

... notice that within-group variation is not used. Therefore, the sigma estimate and residual degrees of freedom may be suboptimal; in the case of replication weights, even wrong. Hence, standard errors and analysis of variance tables should be treated with care.

So I offer the following work-around with a good deal of trepidation: instead of providing weights directly to lm(), do the equivalent by weighting the data. Form the diagonal matrix whose elements are the square roots of the weights, and premultiply both the design matrix and the response matrix by it. Then do an unweighted multivariate regression on the weighted data.
As the Wikipedia page for weighted least-squares regression shows, that provides the same coefficient estimates as a regression with weights. The coefficient estimates for a multivariate regression are the same as those from individual regressions for each of the responses. Thus you will effectively have done a weighted regression of the type that you specified, but lm() and the downstream functions won't know that. Whether the results of that weighting make sense is another question.
Once an unweighted "mlm" object is obtained for an imputed data set, then you get the matrix of regression coefficients and their covariance matrix with coef() and vcov(), respectively. As there now is a whole error-covariance matrix instead of the single error variance estimate obtained in ordinary least squares, and there are intercepts and regression coefficients for each of the outcomes, the covariance matrix among the coefficient estimates will be large. With 4 outcomes, 4 predictors, and intercepts as in this question it will be 20 x 20.
With all of the imputed data sets analyzed, Rubin's rules can then be applied. Marshall et al. provide a useful (and freely available) summary. You use the coefficient matrices to get (1) the average coefficient matrix, and (2) the between-imputation variance of coefficient estimates, called $B$. You average the coefficient-covariance matrices (as large as those might be) for an estimate of within-imputation variance, called $\bar U$.
Table 2 of that paper provides formulas for combining the results on the individual imputations for scalar or  multivariate tests when multivariate normality is assumed to hold, for $\chi^2$ statistics, and for likelihood-ratio tests. For multivariate tests like the desired "joint test of whether all of the X1-related betas are equal to 0"  this will be a Wald test producing an $F$ statistic with $k$ and $\nu$ degrees of freedom. Here, $k$ is the number of coefficients being tested (4 for the joint test specified in the question) and $\nu$ is a value corrected for the relative increase in variance due to the imputation. The correction is a matrix generalization of the method shown for scalar t-tests in another question.
A: have you figured out a workable solution yet? unfortunately it seems there is currently no support based on mice's dependency on broom here. ideally, one would come up with a method for broom's glance method. you can find the source code of how mice pools datasets here, right at the bootom where 'pool' is defined.
regarding just estimates: estimates of univarte and multivarite analysis should be the same. so you could either loop over your outcome variables and run univariate analysis, or use this dirty workaround i came up with this simple script. however it does nothing more than separating out the multivariate models into univariate models and then pooling each model using the mice::pool method. disclaimer: i am neither a biostatitician nor good at programming, this does not account for the multivariate nature but only gives you the pooled estimates for each outcome variable.
pool_multi <- function(object) {
    m <- length(object$analyses)
    if (length(dimnames(object$analyses[[1]]$coefficients)[[2]]) > 1) {
        names <- dimnames(object$analyses[[1]]$coefficients)[[2]]
        pools <- list()
        analysis <- list()
        analyses <- list()
        for (name in names) {
            for (k in seq(1,m)) {
                coefficients <- object$analyses[[k]]$coefficients[, name]
                fitted.values <- object$analyses[[k]]$fitted.values[, name]
                effects <- object$analyses[[k]]$effects[, name]
                residuals <- object$analyses[[k]]$residuals[,name]
                rank <- object$analyses[[k]]$rank
                assign <- object$analyses[[k]]$assign
                qr <- object$analyses[[k]]$qr
                df.residual <- object$analyses[[k]]$df.residual
                contrasts <- object$analyses[[k]]$contrasts
                xlevels <- object$analyses[[k]]$xlevels
                call <- object$analyses[[k]]$call
                terms <- object$analyses[[k]]$terms
                model <- object$analyses[[k]]$model #[, name]
                analysis <- list(residuals = residuals, effects = effects, rank = rank,
                                 fitted.values = fitted.values, assign = assign,
                                 qr = qr, df.residual = df.residual, contrasts = contrasts,
                                 xlevels = xlevels, call = call, terms = terms, model = model,
                                 coefficients = coefficients
                                )
                class(analysis) <- c('lm')
                analyses[[k]] <- analysis
            }
            poolobj <- list(call = object$call, call1 = object$call1,
                            nmis = object$nmis, analyses = analyses)
#            poolobj <- new('mira', call = object$call, call1 = object$call1,
#                           nmis = object$nmis, analyses = analyses[1])
            class(poolobj) <- c("mira", "matrix")
            pools[[name]] <- mice::pool(poolobj)
        }
        pools
    }
    else {
        print('Not a multivariate linear regression object. Use mice::pool instead.')
    }
}

