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