optimize()

[Matt Piercy]

Hi Stata users,
I am trying to get mata to use a formula for REML from a textbook using the K commutation matrix reml.do which is different from the process used in mixed. I have got it to work on a simple dataset but when I test it with optimise() I keep getting errors due to an argument mismatch. Do file attached.
Appreciate any help. GitHub couldn't solve it.
Thanks,
Matt

[Choonjoo Lee]

Hi Matt Piercy ,

The argument mismatch is in the evaluator interface. With

 optimize_init_argument(S, 1, D)

a d0 evaluator must have the arguments

 (todo, b, D, value, gradient, Hessian)

in that order. optimize() does not pass an optimizer-state argument M at the beginning. For d0, the gradient and Hessian arguments are still included, although your function only needs to fill in the objective value. See help mf_optimize, particularly the sections on type d evaluators and extra arguments.

[Choonjoo Lee]

Hi Matt Piercy
Here is a complete example using your data and likelihood expression. It replaces the data in memory, so save your current data first. I have omitted the grid search and Kmat() demonstration, which are not used in the optimization below.

version 16.1
clear all
set more off

input double y x id
10 1 1
12 2 1
11 3 1
20 1 2
22 2 2
21 3 2
30 1 3
31 2 3
33 3 3
end

mata:
real scalar reml_ll(real scalar su2, real scalar se2,
                    real colvector y, real matrix X, real matrix Z)
{
    real matrix V, M, P
    real scalar n, p

    n = rows(X)
    p = cols(X)
    V = su2*Z*Z' + se2*I(n)
    M = invsym(V)
    P = M - M*X*invsym(X'*M*X)*X'*M

    return(-0.5*((n-p)*ln(2*pi()) + ln(det(V)) +
        ln(det(X'*M*X)) + y'*P*y))
}

void reml_obj_opt(real scalar todo, real rowvector b,
                  real matrix D, real scalar fv,
                  real rowvector g, real matrix H)
{
    real colvector yy
    real matrix XX, ZZ

    yy = D[,1]
    XX = D[,2..3]
    ZZ = D[,4..cols(D)]
    fv = -reml_ll(exp(b[1]), exp(b[2]), yy, XX, ZZ)
}

y = st_data(., "y")
X = (J(rows(y),1,1), st_data(., "x"))
id = st_data(., "id")

levels = uniqrows(sort(id,1))
Z = J(rows(id), rows(levels), 0)
for (j=1; j<=rows(levels); j++) Z[,j] = (id:==levels[j])

D = (y, X, Z)
S = optimize_init()
optimize_init_evaluator(S, &reml_obj_opt())
optimize_init_evaluatortype(S, "d0")
optimize_init_argument(S, 1, D)
optimize_init_params(S, (ln(10), ln(1)))
optimize_init_which(S, "min")

bhat = optimize(S)
assert(optimize_result_converged(S))
su2 = exp(bhat[1])
se2 = exp(bhat[2])
llopt = reml_ll(su2, se2, y, X, Z)

printf("su2 = %12.6f\n", su2)
printf("se2 = %12.6f\n", se2)
printf("REML log likelihood = %14.9f\n", llopt)
end

mixed y x || id:, reml variance

In Stata/MP 16.1, the Mata calculation gives approximately:

su2 = 103.070253
se2 =   0.900000
REML log likelihood = -17.400591123

The variance estimates agree with mixed, reml to numerical tolerance, and the restricted log likelihood agrees. Thus, for this random-intercept example, your likelihood expression and mixed lead to the same REML solution.

Two other points may help. Your grid stops at su2 <= 100, whereas the optimum is approximately 103.07, so that grid cannot reach it. Also, use mata: for these blocks and check convergence before reporting results. In the original file, execution continued after the optimize() error, and the final printed su2 and se2 were values left over from the grid loop, not optimization estimates.

Kmat() is not called by reml_ll() or the evaluator in the attachment. If your aim is to reproduce a particular textbook formula involving the commutation matrix, the book reference and equation would help clarify that part.

The example above retains your two-column X; the XX and ZZ column ranges would need changing if you add fixed-effect regressors. For larger problems, Cholesky-based solves and log determinants would also be preferable to explicitly forming inverses and determinants.

[Matt Piercy]

Hi Choojoo,
That's great it worked. Thanks so much, Matt (btw K is not commutation but cancelling matrix in this case - my mistake)

[Matt Piercy]

P = K(K'VK)^-1K' in the book (p177 Generalized, Linear, and Mixed Models (2008) McCullock, Searle and Neuhaus. You have used P = M - MX(X'MX)X'M with M=V^-1. These are the same (p350 same reference).

[Matt Piercy]

I got the code to work for sparse K using the formula for P above. reml2.do

[Choonjoo Lee]

Hi Matt Piercy ,

Thanks for posting the revised example. I reproduced its agreement with mixed, reml. One additional condition on K is worth checking: besides K'X=0, K must have rank n-p.

In reml2.do, K=(C*M)' satisfies the orthogonality condition, but the projection can reduce its rank. For example, swapping just the last two input rows, 31 2 3 and 33 3 3, leaves the dataset unchanged apart from its order, yet rank(K) falls from 7 to 6. The initial likelihood then becomes missing and optimize() reports "initial values not feasible".

Here is a small check using the same x values. It leaves the Stata dataset unchanged:

version 16.1
capture mata: mata drop reml_kcheck()
mata:
void reml_kcheck()
{
    real colvector x
    real matrix X, M, C, K, Q, R, Kqr
    real scalar n, p, i, trial

    x=(1\2\3\1\2\3\1\2\3)
    n=rows(x)
    p=2
    for(trial=1; trial<=2; trial++) {
        if(trial==2) x=x[(1\2\3\4\5\6\7\9\8),.]
        X=(J(n,1,1),x)
        M=I(n)-X*invsym(X'*X)*X'
        C=J(n-p,n,0)
        for(i=1; i<=n-p; i++) {
            C[i,i]=1
            C[i,i+1]=-1
        }
        K=(C*M)'
        Q=R=.
        qrd(X,Q,R)
        Kqr=Q[,p+1..n]
        printf("order %g: rank(X)=%g, rank(K)=%g, rank(Kqr)=%g\n",
            trial,rank(X),rank(K),rank(Kqr))
    }
}
reml_kcheck()
end

In Stata/MP 16.1 this gives:

order 1: rank(X)=2, rank(K)=7, rank(Kqr)=7
order 2: rank(X)=2, rank(K)=6, rank(Kqr)=7

For this small example, one option is to construct K from the full QR decomposition. Add Q and R to the real matrix declaration inside reml_ll(), and replace the C/M/K construction with:

assert(n>p & rank(X)==p)
Q=R=.
qrd(X,Q,R)
K=Q[,p+1..n]

Apply the same change to the separate K-display block. The likelihood expression and evaluator can otherwise remain as they are. With this change, both observation orders give su2 approximately 103.0704, se2 approximately 0.9 and restricted log likelihood approximately -17.400591123, agreeing with mixed, reml to numerical tolerance.

I would also start the block with mata: rather than mata, so that a Mata error stops the do-file. In the reordered example, the current interactive block continues to the later mixed command after the optimization error.

The QR option constructs a dense basis; a sparse construction would also need to retain all n-p independent contrasts. See help mf_qrd for the full QR decomposition.