How to interpret the output of a Generalized Linear Model with R lmer

Question

I am having trouble understanding the output of a GLM I am trying to run with R package lme4. Here is an example of what I would like to achieve with some mock data:

Suppose that I have 3 geographical regions (north, south and east) in which I have planted some tomato plants, apple trees and orange trees. At each location, I have planted 5 of each. At harvest time, I have measured the number of fruits given by each plant and also the height of each plant. I have repeated this process during 4 years, and I am treating each year as a different replicate of the same experiment. My data thus looks like this:

plant region      year height n_fruits
1 tomato  north year_2001   14.5      138
2 tomato  north year_2003   11.5      107
3 tomato  north year_2005   12.0      116
4 tomato  north year_2007   13.0      110
5 tomato  north year_2001   15.0      132
6 tomato  north year_2003   11.0      114

Here is the full data frame for reproducibility:

data <- structure(list(plant = c("tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", "tomato", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "apple", "apple", "apple", 
                                 "apple", "apple", "apple", "apple", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange", "orange", "orange", "orange", "orange", "orange", "orange", 
                                 "orange"),
                       region = c("north", "north", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "south", "south", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "east", "east", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "north", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "north", "south", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "south", "east", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "east", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "north", "north", "north", "north", "north", "north", 
                                  "north", "north", "north", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "south", "south", "south", "south", "south", "south", "south", 
                                  "south", "south", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "east", "east", "east", "east", 
                                  "east", "east", "east", "east", "east", "east"),
                       year = c("year_2001",
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007", "year_2001", 
                                "year_2003", "year_2005", "year_2007", "year_2001", "year_2003", 
                                "year_2005", "year_2007", "year_2001", "year_2003", "year_2005", 
                                "year_2007", "year_2001", "year_2003", "year_2005", "year_2007", 
                                "year_2001", "year_2003", "year_2005", "year_2007"),
                       height = c(14.5,
                                  11.5, 12, 13, 15, 11, 14.5, 15.5, 13, 13, 10, 11, 11, 13.5, 12, 
                                  14, 12.5, 13.5, 16.5, 12, 14, 15, 11, 13, 10.5, 11.5, 12, 8.5, 
                                  12, 14.5, 12, 12.5, 13, 12.5, 11, 14, 13.5, 14, 10.5, 13.5, 12, 
                                  14, 13, 14, 12.5, 12.5, 14, 9, 12.5, 13.5, 13.5, 12.5, 14.5, 
                                  12, 11.5, 10, 10, 11.5, 12.5, 13, 24, 30, 23, 25, 23, 26, 23, 
                                  25, 28, 20, 29, 24, 29, 24, 23, 25, 29, 29, 24, 28, 31, 24, 27, 
                                  24, 23, 27, 22, 27, 21, 23, 21, 22, 19, 26, 29, 28, 28, 25, 24, 
                                  28, 26, 26, 24, 23, 34, 26, 22, 21, 25, 30, 26, 32, 27, 24, 24, 
                                  21, 17, 27, 21, 25, 52, 66, 50, 50, 44, 54, 50, 50, 44, 44, 52, 
                                  52, 40, 38, 52, 60, 52, 52, 50, 66, 50, 54, 52, 44, 46, 54, 50, 
                                  44, 54, 42, 58, 52, 52, 46, 50, 50, 44, 50, 40, 46, 44, 46, 58, 
                                  50, 56, 58, 48, 40, 46, 54, 48, 52, 58, 58, 48, 48, 58, 52, 54, 
                                  52),
                       n_fruits = c(138, 107, 116, 110, 132, 114, 120, 143, 121, 
                                    129, 116, 107, 108, 121, 106, 128, 105, 137, 123, 131, 105, 103, 
                                    102, 106, 92, 109, 97, 98, 107, 110, 103, 106, 99, 100, 96, 119, 
                                    117, 110, 95, 107, 163, 142, 158, 133, 153, 128, 128, 114, 129, 
                                    132, 160, 144, 147, 124, 123, 120, 116, 131, 150, 126, 148, 91, 
                                    101, 179, 127, 110, 81, 129, 126, 82, 141, 95, 217, 153, 124, 
                                    133, 113, 159, 99, 94, 137, 55, 75, 61, 53, 93, 41, 77, 54, 45, 
                                    62, 59, 49, 97, 99, 82, 87, 61, 68, 107, 194, 203, 271, 142, 
                                    316, 282, 121, 158, 169, 197, 144, 186, 267, 207, 215, 219, 116, 
                                    244, 155, 164, 124, 130, 63, 86, 116, 137, 52, 95, 100, 39, 138, 
                                    81, 76, 61, 143, 65, 93, 66, 104, 211, 27, 60, 55, 71, 32, 71, 
                                    87, 42, 47, 40, 78, 38, 89, 49, 1, 34, 72, 37, 38, 50, 79, 174, 
                                    220, 255, 244, 202, 156, 153, 114, 188, 176, 202, 136, 277, 192, 
                                    273, 95, 146, 96, 154)),
                  row.names = c(NA, -180L),
                  class = "data.frame")

