I am trying to estimate Beta regression model with Matlab based on Ferrari & Cribari-Neto (2004) paper (see https://www.jstatsoft.org/article/view/v034i02/v34i02.pdf).
I have encountered a severe issue with the computation of the standard errors (SE) for maximum likelihood estimates (MLE).
In my data set I have 7 proportion (0-1) variables for 1,008 individuals. I first estimated an "empty" beta model (i.e., constant only for the mean and dispersion parameters of the Beta distribution) for each variable separately. When I compared my results with those obtained from "betareg" package (R software), there are very similar (The SE are slightly different but the optimisation strategies are also slightly different).
Mathematically speaking, "estimating separately these 7 models" and "specifying a model in which the likelihood of the 7 items are jointly estimated" should lead exactly to the same results (In the "joint" model, the assumption of independence of the observations is preserved). I tried to verify this hypothesis by comparing the results between these two modelling strategies. The good news is that, as expected, I have obtained exactly the same results in terms of MLE and log-likelihood (The sum of the 7 separate log-likelihood values = log-likelihood of the "joint" model), but very different results for the SE!
I have tried to figure out where this could come from but I am running out of ideas ... Your help would be much appreciated!
MATLAB code for the "joint" model
function loglikelihood = Betareg(b, N, I, data)
% (b) is the vector of parameters to estimate
% (N) is the number of individuals (N=1008)
% (I) is the number of variables (I=7)
% data ~ Beta(mu, pr), where "mu" stands for the mean and "pr" for dispersion ("precision")
likelihood = zeros(N, I);
for i = 1:I % Loop specifying the likelihood for each item
x1 = b(2*i-1); % Specification of MU ==> constant only (i.e., "empty" model)
x2 = b(2*i-0); % Specification of PR ==> constant only (i.e., "empty" model)
mu = 1 ./ (1 + exp(-x1)); % Logit link function for the mean
pr = exp(x2); % Exp(.) transformation to prevent the precision parameter to be <= 0
a1 = gamma(pr) ./ (gamma(mu .* pr) .* gamma((1 - mu) .* pr)); % 1st part of the Beta PDF
a2 = bsxfun(@power, data(:,i)', mu .* pr - 1); % 2nd part of the Beta PDF
a3 = bsxfun(@power, 1-data(:,i)', (1-mu) .* pr - 1); % 3rd part of the Beta PDF
likelihood(:,i) = a1 .* a2 .* a3;
end
loglikelihood = -sum(log(prod(likelihood,2))); % In total the llik function includes 14 parameters (2 for each item)
end