How can I calculate standardized weighted residuals 2 for a beta family glmm in R? Hello and thank you for your time!
I’m trying to calculate the "sweighted2" residuals for a beta family glmm with varying dispersion (aka precision parameter aka phi) fit using glmmTMB in R. Unfortunately only the response residuals are provided when using residuals(), so I’m attempting to do this by hand - sort of...
The standardized weighted residual 2 formula is equation 6 in the betareg vignette, and number 2 in the image here.

First, does anyone know if the equation used for the "sweighted2" residuals of a glm with varying dispersion will still work for a glmm with varying dispersion? Does the addition of random effects in the overall model require any change in the equation?
If it is directly applicable, I think I have everything figured out except for getting the weights to calculate the hat matrix (no hat matrix is provided in initial output of model). I’ve been using the data provided in section 4.2 Variable dispersion model of the betareg vignette to test this out and compare with the residuals I’ve calculated manually to those produced by the function residuals().
Second, does anyone know where I went wrong in the following failed attempt to calculate h? I suspect its either the derivative of the link function or incorrect coding for H.
gy2 <- betareg(yield ~ batch + temp | temp, data = GasolineYield)
h = hatvalues(gy2) ## for comparison to the ones I calculate
mu<-fitted(gy2)  ### mu
eta = predict(gy2, type ="link")
ms <-(digamma(mu*phi))-(digamma((1-mu)*phi))
yi<-gy2$y  ### y_i
ys<-log(yi/(1-yi))
vt<-trigamma((mu*phi))+trigamma((1-mu)*phi)
r<-(ys-ms)/sqrt(vt*(1-h))  #### sw2 residuals that should match ‘residuals(gy2, type = "sweighted2")’
X<-cbind(gy2$model$batch, gy2$model$temp)  ## covariate matrix
Xt<-t(X)     ## transpose of X
p<-diag(phi)
gd = 1/(mu*(1-mu))  #trying to get derivative of the logit link function
wt = phi*vt*(1/((gd)^2))
W <- diag(wt)

H = ((W %*% p)^(1/2)) %*% X %*% ((Xt %*% p %*% W %*% X)^-1) %*% Xt %*% ((p %*% W)^1/2)
diag(H)
h    ### my H doesn’t match their h

 A: As for your first question: I'm fairly certain that you have to account for the random effects in the beta GLMM and cannot simply use the equations for standard beta regression. At least the hat matrix should look different.
Regarding the second question: The R code below shows how to compute
the residuals almost from scratch.
## data and model
library("betareg")
data("GasolineYield", package = "betareg")
gy2 <- betareg(yield ~ batch + temp | temp, data = GasolineYield)

## basic quantities
y <- GasolineYield$yield
x <- model.matrix(gy2)
mu <- predict(gy2, type = "response")
phi <- predict(gy2, type = "precision")

## transformations
ystar <- qlogis(y)
mustar <- digamma(mu * phi) - digamma((1 - mu) * phi)
v <- trigamma(mu * phi) + trigamma((1 - mu) * phi)

## derivative of inverse link function (for logit link)
mu.eta <- dlogis(qlogis(mu))
## more generally you could use:
## link <- gy2$link$mean
## mu.eta <- link$mu.eta(link$linkfun(mu))

## hat values
w <- phi * v * mu.eta^2               ## weights W
xwx1 <- solve(crossprod(sqrt(w) * x)) ## (X'W X)^(-1)
h <- w * x %*% xwx1 %*% t(x)          ## hat matrix
h <- diag(h)                          ## diagonal

## put together residuals
res <- (ystar - mustar) / sqrt(v * (1 - h))

## check
all.equal(res, residuals(gy2))

