Score test after MLE with normal distribution

[Mina Sami]

Hi all,
I am trying to run score test after running a model with MLE norma distribution.
I composed a program:
"
capture program drop prog3
program define prog3
version 13
args lnf xb sigma
quietly replace `lnf'= -.5*ln(2*_pi) - ln(`sigma')- .5* ( ($ML_y1 - `xb')/`sigma')^2
end
"

Then I run my command:
ml model lf prog3 (xb: Y= X1 X2 X3) (sigma: )
ml maximize


Now How can I run a score test?

I saw from STATA that there was a command called testomit but it doesn't work.


What I did exactly:
ml model lf prog3 (xb: Y= X1 X2 X3) (sigma: ), vce(oim)

. predict double score, score

But after this step I failed to find score test like I found lrtest and wald tests (by est store restricted vs unrestricted model).

Thanks in advance for your help and suggestions

[Maarten Buis]

Score tests aren't implemented in Stata. However, you can implement them yourself.

// create our maximization program
clear all
program define prog3
version 13
args lnf xb sigma
quietly replace `lnf'= ln(normalden($ML_y1, `xb', `sigma'))
end

// check with some example data
sysuse auto
replace price = price / 1000

ml model lf prog3 (xb: price= mpg foreign) (sigma: )
ml maximize

reg price mpg foreign

//================================ test hypothesis that b(foreign) = 0
//--------------------------------------------------------------- Wald

qui ml model lf prog3 (xb: price= mpg foreign) (sigma: ), maximize

test foreign = 0

//--------------------------------------------------- likelihood ratio
qui ml model lf prog3 (xb: price= mpg foreign) (sigma: ), maximize
est store unc

qui ml model lf prog3 (xb: price= mpg) (sigma: ), maximize

est store c

lrtest unc c

//-------------------------------------------------------------- score
// store constrained parameters in matrix
qui ml model lf prog3 (xb: price= mpg) (sigma: ), maximize
matrix b0 = e(b)

// use those parameters in the unconstrained model
// but don't maximze (iter=0)
qui ml model lf prog3 (xb: price= mpg foreign) (sigma: ), ///
             maximize init(b0) iter(0)

// use the resulting gradient and var-cov matrix to compute the test statistics
matrix S = e(gradient)*e(V)*e(gradient)'
scalar s = el(S,1,1)
scalar p_s = chi2tail(1, s)

// display the result
di as txt "chi2(1)     = " as result %9.2g `=s'
di as txt "Prob > chi2 = " as result %9.4g `=p_s'

However, as you can see it is quite a bit of work, and more importantly many opportunities to include a bug. So I would not do that, and just stick to the Wald and Likelihood ratio test.

[Mina Sami]

Dear Maarten, thank you for your help. It works!