Created attachment 2035 [details]
a "fixed" version.
This probably doesn't fix the bug in all cases but at least is strictly better than the current code.
There is a bug in `getVarCov` which rears its head when you use the `weights` argument to either lme or gls. It has to do with the order of the rows of the data set, something which obviously should not matter; and it means that you get simply wrong numbers for marginal variances, though correlations should still be correct. I've written a fix for this function that works in at least some additional cases (attached).