How do I calculate the effect size of a fixed effect in a negative binomial mixed effect model?

Question

I am examining a certain fixed effect (effect3 here) by fitting two negative binomial mixed models using the glmmTMB function: one with and one without the fixed effect. Next, I perform a likelihood ratio test to test the significance of effect3 (example data):

eff1 <- c(0.026, 0.003, -0.008, -0.057, -0.022)

eff2 <- c(-0.002, -0.013, 0.036, 0.005, 0.074)

eff3 <- c(0.027, 0.021, -0.015, 0.008, 0)

counts <- c(18317, 19899, 11048, 23920, 20656)

data.eff <- data.frame(effect1 = eff1, effect2 = eff2, effect3 = eff3, Val = counts)

modeling <- function(m.data, vars){
  # This is how the model should look like
  f <- reformulate(termlabels = vars, response = outcome)
  model <- glmmTMB(formula = f, data = m.data, family = nbinom2())
}

outcome <- "Val"
variables <- colnames(data.eff[,-which(colnames(data.eff) == "Val")])
variableswo <- colnames(data.eff[,-which(colnames(data.eff) == "Val")])[-3]

model1 <- modeling(m.data = data.eff, vars = variables)
model2 <- modeling(m.data = data.eff, vars = variableswo)

res <- lrtest(model1, model2)

Sadly, I cannot find a way to evaluate the effect size; or rather none of the numerous potential approaches seem to work. One of the ideas that failed was calculating Cohen's f-squared:

r22 <- r2(model1)
r21 <- r2(model2)
f2 <- (r22-r21)/(1-r22)

However, r2(model) always results in NA values when using the above code. I think the bug is in the function, because the commands within the function work when run one-by-one but calling the arguments seems to crash it. Maybe it has got to do with the environment?

Therefore, my questions are:

1. Can I use glmmTMB() for a fixed effects model so a model without random effects?
2. How do I calculate such an effect size? Is Cohen's f2 the right choice?
3. How do I implement that in R given the code above? How can I fix that function?

Thank you!