My hypothesis is that the number of fruits given by a plant will be largely determined by that plant's height. The specific relationship between height and number of fruits, however, may change across plant types and regions. Here are some plots to illustrate what I mean: if I just plot the number of fruits versus the height of the plant, I see no correlation:

But if I look at each plant type individually, I see positive correlations:

And if I additionally split by region, I see that these correlations are variable, for instance for apple trees the correlation is stronger in the south than in the north:

Also note that my hypothesis is that the year will have no effect, however sometimes it might be the case that the apparent correlations are in fact driven by yearly variation in height and number of fruits. I have been told that I should use a GLM in order find in which cases (i.e., for which plants and in which regions) the height (and not the year) is the main determinant of the number of fruits.

I have no experience with GLMs so I am a bit lost here. I have tried using the lme4 package and running:

library(lme4)
mymodel <- lmer(n_fruits ~ plant + region + (1|height) + (1|year),
                data = data)
summary(mymodel)

Which, first of all, gives a warning that I am not sure how to interpret:

boundary (singular) fit: see help('isSingular')

I think this has to do with the random effects being too small, but I am not sure I am understanding this correctly. Then, the output is:

Linear mixed model fit by REML ['lmerMod']
Formula: n_fruits ~ plant + region + (1 | height) + (1 | year)
   Data: data
REML criterion at convergence: 1795
Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0645 -0.6951 -0.0776  0.5547  3.1917
Random effects:
 Groups   Name        Variance Std.Dev.
 height   (Intercept)  210.6   14.51

 year     (Intercept)    0.0    0.00

 Residual             1347.2   36.70

Number of obs: 180, groups:  height, 44; year, 4
Fixed effects:
            Estimate Std. Error t value
(Intercept)  180.535      7.521  24.004
plantorange  -21.498      9.229  -2.329
planttomato  -11.643      8.999  -1.294
regionnorth  -53.715      7.161  -7.501
regionsouth  -91.735      6.998 -13.109
Correlation of Fixed Effects:
            (Intr) plntrn plnttm rgnnrt
plantorange -0.583

planttomato -0.586  0.486

regionnorth -0.464 -0.002 -0.022

regionsouth -0.471  0.009 -0.010  0.520
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

And I honestly have no idea how to interpret that. Ideally I was expecting to get some metric of significance for both the height and the year, under each combination of plant type and region. My hypothesis is that plant height should be at least sometimes significant, and year should not be significant in any context. But I don't know how to extract this information, and in fact I am not even sure that this is something a GLM will give?

I would appreciate any guidance, since as you can see my lack of experience with this makes me be quite lost.

Thanks in advance!

