# Why did I get inconsistent result from ANOVA and wald test

I tested the effect of one factor X on my response Y. The factor has four levels A, B, C, and D and the null hypothesis I want to test is A=B=C=D.

# my dataset
mydata = data.frame(res=c(5,3,1,2,5,2,6,2),fac = c("A","A","B","B","C","C","D","D"))

# 1. using ANOVA
anova(lm(mydata$$res~mydata$$fac))

# output
Analysis of Variance Table

Response: mydata$$res Df Sum Sq Mean Sq F value Pr(>F) mydata$$fac  3    8.5  2.8333  0.7556 0.5742
Residuals   4   15.0  3.7500

# 2. using Wald test
library("aod")
lr = summary(lm(mydata$$res~mydata$$fac))
wald.test(Sigma = (lr$$cov.unscaled), b= (lr$$coefficients)[,1], Terms = 2:4)
# output
Wald test:
----------

Chi-squared test:
X2 = 8.5, df = 3, P(> X2) = 0.037


As you can see, although both test the same null hypothesis, the p-value from ANOVA is 0.5742 while the p-value from the wald test is 0.037. I am not sure if my calculation process is correct.

You're using the wald.test() function incorrectly (you should also mention what package the wald.test function is in since it's not in base R). You need to supply it with the variance-covariance matrix of the coefficients and the residual degrees of freedom. These are both in the lm object but not the summary.lm object. Here's how you can get the same results using both methods:
fit <- lm(mydata$$res ~ mydata$$fac)
wald.test(Sigma = vcov(fit), b = coef(fit), Terms = 2:4, df = fit$df.residual) ## Wald test: ## ---------- ## ## Chi-squared test: ## X2 = 2.3, df = 3, P(> X2) = 0.52 ## ## F test: ## W = 0.76, df1 = 3, df2 = 4, P(> W) = 0.57  You need to supply an argument to df to get the F-test results. Doing otherwise only produces the Chi-square results, which are only valid in large samples. • Thank you for pointing out my mistake. Your answer also helps me figure out another question posted by me. stats.stackexchange.com/questions/461365/… – purod Apr 24 '20 at 16:21 • An alternative one is to calculate the Wald statistics with the formula "t(coef(fit)[2:4]) %*% solve(vcov(fit)[2:4,2:4], coef(fit)[2:4])" and then apply pchisq(2.3, df=3, lower.tail=T) in case df$residual is difficult to obtain in some complicated model – purod Apr 24 '20 at 17:41