Compute log likelihood of data after fitting the GLM? Please forgive my inexpertise in stats. I am playing around with a toy general linear model (GLM) with MATLAB glmfit.
X=[1;2;3;4];
y=[1.1;2.2;3.5;4.9];
[b, dev, stats] = glmfit(X, y);

How should I compute the log likelihood of the data? I want to do
y_fit = glmval(b, X, 'identity');
ll = nansum(log(normpdf(y, y_fit, sigma)));

But how to find sigma, the standard deviation of the Gaussian?
 A: In order to compute the maximized likelihood, one needs the maximum likelihood estimate (MLE) of the error standard deviation. It can be computed from the GLM residuals, stats.resid:
sigma = sqrt(mean(stats.resid .^ 2));

Using nansum should not normally be necessary, sum should do. To avoid problems with arithmetic underflow, it would be better though to work with the log-density of the normal distribution:
ll = sum(-0.5 * (stats.resid / sigma) .^ 2 - log(sqrt(2*pi) * sigma));


In a previous version of this answer I wrote:
sigma = sqrt(sum(stats.resid .^ 2) / stats.dfe);

This is not correct, since it gives the square root of the unbiased estimator of the residual variance, not the MLE of sigma.

In response to the poster's comment: The standard deviation can be computed this way because the MLE of the standard deviation of a normal distribution is the root-mean-square of deviations. A nice derivation for the variance of a simple normal distribution model (a special case of the GLM) can be found at statlect. Maximum likelihood estimation is invariant under transformations, and therefore the MLE of the standard deviation is just the square root of the MLE of the variance. For the full GLM the derivation is slightly more complicated, but with the same result: Instead of using the deviation from the mean one uses the deviations from the linear fit, and these are stored by glmfit in stats.resid.
You can check the correctness of the computation numerically in Matlab like this. Implement a function that computes the (negative) log likelihood for specific parameters b and sigma, and then use an optimization function to minimize it (find the MLE estimate of the parameters):
function mleGLM

X=[1;2;3;4];
y=[1.1;2.2;3.5;4.9];

[param, fval, exitflag] = fminunc(@nll, [0 0 1], optimset('TolX', 1e-10));

    function nll = nll(param)
        b = param(1 : end - 1);
        sigma = param(end);
        resid = X * b(2 : end) + b(1) - y;
        nll = - sum(-0.5 * (resid / sigma) .^ 2 - log(sqrt(2*pi) * sigma));
    end
end

You will find that param(3), which corresponds to sigma, is identical (up to numerical error) to the result of the formula above; and param(1 : 2) correspond to the b output of glmfit.
