Is it possible to calculate x-intercept from a mixed model?

Question

I understand that the x-intercept can be calculated using $y = mx + b$ for a linear model. I am unsure if this is statistically appropriate for a mixed model with count data, given that counts cannot be negative and there are random effects to consider. I have seen examples of x-intercept calculations for count data with simple linear regressions, but I'm unsure if this method can be extended to mixed models.

Here is my model:

mod_6 <-
  glmmTMB(total_count ~  mean_temp + (1|month) + (1|spread_event), 
    family = nbinom1, data = dat_nc_ncb)
summary(mod_6)

Here is the output.

 Family: nbinom1  ( log )
Formula:          total_count ~ mean_ws + (1 | month) + (1 | spread_event)
Data: dat_nc_ncb
 AIC      BIC   logLik deviance df.resid 

1399.1   1415.6   -694.5   1389.1      194
Random effects:
Conditional model:
 Groups       Name        Variance Std.Dev.
 month        (Intercept) 0.3671   0.6059

 spread_event (Intercept) 0.3279   0.5726

Number of obs: 199, groups:  month, 10; spread_event, 26
Dispersion parameter for nbinom1 family ():  177
Conditional model:
            Estimate Std. Error z value Pr(>|z|)

(Intercept)   3.4928     0.3515   9.936   <2e-16 ***
mean_ws      -1.1099     0.5126  -2.165   0.0304 *

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Is it statistically accurate if extract the fixed effects coefficients using coefficients <- fixef(mod_6), identify the coefficient for the intercept using intercept <- coefficients[1], extract the slope using slope <- coefficients[2] and finally extract x-intercept using x_intercept <- -intercept/slope?

Or would be it more appropriate to use a simple glm with quassipoisson family, and then calculate x-intercept. That way, I won't have to worry about random effects?

Details about the experiment

I left out my potted plants in the field for a week, took them back to the glasshouse and counted the number of infected leaves per plant after two weeks. Plants are infected in ideal condition of temperature.

Analysis goal

