Generalised Linear Model
In a previous module, we showed you had to create a simple ANOVA model using R.Using anova we created a model to assess how yields varied amongst treatments within an experimental trial. However, we have identified the assumption that all of our observations are independent may not be reasonable within our design.
This is very common in agricultural analyses which will have specific layouts or involve operating with multiple farmers - this will invalidate the assumption that all observations are independent. In these sort of cases we need to move from a "general linear model" into a "generalised linear model" which can account for these structures. This is also known as a mixed model since it contains "fixed effects" (variables we are interested in) and "random" effects (variables outlining the structure or clustering of observations in the data) You can some materials on mixed models in the following two links, we recommend the latter if you wish to see a more agricultural focused example
Mathematically this is quite a bit more complicated - but thankfully in R it is not very different at all to move towards this model, you will see many of the same functions as in the previous sections. We will use the same libraries we have seen so far in other parts of this module, but now also add the lmertest library to be able to fit a "linear mixed effects regression" (lmer()) model. This approach is incredibly flexible and can let us take into account any of the standard agricultural experimental designs as long as we are able to code the random effects appropriately.
But it also is the exact same approach we need for data which comes from less structured designs - from farm trials and surveys - where we have different forms of clustering (farmers within the same village) or pseudo-replication (the same farmer with multiple plots). It can be tricky to work out exactly how to structure the code for your model, don't be afraid to ask for help if you are unsure! There are resources at the end of this document which cover many of the more common agricultural designs - but the functions used will all be the same. This even applies to cases where we start to use non-Normal models - such as Poisson regression, for count data, or logistic regression, for binary data. In these cases we use again more or less the same syntax, and all the same functions, but now using the glmer() function, which can account for these modifications to the type of variable we use as a response.
It does not matter at all if we have continuous predictor variables, or categorical predictor variables, or a mixture of the two. All these models are fitted and interpreted in R in more or less exactly the same way.
You can check the next two pages if you wish to remind yourself of the data and packages we have been using in this statistical analysis course.
Once done with this tutorial, you can move onto our second set of statistical analysis exercises
Data
The data used in this example is from a study was conducted in Eastern Zambia with the main objective of improving the efficiency of the natural fallows by using appropriate trees that may have relevance in soil fertility regeneration within permissible fallow periods.
The data has been read into R and assigned to an object called fallow.
The design was a RCBD experiment (randomised complete block design) with 4 blocks and 9 treatments. If you take a look at the data below you will see that we have data concerning the measurement of the yield (the column yield within the data) and the striga count (striga) of a plot (the column plot), within a block (the column rep), after a specific treatment (the column treat).
Packages used
In this tutorial we will use four new packages: emmeans, multcomp, lmerTest and MuMIn; in addition we will use some of the packages we have used regularly so far in this course dplyr and ggplot2. Make sure these are all installed and loaded before running any of the code on your own computer.
library(ggplot2)
library(dplyr)
library(lmerTest)
library(multcomp)
library(multcompView)
library(emmeans)
library(MuMIn)
Fitting the model
In this design, an RCBD (randomised complete block design), we have one treatment factor, "treat", and one layout factor "rep". If we want to take into account this repeated blocks layout structure, we can’t use the “lm” function anymore and instead we need to include the “rep” blocking factor into the specification of a linear mixed effect model. Using the lmer() function to include a blocking factors, or "random effects" as they are known in statistics, we use a slightly different style of code. We use brackets and then within the brackets we write 1|variablename , so in this case (1|rep). Then the rest of the model code is the same as before - this time I will store as an object called gen_lin_mod.
gen_lin_mod<-lmer(yield~treat+(1|rep),data=fallow)Again - to get anything useful out of the model we need to use additional functions - we can start with summary().
summary(gen_lin_mod)The output looks a little different after we summarise a mixed effects model, although most things remain very similar. The table of coefficients is the same as before and has the same interpretation. The additional output showing the correlation between fixed effects is not very useful for us here, as in this designed experiment all treatments are replicated the same number of times.
The residual error term is still presented but in a different place - this value is 0.4843 which can be found in the "Random Effects" table. This sort of model breaks down the variance into that which can be explained by our design factors (the rep term) and that which is left over as residual variance. Comparing the amount of the variance explained by each component can be useful to give us an understanding of how much our design factors are influencing our outcome variable. In this case the amount of variance explained by the rep (standard deviation = 0.1175) is fairly small in relation to the residual variance (standard deviation = 0.4843).
One thing which is missing when we compared to output of running summary after an lm is the R-squared value. This can be obtained by using the function r.squaredGLMM which comes from the MuMIn library.
r.squaredGLMM(gen_lin_mod)Two versions of R-squared as shown - the marginal R square R2m- which shows only the variance explained by the "fixed effects" in this case the treatment term. The conditional R-squared R2c shows the effect of both the fixed and random effects - so includes both rep and treat. So we can see 81% of the variability in yield is explained by the different treatments, and this increases to 82.3% is we include rep as well.
To obtain an ANOVA table - we can obtain this as before using anova(model)
anova(gen_lin_mod)The ANOVA table only includes the fixed effects - showing us a highly significant p-value for treatment, providing evidence of some overall differences between our treatments. We could also test whether our random effects are significant using ranova.
ranova(gen_lin_mod)This does not show us a significant effect of the different reps we have in the data. However because this variable is describing the structure of our data, this is usually not something we are particularly interested in. With structural variables, we generally want to account for any impact they may have on the analysis whether or not they could be considered to be statistically significant. But to help your understanding, note that this non-significance would suggest that the results we got from the last section using the “lm” function was probably somewhat reliable, despite not having taken into account the layout structure in blocks.
Check the model
We have many of the same assumptions to check with this model as well, although the functions here work a little differently to before.
The function plot() will only plot the fitted values from the model against the expected values.
plot(gen_lin_mod)We can check most of the key assumptions from this one plot. The residual Vs fitted plot is a scatter plot of the Residuals on the y-axis and the fitted on the x-axis and the aim for this plot is to test the assumption of equal variance of the residuals across the range of fitted values. Since the residuals do not funnel out (to form triangular/diamond shape) the assumption of equal variance is met.
We can also see that there are no extreme values in the residuals which might be potentially causing problems with the validity of our conclusions (leverage)
To assess the assumption of approximate normality we can produce the same plot as before if we extract the residuals and use the function qqnorm. This shows us how closely the residuals follow a normal distribution - if there are severe and systematic deviations from the line then we may want to consider an alternative distribution. The most elegant way of doing this is to use pipes - firstly to obtain the residuals and then to pipe those residuals into the qqnorm function..
gen_lin_mod %>%
resid() %>%
qqnorm() In this case the residuals seem to fit the assumption required for normality.
Post-hoc tests
The functions we saw in the previous section, emmip and emmeans and cld all work in exactly the same way. All we need to change is the input model object, and we obtain the plots with confidence intervals and the mean separation analysis.
emmip(gen_lin_mod,~treat,CIs = TRUE)emmeans(gen_lin_mod, ~treat) %>%
cld(Letters=letters)In the output, groups sharing a letter in the .group are not statistically different from each other as before.
Using the example from the previous section, see if you can write the code to make a plot showing the results from the lmer model including the letter displays and the error bars
You can see the data visualisation modules here for a reminder
emmeans(gen_lin_mod, ~treat) %>%
cld(Letters=letters) %>%
data.frame() %>%
ggplot(aes(y=emmean,x=treat,ymin=lower.CL,ymax=upper.CL,label=.group))+
geom_point()+
geom_errorbar()+
geom_text(nudge_x=-0.2)Model choices
Ultimately our model fitted using a generalised linear mixed model, with a random effect to account for block, using lmer provided almost identical results to the standard linear model fitted using lm. As seen below.
summary(gen_lin_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yield ~ treat + (1 | rep)
## Data: fallow
##
## REML criterion at convergence: 51.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.98578 -0.40272 -0.00098 0.47492 1.88450
##
## Random effects:
## Groups Name Variance Std.Dev.
## rep (Intercept) 0.01382 0.1175
## Residual 0.23458 0.4843
## Number of obs: 36, groups: rep, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.5475 0.2492 26.3478 22.262 < 2e-16 ***
## treat2 G.sepium -1.7025 0.3425 24.0000 -4.971 4.47e-05 ***
## treat3 L.leuco -1.8725 0.3425 24.0000 -5.468 1.28e-05 ***
## treat4 F.congesta -2.0875 0.3425 24.0000 -6.095 2.70e-06 ***
## treat5 C.siamea -3.4425 0.3425 24.0000 -10.052 4.44e-10 ***
## treat6 C.calo -2.9600 0.3425 24.0000 -8.643 7.82e-09 ***
## treat7 groundnut -2.4650 0.3425 24.0000 -7.198 1.94e-07 ***
## treat8 maize -3.5725 0.3425 24.0000 -10.431 2.14e-10 ***
## treat9 nat fallow -2.7375 0.3425 24.0000 -7.993 3.20e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trt2G. trt3L. trt4F. trt5C. trt6C. trt7gr trt8mz
## treat2G.spm -0.687
## treat3L.lec -0.687 0.500
## trt4F.cngst -0.687 0.500 0.500
## treat5C.sim -0.687 0.500 0.500 0.500
## treat6C.cal -0.687 0.500 0.500 0.500 0.500
## tret7grndnt -0.687 0.500 0.500 0.500 0.500 0.500
## treat8maize -0.687 0.500 0.500 0.500 0.500 0.500 0.500
## tret9ntfllw -0.687 0.500 0.500 0.500 0.500 0.500 0.500 0.500summary(anova_mod)##
## Call:
## lm(formula = yield ~ treat, data = fallow)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02250 -0.25000 0.02625 0.21125 0.99250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5475 0.2492 22.262 < 2e-16 ***
## treat2 G.sepium -1.7025 0.3524 -4.831 4.80e-05 ***
## treat3 L.leuco -1.8725 0.3524 -5.313 1.31e-05 ***
## treat4 F.congesta -2.0875 0.3524 -5.923 2.59e-06 ***
## treat5 C.siamea -3.4425 0.3524 -9.768 2.34e-10 ***
## treat6 C.calo -2.9600 0.3524 -8.399 5.20e-09 ***
## treat7 groundnut -2.4650 0.3524 -6.995 1.61e-07 ***
## treat8 maize -3.5725 0.3524 -10.137 1.06e-10 ***
## treat9 nat fallow -2.7375 0.3524 -7.768 2.36e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4984 on 27 degrees of freedom
## Multiple R-squared: 0.8493, Adjusted R-squared: 0.8047
## F-statistic: 19.03 on 8 and 27 DF, p-value: 3.144e-09However, this is because the blocking in this experiment had very little impact on the results - it was an on-station trial with very controlled conditions and similar environments across the four blocks.
In increasingly complex situations, particularly those with farmers conducting multiple treatments or with clusters of farmers in very different situations, very different results would be obtained from these two modelling approaches. Plots on one farmers field will be similar to other plots on the same field but likely extremely different to another farmer's plots.
In cases where very different results are found then it is usually worth trying to understand why! You will see an example of this in the exercises.
Methodological Principles
There are always many different ways of doing all that we have done here in R. If you look for guidance online about how to analyse particular designs you may find lots of different options for analysing "split plot" designs or "survey" data. It is always a good idea to ask for help when starting to analyse data if you are not sure about the way it should be analysed - you should make sure however you analyse the data it is done in a way which lets you address the questions that are most relevant to you. This is the big advantage of using the approach shown in this tutorial of the lmer and lm functions - they are incredibly flexible and can be used in many different situations. For instance, we demonstrated fitting our model as a linear mixed effect model rather than traditional ANOVA. Although the analysis in this instance was not effected too much by this change the lmer model has a number of other advantages as well:
- They are very flexible especially where we have repeated measures, for instance you don’t need to have the same number of observations per subject/treatment.
- Ability to account for a series of random effects. Not only are farms/farmers/plots…. different from each other, but things with in farms/plots….. also differ . Not taking these sources of variation into account will lead to underestimations of accuracy.
- Allows for generalization of non-normal data.
- Handling missing data: If the percentage of missing data is small and that data missing is a random sample of the data set,data from the observations with missing data can be analysed with
lmer(unlike other packages that would do list-wise deletion). - Takes into account variation that is explained by the predictor variables of interest i.e. fixed effects and variation that is not explained by these predictors i.e. random effects.