incorrect labelling of variance components in summary() for 3+ trait models fitted with AI-REML

[famuvie]

Nobs <- 1e4

## residual covariance matrix
S_res <- matrix(c(
  9, 3, -3,
  3, 9, 9,
  -3, 9, 14
), nrow = 3, ncol = 3)

## simulated residual-only dataset
testdat <- data.frame(breedR.sample.ranef(3, S_res, Nobs, vname = 'e'))

## fitted model
res <- remlf90(
  cbind(e_1, e_2, e_3) ~ 1,
  data = testdat,
  method = "ai"
)  
#> Using default initial variances given by default_initial_variance()
#> See ?breedR.getOption.

summary(res)
#> Formula: cbind(e_1, e_2, e_3) ~ 0 + Intercept 
#>    Data: testdat 
#>    AIC   BIC logLik
#>  83966 84031 -41974
#> 
#> 
#> Variance components:
#>                           Estimated variances    S.E.
#> Residual.e_1                            9.021 0.12758
#> Residual.e_1_Residual.e_2               3.098 0.09474
#> Residual.e_2                           -2.900 0.11476
#> Residual.e_1_Residual.e_3               8.886 0.12567
#> Residual.e_2_Residual.e_3               8.787 0.14095
#> Residual.e_3                           13.666 0.19327
#> 
#> Fixed effects:
#>                   value   s.e.
#> Intercept.e_1 -0.042498 0.0300
#> Intercept.e_2  0.020878 0.0298
#> Intercept.e_3  0.045872 0.0370

res$var[["Residual", "Estimated variances"]]
#>         e_1    e_2     e_3
#> e_1  9.0205 3.0978 -2.9004
#> e_2  3.0978 8.8856  8.7875
#> e_3 -2.9004 8.7875 13.6660

Notice how the value summarized as Residual.e_2 corresponds in fact to the residual covariance between e_2 and e_3. In particular, it is negative!. The values are fine, but the labeling is incorrect.

[famuvie]

Thanks to Vincent Segura for reporting.