If you want to model the association of height with n_fruits then you need to have height as a "fixed-effect" predictor in your model, not a random effect as you have coded. Also, if you want to know specific combinations of region and plant that alter the association between height and n_fruits, they also should be treated as fixed effects. Things modeled as "random effects" aren't really evaluated individually; their distribution is modeled. See this page for much guidance and discussion. — EdM, Nov 03 '22 at 20:41
it's been a while until I found this history and I am currently struggled with lmer as well. The answers help me a lot understand a better model, however, what I really want to know how can I correctly understand each row of data displayed for fixed effects in summary. For example: plantorange 1.776 53.823 0.033 I probably know that the number of fruit for orange is significant when in default height of this model. However, if I want to know if the plant of orange is significant for the whole model how could I know the answer? Or someone could correct my understanding of each row? Thanks a lot — otiliaovo, Nov 18 '23 at 22:02

score 3 · Accepted Answer · answered Nov 03 '22 at 20:55

thanks for the reproducible example.
first, this is a linear mixed-effects model (lmm), not a glm (generalized linear model). A lmm does seem appropriate, given your explanation.

based on your hypothesis, height needs to be a fixed effect. The general rule of thumb is: use fixed effects for estimating things you want to know (in this case the effect of height on number of fruits), and use random effects to account for 'nuisance parameters (i.e., remove known variation in a different covariate, so you can better estimate fixed effects).

The model you give above does not test for any effect of plant height. including plant height as a random intercept term (1|height) in effect removes variation in n_fruits with variation in plant height then compares differces among plant groups in an ANOVA-like lm fashion.

height is also confounded a little bit with plant type (there are different ranges on heights for oranges, apples, and tomatoes).

here are a few models to consider

## convert to factors
data$plant <- as.factor(data$plant)
data$region <- as.factor(data$region)
data$year <- as.factor(data$year)
fit models
this m1 model just tests for difference in fruit_n with height, controlling for region and plant type
m1 <- lmer(n_fruits~height + (1|region) + (1|plant), data = data)
m2 adds an interaction with height and plant type
m2 <- lmer(n_fruits~height*plant + (1|region), data = data)
m3 adds an fixed effect for year -- if you are interested to show that year has no effect, include it as a fixed effect then intepret its coefficient
m3 <- lmer(n_fruits~height*plant + year + (1|region), data = data)

