I have a follow-up to the post Regression (logistic) in R: Finding x value (predictor) for a particular y value (outcome).
I have a more complex model:
b1 <- glmmTMB(lambda ~ species * period * ffd_stan + (1|SITE_COUNTY),
data = fives)
where species
(5 species) and period
(3 periods) are categorical and ffd_stan
is numeric. I want to calculate the predicted values of the predictors when lambda = 0 (essentially at what value of ffd_stan
does the population growth not change for each of the species and period combinations).
I know I can build a data frame to predict lambda for many values of ffd_stan
but I would love to be able to do this in a more efficient manner. I'm not sure either how to adapt the following code or whether there is a better solution for my data
findInt <- function(model, value) {
function(x) {
predict(model, data.frame(mpg=x), type="response") - value
}
}
uniroot(findInt(model, .5), range(mtcars$mpg))$root
You could figure out which parameters to combine linearly to get the intercepts and slopes for each species/period combination, but below I use emmeans
to simplify these calculations.
fives <- expand.grid(species = factor(1:2),
period = factor(1:2),
ffd_stan = c(-1,0,1))
set.seed(101)
fives$lambda <- rnorm(nrow(fives))
library(glmmTMB)
b1 <- glmmTMB(lambda ~ species * period * ffd_stan,
data = fives)
library(emmeans)
## intercepts
m1 <- emmeans(b1, ~species*period, at = list(ffd_stan = 0)) |> as.data.frame()
stopifnot(all.equal(m1[1,"emmean"], unname(fixef(b1)$cond[1])))
## slopes
m2 <- emtrends(b1, ~species*period, var = "ffd_stan") |> as.data.frame()
Now we want to solve a + b*x = c
for each group, which gives us x = (c-a)/b
:
cc <- 0.5
(cc-m1[,"emmean"])/m2[,"ffd_stan.trend"]
You should check these solutions against your less efficient answer ...