How to interpret squared term in GLMER model

Question

I'm trying to explain Y (a count variable), in terms of X, Xsquared, and Size, with random effects for Industry/Firm, by year. My hypothesis is that there is a positive relationship between Y and X, and that this relationship has a inverted-U shape (thus the Xsquared term). I normally do this with normal regression models, but I want to know if it would mean the same in a GLMER one (I'm new with these models)

Do my results, here below, support my inverted-U hypothesis, or something else? Any tips on how to show this inverted-U effect graphically in r? Lastly, how can I interpret the random effects part?

My GLMER model looks like this:

model <- glmer(Y ~ Year + X + Xsquared + Size + (1 + Year|Industry/Firm), data = mydata, family = poisson)

Results:

Random effects:
 Groups      Name        Variance    Std.Dev.  Corr 
 Firm:Industry (Intercept) 1.787e+04 1.337e+02      
             Year          4.434e-03 6.659e-02 -1.00
 Industry      (Intercept) 7.749e-01 8.803e-01      
             Year          1.923e-07 4.385e-04 -1.00
Number of obs: 436, groups:  Firm:Industry, 109; Industry, 37

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 51.697639   9.178129   5.633 1.77e-08 ***
Year        -0.027132   0.004535  -5.983 2.19e-09 ***
X            1.322702   0.334092   3.959 7.52e-05 ***
Xsquared    -0.277335   0.141129  -1.965   0.0494 *  
Size         0.321026   0.043825   7.325 2.39e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr) Year   DOI    DOI2  
Year     -0.999                     
DOI       0.004 -0.019              
DOI2     -0.009  0.020 -0.953       
LNAssets -0.162  0.119  0.003 -0.007

Updated results after changing Year to cYear

Random effects:
 Groups        Name        Variance  Std.Dev. Corr
 Firm:Industry (Intercept) 5.840e-01 0.764205     
               cYear       5.546e-03 0.074474 0.38
 Industry      (Intercept) 8.243e-05 0.009079     
               cYear       2.011e-05 0.004485 0.54
Number of obs: 436, groups:  Firm:Industry, 109; Industry, 37

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.87753    0.47970  -5.999 1.99e-09 ***
cYear       -0.02303    0.02778  -0.829    0.407    
X            1.32358    0.33632   3.935 8.30e-05 ***
Xsquared    -0.27628    0.14217  -1.943    0.052 .  
FirmSize     0.31789    0.04989   6.371 1.87e-10 ***

Answer 1

Going through your results a bit at a time:

Random effects:
 Groups      Name        Variance    Std.Dev.  Corr 
 Firm:Industry (Intercept) 1.787e+04 1.337e+02      
             Year          4.434e-03 6.659e-02 -1.00
 Industry      (Intercept) 7.749e-01 8.803e-01      
             Year          1.923e-07 4.385e-04 -1.00

The very small standard deviations of Year, and the perfect negative correlation between year and intercept within each group, suggest that your model is overfitted. I would recommend dropping the Year term and reverting to an intercept-only model (although it probably won't hurt anything if you insist on leaving it there). The among-firm-within-industry standard deviation is very large (approx. 130); if this is really a Poisson regression with a log link, it suggests that your response variable consists of very large and variable numbers of counts.

I would definitely check for outliers in your data, and strongly consider checking for overdispersion/fitting a negative-binomial model/if the counts are large enough, consider just a linear model on log-transformed data.

The variation among Industries is much smaller than the variability at the level of Firm-within-Industry -- not negligible (it's larger than some of your effect sizes).

Number of obs: 436, groups:  Firm:Industry, 109; Industry, 37

This is a reasonable-sized (I would call it 'small to medium') data set for mixed modeling.

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 51.697639   9.178129   5.633 1.77e-08 ***
Year        -0.027132   0.004535  -5.983 2.19e-09 ***
X            1.322702   0.334092   3.959 7.52e-05 ***
Xsquared    -0.277335   0.141129  -1.965   0.0494 *  
Size         0.321026   0.043825   7.325 2.39e-13 ***

Everything is highly significant (except for your X-squared term, which is below the magic $p<0.05$ but just barely!). The intercept is huge, but this is an artifact of the fact that you have Year in your model (probably) as years-since-AD 0; see below. When X=0, the effect of X on the log-counts is positive (1.32 log-counts per unit increase in X), but the curve decelerates/effect of X decreases with increasing X (the effect of X is zero, or the predicted curve of log-counts vs. X peaks, at $X=1.32/0.28 \approx 5$).

Correlation of Fixed Effects:
         (Intr) Year   DOI    DOI2  
Year     -0.999                     
DOI       0.004 -0.019              
DOI2     -0.009  0.020 -0.953       
LNAssets -0.162  0.119  0.003 -0.007

You're probably putting year in as a Gregorian-calendar year (e.g. 2015); this puts the intercept of the model at AD 0, which results in a very strong correlation between the intercept and the Year slope. This probably isn't messing anything up, but it would be clearer and would reduce the chances of numerical problems to set cYear=Year-mean(Year) (or subtract the minimum value, or the median, or any other reasonable value that means the zero value of Year is somewhere near your data).

For graphical solutions, try using predict() to plot predicted values from the model vs. X.

Answer 2

The interpretation of a squared term is a difference in differences comparing groups differing in 1 unit of the outcome. For instance if the model is linear, and we fit $E[Y|X] = \beta_0 + \beta_1 X$ then $\beta_1 = E[Y|X=9] - E[Y|X=8]$ as well as $\beta_1 = E[Y|X=1] - E[Y|X=0]$. But if the true model is quadratic, then that difference is not constant among all values.

Case 1: model misspecification, if the true model is $E[Y|X] = \beta_0 + \beta_1 X + \beta_2 X^2$ then $E[Y|X=x+1] - E[Y|X=x] = \beta_1 + \beta_2 (2X+1)$. Kind of complicated, but you see the "linear slope" is now a linear function of the $X$. But if you do a difference in differences for $(E[Y|X=x+2] - E[Y|X=x+1]) - (E[Y|X=x+1] - E[Y|X=x]) = 2\beta_2$. So basically the $\beta_1$ is the tangent slope of the quadratic curve at the origin, and $\beta_2$ is a quadratic slope.

With Poisson GLMs, the coefficients are exponentiated and interpreted as relative rates. So $e^{\beta_2}$ would be called a ratio of ratio of rates comparing groups differing by 1 unit differing by 1 unit of $X$. It's basically an interaction term.

The only distinction in mixed models is the interpretation of individual level effects.