Formula with multiple Y and clustered errors: errors in computation?

[MatthieuStigler]

If one computes a formula like y1+y2~x1|fe|0|fe, the standard errors will be correct (i.e. the same as from the formula with only one response) only for the first response, not the second!?

See:

library(lfe)
#> Loading required package: Matrix
library(plm)
#> 
#> Attaching package: 'plm'
#> The following object is masked from 'package:lfe':
#> 
#>     sargan

data(Produc)

coef(summary(felm(pcap~hwy|state|0|state, data=Produc)))
#>     Estimate Cluster s.e.  t value    Pr(>|t|)
#> hwy  3.05037    0.3654909 8.345954 7.79537e-11
coef(summary(felm(pcap+util~hwy|state|0|state, data=Produc), lhs="pcap"))
#>     Estimate Cluster s.e.  t value    Pr(>|t|)
#> hwy  3.05037    0.3654909 8.345954 7.79537e-11
coef(summary(felm(util+pcap~hwy|state|0|state, data=Produc), lhs="pcap"))
#>     Estimate Cluster s.e. t value   Pr(>|t|)
#> hwy  3.05037     1.075536 2.83614 0.00671554

coef(summary(felm(util~hwy|state|0|state, data=Produc)))
#>     Estimate Cluster s.e.  t value     Pr(>|t|)
#> hwy 1.522517     0.236826 6.428842 6.021228e-08
coef(summary(felm(pcap+util~hwy|state|0|state, data=Produc), lhs="util"))
#>     Estimate Cluster s.e.  t value   Pr(>|t|)
#> hwy 1.522517    0.6969115 2.184663 0.03393527
coef(summary(felm(util+pcap~hwy|state|0|state, data=Produc), lhs="util"))
#>     Estimate Cluster s.e.  t value     Pr(>|t|)
#> hwy 1.522517     0.236826 6.428842 6.021228e-08



coef(summary(felm(pcap~hwy|state, data=Produc)))
#>     Estimate Std. Error  t value      Pr(>|t|)
#> hwy  3.05037 0.08476667 35.98549 7.128369e-167
coef(summary(felm(pcap+util~hwy|state, data=Produc), lhs="pcap"))
#>     Estimate Std. Error  t value      Pr(>|t|)
#> hwy  3.05037 0.08476667 35.98549 7.128369e-167
coef(summary(felm(util+pcap~hwy|state, data=Produc), lhs="pcap"))
#>     Estimate Std. Error  t value      Pr(>|t|)
#> hwy  3.05037 0.08476667 35.98549 7.128369e-167

^{Created on 2021-01-15 by the reprex package (v0.3.0)}