I need to find lower temperature thresholds. More details can be found in figures 1-4 [here], (http://uspest.org/wea/Boxwood_blight_risk_model_summaryV21.pdf), but the basic idea is that we want to find out temperature at which no disease was observed (lower temperature threshold for disease). Since the goal is to find thresholds, I am happy to let go of the random effects if this allows me to calculate x-intercept for the mean_temp.

Answer 1

Unfortunately, this is difficult, for conceptual (rather than computational) reasons. Count models (negative binomial, Poisson, etc.), almost always use a logarithmic link function, meaning that they model the expected number of counts as an exponential function of some linear combination of covariates, in your case (because you only have one covariate: I'm ignoring the random effects here for simplicity, but they won't affect the argument)

$$ \textrm{expected counts} = \exp(\beta_0 + \beta_1 x_1) $$

The models are set up this way primarily because otherwise it would be possible for the expected number of counts to be negative, which would cause both computational and conceptual problems. However, this setup means that the expected-counts curve never intersects the $x$ axis — it just gets closer and closer to zero as the linear predictor ($\beta_0 + \sum \beta_i x_i$) becomes more and more negative.

One possibility would be to pick a small expected number of counts $C$ (say, 0.1 or 0.01) and define that as the effective minimum value; then you'd solve the equation $\log(C) = \beta_0 + \beta_1 x$ for $x$.

For example, if you have a single numeric predictor, then you would use coef() or fixef() (depending on the model type you're using) to extract the fixed-effect coefficients and then compute x_intercept = (log(C) - b[1])/b[2]. (Once you've done that, check a plot of the predicted values to confirm that the answer makes sense.)

Another possibility, which I don't recommend, would be to fit an identity-link model (family = nbinom1(link = "identity")); this would make expected counts equal to $\beta_0 + \beta_1 x$ rather than $\exp(\beta_0 + \beta_1 x)$.

Answer 2

While specifying a threshold in the expected value of the count response is one approach to identifying the required predictor values, I would focus on the discrete-event aspect of the response: Each nonnegative integer value correspondents to a predicted probability mass according to the fitted negative binomial model. For example,

Pr(y = 0 | x = 32) = .85
Pr(y = 1 | x = 32) = .10
Pr(y = 2 | x = 32) = .02
...

I recommend translating your research question into one that aims to find the range of x that makes Pr(y = 0 | x = x) > theta. Here theta is an exogeneous, fixed parameter of desired noninfection probability typically obtained by expertise. Then, you can interpret that "to control the chance of infection at or below 1 - theta, it requires a temperature no higher than x degrees; further climate changes that raise the temperature from A to B degrees will increase the risk of infection from 10% to 30%." A larger Pr(y = 0 | x) correspondents to a smaller E(y | x) = sum (k * Pr(y = k | x)) = exp(b * x). Therefore, while Ben Bolker's approach is reasonable, mine would be clearer and easier to comprehend since achieving no infection appears particularly interesting to you.

In addition to making model prediction wisely, your model could benefit from several other aspects:

0. Plot the data. Like any data analysis process, conduct descriptive analysis before modeling. To help you attack your research question, provide a plot of y over x, color the points by month, alter the shape by field, and adjust the size by the number of cases (if you have multiple plants sharing the same infections under the same temperature, so several y-x coordinates overlap). This gives you an idea what mean structure specifications might fit the data well. Examine carefully if you have enough cases of y = 0 to fit and estimate its probability accurately.

1. Check linearity and interactions in fixed effects. Like other discrete choice modeling such as logistic regression, the error term distribution in count models is assumed rather than estimated. Unlike in linear regression, you won't get consistent coefficient estimates unless your mean structure is correctly specified. Therefore, it is very important to check if the model formula is set correctly. Temperature in many subjects show nonlinear effects. In your case, set a lowest possible value or a reference temperature t0, so the temperature predictor is a difference from this reference level t - t0, consider adding at least a square, cubic, or logarithm term of temperature, and interact temperature with other numeric and categorical predictors like month. You may not want to use semiparametric approaches such as spline for nonlinear effects, as you do care about the predictor's functional form and need it to solve for desired temperatures.

2. Test other random effect specifications. Random effects are usually used in repeated or clustered measurements, but I cannot tell from your variable names whether any scenario should apply. Was the same pot of plant exposed in a field for multiple times? Was the same set of fields used for multiple plants? I would use months as a categorical or numerical predictor in the regular, fixed effect component instead of a random component, since there are 12 and only 12 months in a year that not a random sample from a larger pool despite only 10 appearing in your sample. It is also quite meaningful to see how the infection varies over a year even the temperature is held constant (e.g., 80 degrees in both April and August may lead to very different infection problems.) I will add months first as dummy variables and then try to simplify the pattern with numeric values and functional forms. I would also consider representing the field effect by neither a random nor a fixed indicator but a set of variable describing its characteristics (e.g. soil water content, pH, and humus) to increase the research impacts and applicability, assuming you used many different fields if 26 denotes the number of fields.

3. Evaluate exposure and offset. Reexamine your experiment design and see whether each observation correspondents to the same amount of exposure (e.g., duration of field exposure that affects the maximum possible value of counts). You mentioned that you counted the number of infected leaves per plant. I would consider using the total number of leaves as the exposure, as a plant with five leaves cannot possibly get infected in six leaves, whereas a plant with six or more leaves can. So add a term + offset(log(leaves)) into the formula to capture this variation in top infection potential.

4. Compare different error-term and dispersion assumptions. You presented a model in family = nbinom1, where the extremely large value of 177 of the dispersion parameter means that overdispersion is huge, as Var(y) = E(y) * (1 + 177) = 178 * E(y). You should also test family = nbinom2 to see if this assumption fits the data better. Anyway, with this much overdispersion family = poisson should outperform neither, and you may wonder whether omitted variables (such as temperature squared) in the mean structure contributed to this overdispersion and whether the dispersion parameter is associated with any predictors. There is a dispersion formula in glmmTMB(disp = ~ ) that you should explore. You can use the same predictors in both measure and dispersion structures, such as glmmTMB(count ~ temperature, disp = ~ temperature). Give explanations to overdispersion. One plausible reason is that infection attracts additional infection-- each leaf's infection event is not independent. Another common overdispersion cause is zero inflation if your sample has two groups of plants, one is infection resistant (always infection free) and the other is not but you did not assay which carries such genes and which does not. You can test zero inflation through glmmTMB(count ~ temperature, zi = ~ ...) with an intercept and optionally some predictors of infection resistance.

5. Consider alternative model types. To assess how sensitive your conclusions are to your model assumptions, you need to present alternative model results. In addition to check different combinations of mean structure, dispersion structure, and zero inflation structure, you may also manipulate the response in several ways: (1) Model the number of infected areas (possibly multiple on one leaf) using a count model. (2) Predict the proportional of leaves infected using a beta regression; there is also zero inflated beta regression. For this, glmmTMB has the beta family. (3) Predict infection severity with ordinal regression, especially if your observed response has limited unique values (such as only 0, 1, 2, 3, 4 leaves infected). For this, you wan to use ordinal package and clmm function with scale and nominal effects.

6. Deliver model interpretation with visualization. As the transformed research question suggests, you should plot predicted E(y), Pr(y = 0), and E(y | y >= 1) over x with confidence intervals to capture uncertainty, to see how infection and noninfection change with temperature. You may overlay curves created in different months, fields, ect. Make sure that you understand what a confidence interval includes or excludes a reference value means and does not mean. Note that nonsignificance does not establish equivalence. In your case, you may want to refer to noninferiority analysis popular in medical studies, as you will test "H0: Pr(y = 0 | x) <= theta." Essentially, you will build 90% two-sided confidence interval if you test this one-sided hypothesis at a 5% significance level. Then, your plot shows what temperatures keep noninfection probability above a certain threshold at 95% confidence. As the chance of false alarm will be juxtaposed with the probability of noninfection, ensure a good choice of significance level and a smart interpretation. You will find the R package marginaleffects useful because it produce predictions with standard errors, confidence intervals, and noninferiority p values. My understanding is that random effects are not included in these predictions that are population means instead of individual or group means. marginaleffects should also give Pr(y = 0) with uncertainty measures, but it is yet to be verified. You can also try the package ggeffects that does account for random effects (larger confidence intervals). The last time when I tried with nlme random-effect models, the confidence intervals with random effects are calculated incorrectly in ggeffects. I have not tried ggeffects with random-effect count models.