Covariance of parameter estimates in nonlinear mixed effects fit

I'm using nlmefit, and am trying to test the hypotheses using linhyptest. For that I need the covariance of the parameter estimates. The documentation suggests that the "stats" output structure should have a field called "covb". My output structure does not contain this!

[beta,psi,stats,re_subj] = nlmefit( X , Y, [], 
        @(b,X)X*b', [0,0,0,0], 'REParam',[1] ); 

stats = 

           dfe: 4844
          logl: -1.1235e+04
           mse: 5.9078
          rmse: 2.4274
    errorparam: 2.4306
           aic: 2.2482e+04
           bic: 2.2488e+04
        sebeta: [0.4017 0.0919 0.0426 0.0426]

How can I get the required matrix? (I was expecting a 4x4, with 'sebeta.^2' as diagonals maybe?)