Standard errors in beta regression

Question

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

Answer 1

Well ... I have made some progresses ... The issue of joint vs. separate estimation of the Beta models is now fixed - But only when I switch to R (instead of Matlab) so now I need to figure out what is wrong with computation of SEs using Matlab

In case you want to see my R codes (to reproduce what I have been doing so far): This first function is used to estimate the Beta model separately I have used "tricks" to "stabilise" the optimisation routine (i.e., to avoid computational issues with NaN/Inf values) but Achim Zeileis gave me some useful tips to avoid this issue (Thanks again!) which I need to implement now.

betalik = function(b, Y){
  # "Empty" model (i.e., intercept only)
  x1 = b[1]
  x2 = b[2]
  # Trciks to avoid estimation issues
  trick1 = 10 ** -15
  trick2 = 100
  # model log-likelihood
  mu = max(min(1 / (1 + exp(-x1)), 1-trick1), trick1) # Logit link
  pr = max(min(exp(x2),trick2), trick1)
  a1 = gamma(pr) / (gamma(mu*pr) * gamma((1-mu)*pr))
  a2 = Y ** (mu*pr - 1)
  a3 = (1-Y) ** ((1-mu)*pr - 1)
  llik = -sum(log(a1 * a2 * a3))
  return(llik)}

I used this second function to "jointly" estimate the Beta models - I am not talking about bivariate modelling - The proportion variables remain independent in this "joint" model (and then results from separate and joint models should be exactly the same! see my previous post)

betalik_V2 = function(b, Y){
  # Trciks to avoid estimation issues
  trick1 = 10 ** -10
  trick2 = 50
  # "Joint" log-likelihood
  lik = matrix(NA,nrow=nrow(Y),ncol=ncol(Y))
  for(i in 1:ncol(Y)){
    x1 = b[2*i-1]
    x2 = b[2*i-0]
    mu = max(min(1 / (1 + exp(-x1)), 1-trick1), trick1)
    pr = max(min(exp(x2),trick2), trick1)
    a1 = gamma(pr) / (gamma(mu*pr) * gamma((1-mu)*pr))
    a2 = Y[,i] ** (mu*pr - 1)
    a3 = (1-Y[,i]) ** ((1-mu)*pr - 1)
    lik[,i] = a1 * a2 * a3}
  llik = -sum(log(prod(lik,2)))
  return(llik)}

Here is a function to compute SEs, Pval and reshape the results:

NORSE = function(Model){
  Hess = Model$hess
  MLE = as.numeric(Model$par)
  Para = names(Model$par)
  SE = as.numeric(sqrt(diag(solve(Hess))))
  Pval = 2*(1-pnorm(abs(MLE/SE)))
  Pval = ifelse(Pval < .001, '< 0.001', round(Pval, digit=3))
  cinf = MLE - 1.96*SE
  csup = MLE + 1.96*SE
  data.frame(Para, MLE, SE, Pval, cinf, csup)}

We are getting there! Now few lines of code to simulate two Beta variables

Y1 = rbeta(n=1000, shape1=1, shape2=0.5)
Y2 = rbeta(n=1000, shape1=0.5, shape2=1)
Y = cbind(Y1,Y2)
plot(hist(Y1))
plot(hist(Y2))

And to maximise their log-likelihood:

sv = matrix(0, ncol=2)
names(sv) = c('mu','pr')
res1 = optim(par=sv, Y=Y1, fn=betalik, method='BFGS', control=list(trace=T, REPORT=1), hessian=T)
res2 = optim(par=sv, Y=Y2, fn=betalik, method='BFGS', control=list(trace=T, REPORT=1), hessian=T)
res1 = NORSE(res1)
res2 = NORSE(res2)

I have also compared results with those obtained from "betareg" - They are the same, except for phi (or what I call "pr" in my codes) because we use different link functions (identity vs. log)

Finally the last piece of code to "jointly" estimate the 2 Beta models:

sv = matrix(0, ncol=4)
names(sv) = c('mu1','pr1','mu2','pr2')
res3 = optim(par=sv, Y=Y, fn=betalik_V2, method='BFGS', control=list(trace=T, REPORT=1), hessian=T)
res3 = NORSE(res3)

As expected results are similar for both MLE and SE - With Matlab it is not true for the SEs ... why? ... MY Matlab version of the NORSE function seems identical:

function res = NORSE(Model)
  Hess = Model.h;
  MLE = Model.b';
  SE = sqrt(diag(inv(Hess)));
  P = round(2*(1-normcdf(abs(MLE ./ SE))),3);
  SV = [Model.sv; Model.nbid; Model.obs; Model.para; round(Model.ll,1); round(Model.bic,1)];
  PARA = [Model.names'; 'indiv'; 'obs'; 'para'; 'llik'; 'bic'];
  MLE = [MLE; zeros(5,1)];
  SE = [SE; zeros(5,1)];
  P = [P; zeros(5,1)];
  res = table(PARA, SV, MLE, SE, P);
end

Hope this helps! (and thanks again to Achim Zeileis for his help)

Answer 2

End of the story!
"fmincon" is a command for constrained optimisation, and it only provides a BFGS approximation to the Hessian.
The resulting Hessian matrix should not be used for any comparisons or computations outside of the optimisation routine.
Say otherwise this approximate Hessian matrix should not be used to derive SEs (as I initially did).
Two options:
1/ To switch to a command for unconstrained optimisation (e.g., "fminunc")
2/ To first use "fmincon" to obtain the final solution and then to re-estimate the model with "fminunc" with estimates of the final solution as starting values
I think the interest of using "fmincon" in the first place is to decrease computation time (I am new to matlab so not 100% sure of this).