I have three outcomes of interest, each of which is measured on multiple occasions (before, during, and after an intervention). Given my data structure (and because I’m not necessarily interested in random effects), I’m attempting to apply a cluster-robust approach to multivariate regression. My code looks something like this:

```
library(estimatr)
y <- cbind(a, b, c)
lm.robust <- lm_robust(y~factor(Timepoint), clusters = ID, data = outcome.data)
summary(lm.robust)
```

As expected, when I run this code, the summary breaks out each outcome variable separately. Based on past experience using multivariate models using lm(), I would expect this portion of the output to be identical to running three separate models, each with variables “a”, “b”, and “c” as the dependent variables. Instead, however, for every outcome variable other than the first one, the multivariate model calculates an estimate, a standard error, and a t value, but returns NA in lieu of p-values and confidence intervals, and the r-squared values that it returns are ludicrously high (e.g., adjusted r2 = .97). This is true regardless of the order of my dependent variables (i.e., if the order is abc, the summary returns correct coefficients for a, but not b or c; if the order is bca, then the summary returns correct coefficients for b, but not c or a).

I’m trying to understand why this might be the case. It may be worth noting that there is a fair amount of missingness in my outcome variables. Any insight is much appreciated!

P.S. I’m also interested in calculating multivariate test statistics based on the cluster-robust SEs. Will this code do the trick?

`summary(manova(lm.robust.test))`