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
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
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.