How to conduct joint significance test in seemingly unrelated regression

Question

I'm trying to conduct a joint test of significance in a seemingly unrelated regression setup with robust standard errors. I have three outcomes Y1, Y2, and Y3 and I want to conduct a joint hypothesis test against the null that the average effect of the treatment Z is zero on all three outcomes.

I think that I have the model set up correctly, but I don't think that I have the hypothesis.matrix set correctly in car::linearHypothesis.

Here's some data:

library(tibble)
library(car)
library(systemfit)

set.seed(343)
N = 800
dat <-
  tibble(
    U = rnorm(N),
    Z = rbinom(N, 1, 0.5),
    Y = 0.2 * Z + U,
    Y1 = Y + rnorm(N, sd = 0.3),
    Y2 = Y + rnorm(N, sd = 0.5),
    Y3 = Y + rnorm(N, sd = 0.5)
  )

Here's the seemingly unrelated regression fit:

sur <- systemfit(list(Y1 ~ Z, Y2 ~ Z, Y3 ~ Z), method = "SUR", data = dat)
summary(sur)

Which is identical to the ols fit in this case:

ols <- lm(cbind(Y1, Y2, Y3) ~ Z, data = dat)
summary(ols)

Which is useful, because I need to estimate robust standard errors for this test:

linearHypothesis(ols, hypothesis.matrix = "Z = 0", white.adjust = "hc2")

This last line is the one that I think is incorrect. I think it's incorrect because the individual coefficients all have lower p-values than the joint test, but I could be wrong?

Answer 1

Looks right to me. You'd get the same result by estimating the null model (ols0 below) and using anova() to test the difference between the estimated and null models.

library(tibble)
library(car)
#> Loading required package: carData

set.seed(343)
N = 800
dat <-
  tibble(
    U = rnorm(N),
    Z = rbinom(N, 1, 0.5),
    Y = 0.2 * Z + U,
    Y1 = Y + rnorm(N, sd = 0.3),
    Y2 = Y + rnorm(N, sd = 0.5),
    Y3 = Y + rnorm(N, sd = 0.5)
  )


ols <- lm(cbind(Y1, Y2, Y3) ~ Z, data = dat)

linearHypothesis(ols, hypothesis.matrix = "Z = 0")
#> 
#> Sum of squares and products for the hypothesis:
#>          Y1       Y2       Y3
#> Y1 3.201796 4.693391 3.359617
#> Y2 4.693391 6.879863 4.924734
#> Y3 3.359617 4.924734 3.525216
#> 
#> Sum of squares and products for error:
#>          Y1       Y2       Y3
#> Y1 829.5535 756.1586 770.0808
#> Y2 756.1586 965.5959 770.4636
#> Y3 770.0808 770.4636 980.0664
#> 
#> Multivariate Tests: 
#>                  Df test stat approx F num Df den Df  Pr(>F)
#> Pillai            1 0.0073689  1.96972      3    796 0.11703
#> Wilks             1 0.9926311  1.96972      3    796 0.11703
#> Hotelling-Lawley  1 0.0074236  1.96972      3    796 0.11703
#> Roy               1 0.0074236  1.96972      3    796 0.11703

ols0 <- lm(cbind(Y1, Y2, Y3) ~ 1, data = dat)
anova(ols, ols0, test="Pillai")
#> Analysis of Variance Table
#> 
#> Model 1: cbind(Y1, Y2, Y3) ~ Z
#> Model 2: cbind(Y1, Y2, Y3) ~ 1
#>   Res.Df Df Gen.var.    Pillai approx F num Df den Df Pr(>F)
#> 1    798     0.48198                                        
#> 2    799  1  0.48257 0.0073689   1.9697      3    796  0.117

^{Created on 2022-07-08 by the reprex package (v2.0.1)}