Search code examples
rpredictglmm

GLM: Finding x value (predictor) for a particular y value (response)


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

Solution

  • 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.

    simplified example (construct data and fit model)

    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)
    

    use emmeans to get intercepts and slopes

    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()
    

    finish

    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 ...