Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Small bug in MATLAB R2017B LogLikelihood after fitnlm?

Background: I am working on a problem similar to the nonlinear logistic regression described in the link [1] (my problem is more complicated, but link [1] is enough for the next sections of this post). Comparing my results with those obtained in parallel with a R package, I got similar results for the coefficients, but (very approximately) an opposite logLikelihood.

Hypothesis: The logLikelihood given by fitnlm in Matlab is in fact the negative LogLikelihood. (Note that this impairs consequently the BIC and AIC computation by Matlab)

Reasonning: in [1], the same problem is solved through two different approaches. ML-approach/ By defining the negative LogLikelihood and making an optimization with fminsearch. GLS-approach/ By using fitnlm.

The negative LogLikelihood after the ML-approach is:380

The negative LogLikelihood after the GLS-approach is:-406

I imagine the second one should be at least multiplied by (-1)?

Questions: Did I miss something? Is the (-1) coefficient enough, or would this simple correction not be enough?

Self-contained code:

%copy-pasting code from [1]
        myf = @(beta,x) beta(1)*x./(beta(2) + x);
        mymodelfun = @(beta,x) 1./(1 + exp(-myf(beta,x)));
        rng(300,'twister');
        x    = linspace(-1,1,200)';
        beta = [10;2];
        beta0=[3;3];
        mu = mymodelfun(beta,x);
        n = 50;
        z = binornd(n,mu);
        y = z./n;
%ML Approach
        mynegloglik = @(beta) -sum(log(binopdf(z,n,mymodelfun(beta,x))));
        opts = optimset('fminsearch');
        opts.MaxFunEvals = Inf;
        opts.MaxIter = 10000;
        betaHatML = fminsearch(mynegloglik,beta0,opts)
        neglogLH_MLApproach = mynegloglik(betaHatML);

%GLS Approach
        wfun = @(xx) n./(xx.*(1-xx));
        nlm = fitnlm(x,y,mymodelfun,beta0,'Weights',wfun)      
        neglogLH_GLSApproach = - nlm.LogLikelihood;

Source:

[1] https://uk.mathworks.com/help/stats/examples/nonlinear-logistic-regression.html

like image 667
Alexandre C-L Avatar asked Feb 20 '20 12:02

Alexandre C-L


1 Answers

This answer (now) only details which code is used. Please see Tom Lane's answer below for a substantive answer.

Basically, fitnlm.m is a call to NonLinearModel.fit.

When opening NonLinearModel.m, one gets in line 1209:

model.LogLikelihood = getlogLikelihood(model);

getlogLikelihood is itself described between lines 1234-1251.

For instance:

 function L = getlogLikelihood(model)
            (...)
                L = -(model.DFE + model.NumObservations*log(2*pi) + (...) )/2;
            (...)

Please also not that this notably impacts ModelCriterion.AIC and ModelCriterion.BIC, as they are computed using model.LogLikelihood ("thinking" it is the logLikelihood).

To get the corresponding formula for BIC/AIC/..., type:

edit classreg.regr.modelutils.modelcriterion
like image 105
Alexandre C-L Avatar answered Sep 29 '22 07:09

Alexandre C-L