# Classical Hypothesis tests: Wald, LR, LM

**URL:** <https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575>\
**Category:** Statistics\
**Tags:** hypothesis-tests\
**Created:** [May 16, 2020, 7:19am UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575 "2020-05-16T07:19:18Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 16, 2020, 7:19am UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/1 "2020-05-16T07:19:18Z")

</div>

I was able to code up the [LR & Wald tests](https://github.com/JuliaStats/HypothesisTests.jl/issues/200).  
I need help w/ the score test. Does anyone know how to obtain the loglik functions after `glm`?

```julia
using DataFrames, GLM, Distributions

n=7; k=1; X=randn(n,k); β=[2 ;-4]; σ_e = 10.1;
Y= [ones(n) X]*β + σ_e*randn(n);

d = DataFrame(X1= X[:,1], Y=Y);
#d = DataFrame(X1= X[:,1], X2= X[:,2], Y=Y);
#
f0 = @formula(Y ~ 1)
fA = @formula(Y ~ 1 + X1)
#
mH0 = lm(f0, d)
mHA = lm(fA, d)
#
mH0 = glm(f0, d, Normal())
mHA = glm(fA, d, Normal())

# Likelihood Ratio Test
function HT_LR(mH0, mHA)
    LR = 2.0*(loglikelihood(mHA) - loglikelihood(mH0)) |> abs
    df = dof_residual(mH0) - dof_residual(mHA) |> abs
    pv = 1 - cdf(Chisq(df), LR)
    return LR, pv, df
end
#
HT_LR(mH0, mHA) #better practice
HT_LR(mHA, mH0)
# Wald Test
"https://en.wikipedia.org/wiki/Wald_test"
"H0: R*θ = r"
"HA: R*θ ≂̸ r"
#
function HT_Wald(mHA, R, r, V)
    θ̂ = coef(mHA)
    A = (R*θ̂ - r)
    #V = vcov(mHA) #[2,2]
    n = size(mHA.mm.m,1)
    W = A' * inv(R*V*R') * A #V/n
    df = size(R,1)
    pv = 1 - cdf(Chisq(df), W)
    return W, pv, df
end
R = [0 1]
r = [0]
HT_Wald(mHA, R, r, vcov(mHA))

"LM aka Score test Matlab code"
G = gradp(@nll_lin,theta_r_full,datamat,1);
H_r = HessMp(@nll_lin,theta_r_full,datamat,1);
V_r = inv(H_r);
LM = G*V_r*G';
LM_p = 1-chi2cdf(LM,2);

```

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [May 16, 2020, 10:04am UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/2 "2020-05-16T10:04:07Z")

</div>

You might find the code in [Econometrics/TestStatistics.jl at main · mcreel/Econometrics · GitHub](https://github.com/mcreel/Econometrics/blob/master/src/LinearRegression/TestStatistics.jl) useful, it has the score test for the linear regression model with classical assumptions.

```julia
julia> TestStatistics();
┌───────┬──────────────┬──────────────┐
│ │ Value │ p-value │
├───────┼──────────────┼──────────────┤
│ qF │ 0.93966 │ 0.29772 │
│ Wald │ 1.04407 │ 0.33153 │
│ LR │ 1.02631 │ 0.32581 │
│ Score │ 1.00895 │ 0.32021 │
└───────┴──────────────┴──────────────┘

julia> 

```

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [May 16, 2020, 1:07pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/3 "2020-05-16T13:07:38Z")

</div>

I think you want something like `loglik_obs` [here](https://github.com/JuliaStats/GLM.jl/blob/1844866761deded1d59daf6efe38050dc4906cdc/src/glmfit.jl#L262)?

Also StatsModels has something called HypothesisCoding for a matrix. Not sure what it means but seems like it could be useful.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 16, 2020, 6:48pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/4 "2020-05-16T18:48:14Z")

</div>

@mcreel thanks for the link.  
My understanding is that Wald/LR/LM tests are asymptotically equivalent but:  
Finite sample statistics: Wald \> LR \> LM (as in your table)  
Finite sample p-values: Wald \< LR \< LM (opposite of your table)  
I think your table can be fixed by  
replacing: `pvalues = chisqccdf.(tests,q)`  
with: `pvalues = ccdf.(Chisq(q), tests)`

I think you’re goal is to write out the entire likelihood for educational purposes.  
The LRT can work w/ any Julia object w/ `loglikelihood(m)`

```julia
function HT_LR(mH0, mHA)
    LR = 2.0*(loglikelihood(mHA) - loglikelihood(mH0)) |> abs
    df = dof_residual(mH0) - dof_residual(mHA) |> abs
    pv = 1 - cdf(Chisq(df), LR) # ccdf(Chisq(df), LR)
    return LR, pv, df
end

```

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [May 16, 2020, 8:24pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/5 "2020-05-16T20:24:14Z")

</div>

Thanks! You’re right about the p-values. I have fixed that.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 16, 2020, 8:43pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/6 "2020-05-16T20:43:34Z")

</div>

@mcreel should your qF statistic be distributed FDist(q, n-k) under the H0?  
(not sure about this)

**Update** :  
Bruce Hansen’s [Econometrics](https://www.ssc.wisc.edu/~bhansen/econometrics/Econometrics.pdf) p 240 discusses this. Basically your code is correct.

W(\theta) \sim\_{\text{asymptotically}} \chi^2\_q  
The F version of the test: F = W/q  
It’s conventional to use F\_{q,n-k} p-values instead of \chi^2\_q.  
…  
While there is no formal justification to using the F\_{q,n-k} distribution for non-homoskedastic covariance matrices, the F\_{q,n-k} distribution provides continuity with the exact distribution theory under normality and is a bit more conservative than the \chi^2\_q distribution. (Furthermore, the difference is small when n−k is moderately large.)

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [May 16, 2020, 11:25pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/7 "2020-05-16T23:25:08Z")

</div>

A few comments.

Depending on the structure of the covariance matrix it may be an F-test or a Wald-test.

The hypotheses and restrictions need to take into account the contrasts specified in the model.

The residual degrees of freedom for the test under the null hypothesis should come from the model  
`dof_residual` and there are various options. For example, ridge models would use effective degrees of freedom in the computation and linear mixed models can choose any of several approximations (there isn’t a “correct way” just flavored approximations which are more or less suitable depending on various factors).

For models that are not maximum likelihood (e.g. restricted maximum likelihood) there are a number of caveats such as no true LR-test and you can use the “on the deviance scale” metrics but requires the same structure of fixed effects.

The more general approach is the restrictions based on linear combinations.

The main difference between LR test and Wald is that the Wald test allows for testing any number of linear combinations and hypothesis tests including incorporating the information from the variance covariance estimates all with a single fitted model. The likelihood ratio test is more powerful but requires re-fitting the model for each combination.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 17, 2020, 12:10am UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/8 "2020-05-17T00:10:04Z")

</div>

@pdeffebach thanks. I tried that yesterday but gave up.  
Just tried again. Close, but no cigar (yet):

```julia
using Zygote
function AZll(m, β)
    r = m.model.rr
    y = r.y
    ll = 0.0 #zero(eltype(mu))
    ϕ = deviance(m)/length(y)
    x = m.model.pp.X
    @inbounds for i in eachindex(y)
        ll += logpdf(Normal(x[i,:]'*β, sqrt(ϕ)), y[i])
    end
    ll
end
AZll(mHA, coef(mHA))
AZ(β) = AZll(mHA, β)

AZ(coef(mHA))
AZ([coef(mH0); 0.0])
Sc= Zygote.forward_jacobian(AZ, [coef(mH0); 0.0])[2]
FI= -Zygote.hessian(AZ, [coef(mH0); 0.0])
LM = (Sc' * inv(FI) * Sc)[1]
LM ∼ Χ²(q)

```

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [May 17, 2020, 5:33pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/9 "2020-05-17T17:33:23Z")

</div>

Can you also try `Distribtions.gradloglikpdf`?

---

<div class="post-metadata">

**Author:** ![alexfakos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alexfakos/32/209024_2.png) [@alexfakos](https://discourse.julialang.org/u/alexfakos)\
**Post date:** [May 4, 2024, 4:13pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/10 "2024-05-04T16:13:40Z")

</div>

Did you ever figure out how to get the loglikelihood from GLM as a function of parameters so that you can calculate the score vector and the hessian numerically?

Is there a way to have a method of `GLM.loglik_obs` that maps from the model parameter vector to a vector equal to the sample size?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 4, 2024, 6:34pm UTC](https://discourse.julialang.org/t/classical-hypothesis-tests-wald-lr-lm/39575/11 "2024-05-04T18:34:35Z")

</div>

I haven’t looked at this in years. Sorry