0

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))
J AK
  • 65
  • 5
  • yes and I want coefficients manually. – J AK Mar 19 '23 at 14:08
  • I dont think so my problem is as much complicated....check my referred link. – J AK Mar 19 '23 at 15:55
  • In the link they are interested in the Hat matrix and not the coefficients. You are interested in the coefficients. Two separate problems – Onyambu Mar 19 '23 at 20:03
  • nxt i want hat matrix. But I am stuck on starting point. – J AK Mar 20 '23 at 14:09
  • issue of Hat matrix is solved. But still want to understand how you make a link function? – J AK Mar 20 '23 at 14:16
  • Check the code below and you will see that I used the `logit` link function for beta and the `log` link function for `phi`. – Onyambu Mar 20 '23 at 17:17
  • If I want to change link for betas thn.....? – J AK Mar 21 '23 at 14:27
  • Use a different function. Eg if you want probit, use `pnorm` instead of the `plogis`. Seems to me that the code i gave is of no help to you. Note that you first need to learn the distributions, then learn R extensively. I do have free R for learning R you could attend. I can hell you with R. – Onyambu Mar 21 '23 at 15:27
  • thankq very much. But at which plate form??? – J AK Mar 22 '23 at 02:01
  • Seems you did not even look at the solution. One thing. You need to start learning R. thats all i can say for now. – Onyambu Mar 22 '23 at 03:19

1 Answers1

0

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.

Onyambu
  • 67,392
  • 3
  • 24
  • 53
  • thanks u hve solved my 1st problem bt i didn't understand the code. Its necessary coz nxt I hve to make this ridge by adding shrinkage parameter. Kindly plz define at which stage that will be added. – J AK Mar 20 '23 at 14:04
  • @JAK The code simply computes the beta regressionparameters. First I wrote a log likelihood taking into consideration the logit link then I maximized it to obtain the MLEs. Please note that if at any point the code seems difficult, then you need to learn the theory behind beta regression and even learn more about R. Notice that I only used BASE R meaning the code should be straight forward. – Onyambu Mar 20 '23 at 17:16
  • If I want to make this Ridge thn @ which point I 'll add shirinkge parameter? plus is my 3 liner code for H is correct or not? – J AK Mar 21 '23 at 14:07
  • @Ridge is a linear regression with L2 penalty. Seems to me that you need more of statistics background to understand the difference. Use `glmnet` package to run ridge regression. No you should not mix ridge and beta regressions. – Onyambu Mar 21 '23 at 15:36
  • I want to mix ridge and Beta, which is not available in glmnet. coz my main task is H matrix of beta with ridge. – J AK Mar 22 '23 at 02:02
  • @JAK is this a research thing you are doing? Well i have never heard of that. Probably someone else will help you – Onyambu Mar 22 '23 at 03:18