Hat Matrix of Beta Regression

Question

I am trying to calculate beta regression by using Hat matrix, manually. I am using idea of https://stats.stackexchange.com/questions/8126/how-to-calculate-the-hat-matrix-for-logistic-regression-in-r?noredirect=1&lq=1 answer. But its not working for me. Further, I am not sure its correct or not. Kindly help.
My code

  require(betareg)
df<-data("ReadingSkills")
y<-ReadingSkills$accuracy
n<-length(y)
x1<-rnorm(n,0,1)
x2<-rnorm(n,0,1)

X<-cbind(1,x1, x2)


bfit <- solve(t(X) %*% X + u * diag(3)) %*% t(X) %*% y
bfit1<-betareg(accuracy ~ x1+x2, data = ReadingSkills, link="logit")
bfit
bfit1

v <- 1/(1+exp(-X%*%bfit))
VX <- X*v 
H <- VX%*%solve(crossprod(VX,VX),t(VX))

Answer 1

This problem is best suited for cross validated as it requires math equations for the proofs. But we can still do the code here. Anyway, since your interest is the coefficients, it is worthwile to not that just like logistic regression, betareg does not have a closed form solution for the MLEs. You need to use any iterative methods for optimization. Below I use the readily available optim function from base R.

y<-ReadingSkills$accuracy
set.seed(1)
n<-length(y)
x1<-rnorm(n,0,1)
x2<-rnorm(n,0,1)
X<-cbind(1,x1, x2)



bfit1<-betareg(accuracy ~ x1+x2, data = ReadingSkills, link="logit")

bfit1

Call:
betareg(formula = accuracy ~ x1 + x2, data = ReadingSkills, 
    link = "logit")

Coefficients (mean model with logit link):
(Intercept)           x1           x2  
    1.30730      0.28936     -0.08727  

Phi coefficients (precision model with identity link):
(phi)  
3.409  

Can we do the same directly without depending on the betareg function? Yes of course. we just write the log likelihood function and optimize that for the 4 parameters, ie the 3 betas, and the phi coefficient:

ll <- function(par){
  mu <- c(plogis(X %*% head(par, -1)))
  phi <- tail(par, 1)
  -sum(dbeta(y, mu*phi, (1-mu)*phi, log = TRUE))
}

optim(runif(ncol(X) + 1),ll)

$par
[1]  1.30733153  0.28936064 -0.08723266  3.40954180

$value
[1] -27.84788

$counts
function gradient 
     169       NA 

$convergence
[1] 0

$message
NULL
  

You can see that the results from optim match exactly what we have in beta regression.