# Standard errors in beta regression

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

• For those who are interested in coding /writing their own log-likelihood function for a Beta regression, I found that standard a optimisation routine (i.e., BFGS) could "easily" run into estimation issues because of NaN/Inf values generated by extreme values for the parameters of the Beta distribution (i.e., mean (mu) and precision (pr)). One way to overcome this limitation is to manually add some constraints (lower & upper limits) on possible values for "mu" and "pr". Here is what I did and it worked fine (so far): – Umka May 17 '17 at 22:00
• trick = 10^-15 pr = max(min(exp(x2),100), trick) mu = max(min(1 ./ (1 + exp(-x1)), 1-trick), trick) – Umka May 17 '17 at 22:00

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)

• The issue comes from the Hessian matrix When estimating with R, we obtain expected a block diagonal matrix for the "joint" model (e.g. elements (3:4,1:2) are null) But not the case with Matlab - The "off-block" elements are not null (sorry I cann't add a picture in a comment) For info, I am using "fmincon" command with Matlab – Umka May 18 '17 at 16:53

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).