can run summary(m#) to see the model output. To interpret those see help("pvalues",package="lme4") or many other posts here on stack overflow.

> summary(m2)
Linear mixed model fit by REML ['lmerMod']
Formula: n_fruits ~ height * plant + (1 | region)
   Data: data
REML criterion at convergence: 1776.5
Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3882 -0.7032 -0.0916  0.6210  3.4359
Random effects:
 Groups   Name        Variance Std.Dev.
 region   (Intercept) 2127     46.11

 Residual             1260     35.49

Number of obs: 180, groups:  region, 3
Fixed effects:
                   Estimate Std. Error t value
(Intercept)         -37.673     44.619  -0.844
height                6.747      1.410   4.785
plantorange           1.776     53.823   0.033
planttomato          99.017     51.443   1.925
height:plantorange   -3.875      1.617  -2.397
height:planttomato   -2.070      3.240  -0.639
Correlation of Fixed Effects:
            (Intr) height plntrn plnttm hght:plntr
height      -0.796

plantorange -0.532  0.658

planttomato -0.556  0.687  0.461

hght:plntrn  0.693 -0.871 -0.938 -0.600

hght:plnttm  0.343 -0.431 -0.286 -0.942  0.377

> summary(m3)
Linear mixed model fit by REML ['lmerMod']
Formula: n_fruits ~ height * plant + year + (1 | region)
   Data: data

REML criterion at convergence: 1759.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3877 -0.7201 -0.1014  0.6140  3.3601 

Random effects:
 Groups   Name        Variance Std.Dev.
 region   (Intercept) 2126     46.11   
 Residual             1279     35.77   
Number of obs: 180, groups:  region, 3

Fixed effects:
                   Estimate Std. Error t value
(Intercept)        -35.0741    45.4230  -0.772
height               6.6668     1.4303   4.661
plantorange         -0.6461    54.4880  -0.012
planttomato         97.3854    51.9045   1.876
yearyear_2003       -0.2350     7.5598  -0.031
yearyear_2005       -3.3337     7.5868  -0.439
yearyear_2007        1.2257     7.5550   0.162
height:plantorange  -3.7871     1.6388  -2.311
height:planttomato  -2.0203     3.2683  -0.618

Correlation of Fixed Effects:
            (Intr) height plntrn plnttm y_2003 y_2005 y_2007 hght:plntr
height      -0.797                                                     
plantorange -0.540  0.660                                              
planttomato -0.554  0.686  0.462                                       
yearyr_2003 -0.073 -0.013  0.040  0.010                                
yearyr_2005 -0.157  0.093  0.096  0.036  0.497                         
yearyr_2007 -0.104  0.026  0.045 -0.013  0.499  0.502                  
hght:plntrn  0.699 -0.872 -0.938 -0.599 -0.021 -0.104 -0.041           
hght:plnttm  0.339 -0.427 -0.284 -0.941 -0.018 -0.004  0.028  0.374

lastly, you can compare the models using anova anova(m1, m2, m3)

Data: data
Models:
m1: n_fruits ~ height + (1 | region) + (1 | plant)
m2: n_fruits ~ height * plant + (1 | region)
m3: n_fruits ~ height * plant + year + (1 | region)
   npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
m1    5 1834.1 1850.1 -912.05   1824.1                          
m2    8 1819.4 1844.9 -901.69   1803.4 20.7115  3  0.0001208 ***
m3   11 1825.0 1860.1 -901.48   1803.0  0.4124  3  0.9376655    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 3 seems most like what you were describing you wanted to see, but as the expected year is rather unimportant (t-values look non-significant). model2 is a good model. you can play around more. the model singularity error you got was likely due to including a continuous covariate (height) as a random effect and only factors as fixed effects. that will mess things up. you need the continuous varibale of height as a fixed effect the other nuisance factor groupings as random effect. year might be another random effect that would allow the model to fit better.
hope this helps.

Thanks for contributing! (+1) It's possible that the person who suggested a GLM had a Poisson model in mind, as these are count outcomes. If the real data only have 3 levels corresponding to region, it might be just as good to treat it as a fixed effect, too. Note that you have formally evaluated the values of year all at once by anova() for m3 versus m2, with p-value of 0.94. Looking at individual coefficient t-values for a categorical predictor like year in this case can be misleading, as they depend on the choice of reference category. — EdM, Nov 03 '22 at 21:18
Thank you so much for your detailed answer. This was extremely helpful! — J. Díaz-Colunga, Nov 03 '22 at 23:10
I think, though, that random slopes for height should probably also be considered. For instance, a model such as lmer(n_fruits ~ height * plant * year + (height || region:plant) ) . For the fixed effects part I am using what is possibly the largest model that the author would want to consider (which might need further simplification). The random effects part takes into account that slopes can show random variation by region and by plant, which I think makes biological sense and has the added bonus of it resulting in 9 levels of the random effects. — Ramon Diaz-Uriarte, Nov 04 '22 at 08:25
Instead of "that slopes can show random variation by region and by plant" I should have written: "that slopes can show random variation in each combination of region:plant". (This is not the same has having plant as both a random and a fixed effect; we could think of the random effect of region being allowed to differ between plant. See a somewhat similar case here https://stats.stackexchange.com/a/79392). — Ramon Diaz-Uriarte, Nov 05 '22 at 20:56

How to interpret the output of a Generalized Linear Model with R lmer

1 Answers1

fit models

this m1 model just tests for difference in fruit_n with height, controlling for region and plant type

m2 adds an interaction with height and plant type

m3 adds an fixed effect for year -- if you are interested to show that year has no effect, include it as a fixed effect then intepret its coefficient