# GLMM: Elevation as numerical or factor in model?

I want to model the effect of treatments (fragmented and continuous forest) on species abundance within three different elevational regions in R, using the lme4 package. I'm also interested in differences between elevations, so I included this term in the model as a fixed effect:

model <- glmer(abundance ~ (1|site) + treatment + elevation + treatment*elevation, data=df, family= "poisson)


The three elevations are at 300m, 800m and 1200m, and several kilometre apart from each other. It seems that I have two options of treating elevation:

1. Continuous variable (300, 800, 1200)
2. Factor ("300m", "800m", "1200m")

Both give different results in the models. Which one should be preferred? What do I have to consider when deciding for one of the options?

Thanks!

### Edit

Here is the data frame that I am working with. I have two treatments per elevation, resulting in six sites

         date      chao    chao.se   nSpec  treatment region field.season elevation site count
1  2016-03-06 12.333333  3.9907300       9 fragmented      M         2016       800   CM    19
2  2016-03-07  6.777778  3.3554820       5 fragmented      M         2016       800   CM    17
3  2016-03-08  7.571429  3.7923580       5 continuous      M         2016       800   SM    12
4  2016-03-10 16.090909  9.2056194       7 fragmented      A         2016      1200   CA    21
5  2016-03-11  9.200000  4.3081318       6 continuous      A         2016      1200   SA    10
6  2016-03-12 27.000000 21.0396451      11 continuous      A         2016      1200   SA    29
7  2016-03-14  7.600000  1.2000000       7 continuous      B         2016       300   SB    15
8  2016-03-15 11.076923  3.1599378       9 fragmented      B         2016       300   CB    30
9  2016-03-19  8.615385  1.2243422       8 continuous      M         2016       800   SM    22
10 2016-03-20  8.062500  3.1402752       6 fragmented      M         2016       800   CM    17
11 2016-03-21 20.458333 15.7412343       9 continuous      A         2016      1200   SA    20
12 2016-03-22 16.384615 10.7911864       9 fragmented      A         2016      1200   CA    28
13 2016-03-26 11.384615  2.1400658      10 fragmented      B         2016       300   CB    33
14 2016-03-27  8.312500  2.0453835       7 continuous      B         2016       300   SB    13
15 2016-12-15 27.615385 21.8342329      11 fragmented      A         2017      1200   CA    38
16 2016-12-16  5.909091  2.0170916       5 fragmented      A         2017      1200   CA    22
17 2017-02-02  7.285714  2.0101782       6 fragmented      M         2017       800   CM    10
18 2017-02-03 17.937500 15.0422582       7 continuous      M         2017       800   SM    10
19 2017-02-04 10.142857  6.0421058       5 continuous      M         2017       800   SM    10
20 2017-02-05 18.000000  9.0921211      10 fragmented      M         2017       800   CM    17
21 2017-02-06 19.250000 15.4616461       8 continuous      A         2017      1200   SA    16
22 2017-02-07  4.208333  0.6353313       4 continuous      A         2017      1200   SA    10
23 2017-02-13 12.733333  4.9638695       9 continuous      B         2017       300   SB    20
24 2017-02-14  9.200000  4.3081318       6 continuous      B         2017       300   SB     8
25 2017-02-15 11.400000  6.3673096       6 fragmented      B         2017       300   CB    14
26 2017-02-16 14.769231  7.0027467       9 fragmented      B         2017       300   CB    18
27 2017-02-17  3.200000  0.6164414       3 fragmented      A         2017      1200   CA     6
28 2017-02-18 10.370370  3.1403680       8 continuous      A         2017      1200   SA    14
29 2017-02-27  9.285714  3.0416608       7 fragmented      M         2017       800   CM    11
30 2017-02-28 13.111111 10.4041777       6 continuous      M         2017       800   SM    11
31 2017-04-02 29.437500 27.2692736       8 continuous      B         2017       300   SB    10


The core part of the question isn't specific to mixed models; it would apply to any (generalized) linear modeling approach. However, the small number of sites, and the potential confounding of site with elevation, could raise some more MM-specific issues (see below).

## Part 1 (general)

• If you treat elevation as a continuous variable, you are assuming a strictly linear relationship between (log) abundance and elevation, i.e. the change in log-abundance between 300 and 800 meters is exactly 5/4 of the difference from 800 to 1200 (from the ratio of the differences, (800-300)/(1200-800)).
• If you treat elevation as a factor, you are essentially assuming independent rates of change from 300 to 800 and from 800 to 1200. Depending on the contrasts you use, your parameters can be interpreted in different ways: with the default treatment contrasts, you will get estimates of the difference between 300 and 800 and between 800 and 1200 m, while "successive differences" contrasts (MASS:contr.sdif) will give estimates between 300 and 800 and between 800 and 1200. If you use orthogonal polynomial contrasts (which will happen by default if you make elevation into an ordered factor, ?ordered), you will get one parameter corresponding to linear change and another encapsulating quadratic changes.

In general I would recommend making elevation a factor, as assuming that the change per unit elevation is always the same seems strong.

## Part 2 (mixed models)

updated after clarification from OP

If you (or some future reader) have only one site per elevation, or one site per elevation-treatment combination (the OP has two treatments, each conducted at each elevation, but only one site per treatment-elevation combo) then there are further problems which are mixed-model related:

• in this case, treating elevation as a factor will make it (or the treatment:elevation interaction) completely redundant with the site random effect (knowing the elevation/elevation+treatment of an observation is equivalent to knowing what site it occurs in)
• if you only have three sites, it won't really be practical to treat site as a random effect.
• more fundamentally, in your design the site effect and the elevation effect (or the elevation-treatment combo) are confounded. It's almost impossible to ascribe changes among sites to elevation in a watertight way: "elevation-related" changes could just be site effects.

In that case, my suggestion would be drop the site effect (and fit a regular generalized linear model), treat elevation as a (fixed) factor, but be completely upfront in your presentation that your design is pseudoreplicated, and you can only tentatively ascribe these effects to the influence of elevation/treatment (and go read, or re-read, Hurlbert 1984 ...)

• Thank you for the extensive answer! It solved my question about the elevation. Regarding Part 2 of your answer: I have two sites per elevation (1 per treatment), each consisting of 4-6 sampling days. I initially included site as a random effect, to account for the replicated sampling at the sites. Is the model still fundamentally problematic in that context? May 29 '17 at 16:48
• no, that sounds reasonable. May 29 '17 at 18:31
• actually, I take it back. If you include a treatment-by-elevation interaction in your model (and elevation is treated as a factor), you're back to the same problem, since you have one fixed-effect parameter per site ... May 29 '17 at 18:35
• see edits ...... May 29 '17 at 18:55
• Thank you for the clarification. The pseudoreplicates would be the sampling days per site, which I'm not accounting for anymore, after removing the random effect? May 29 '17 at 19:01