<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/general-linear-models-the-basics/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.seascapemodels.org/rstats/2018/01/19/intro-to-glms.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h1>General Linear Models: The Basics</h1>
<p>General linear models are one of the most widely used statistical tool<br>
in the biological sciences. This may be because they are so flexible and<br>
they can address many different problems, that they provide useful<br>
outputs about statistical significance AND effect sizes, or just that<br>
they are easy to run in many common statistical packages.</p>
<p>The maths underlying General Linear Models (and Generalized linear<br>
models, which are a related but different class of model) may seem<br>
mysterious to many, but are actually pretty accessible. You would have<br>
learned the basics in high school maths.</p>
<p>We will cover some of those basics here.</p>
<h2>Linear equations</h2>
<p>As the name suggests General Linear Models rely on a linear equation,<br>
which in its basic form is simply:</p>
<p><em>y</em><sub><em>i</em></sub> = <em>α</em> + <em>β</em>x<em><sub>*i*</sub> + </em>ϵ*<sub><em>i</em></sub></p>
<p>The equation for a straight line, with some error added on.</p>
<p>If you aren’t that familiar with mathematical notation, notice a few<br>
things about this equation (I have followed standard conventions here).<br>
I used normal characters for variables (i.e. things you measure) and<br>
Greek letters for parameters, which are estimated when you fit the model<br>
to the data.</p>
<p><em>y</em><sub><em>i</em></sub> are your response data, I indexed the <em>y</em> with <em>i</em> to<br>
indicate that there are multiple observations. <em>x</em><sub><em>i</em></sub> is<br>
variously known as a covariate, predictor variable or explanatory<br>
variable. <em>α</em> is an intercept that will be estimated. <em>α</em> has the same<br>
units as y. (e.g. if y is number of animals, then <em>α</em> is expected the<br>
number of animals when x = 0).</p>
<p><em>β</em> is a slope parameter that will also be estimated. <em>β</em> is also termed<br>
the effect size because it measures the effect of x on y. <em>β</em> has units<br>
of ‘y per x’. For instance, if x is temperature, then <em>β</em> has units of<br>
number of animals per degree C. <em>β</em> thus measures how much we expect y<br>
to change if x were to increase by 1.</p>
<p>Finally, don’t forget <em>ϵ</em><sub><em>i</em></sub>, which is the error.<br>
<em>ϵ</em><sub><em>i</em></sub> will measure the distance between each prediction of<br>
<em>y</em><sub><em>i</em></sub> made by the model and the observed value of<br>
<em>y</em><sub><em>i</em></sub>.</p>
<p>These predictions will simply be calculated as:</p>
<p><em>y</em><sub><em>i</em></sub> = <em>α</em> + <em>β</em>x*<sub><em>i</em></sub></p>
<p>(notice I just removed the <em>ϵ</em><sub><em>i</em></sub> from the end). You can<br>
think of the linear predictions as: the mean or ‘expected’ value a new<br>
observation <em>y</em><sub><em>i</em></sub> would take if we only knew<br>
<em>x</em><sub><em>i</em></sub> and also as the ‘line of best fit’.</p>
<h2>Simulating ideal data for a general linear model</h2>
<p>Now we know the model, we can generate some idealized data. Hopefully<br>
this will then give you a feel for how we can fit a model to data. Open<br>
up R and we will create these parameters:</p>
<div>
<div>
<pre><code>n &lt;- 100
beta &lt;- 2.2
alpha &lt;- 30
</code></pre>
</div>
</div>
<p>Where <code>n</code> is the sample size and <code>alpha</code> and <code>beta</code> are as above.</p>
<p>We also need some covariate data, we will just generate a sequence of<br>
<code>n</code> numbers from 0 to 1:</p>
<div>
<div>
<pre><code>x &lt;- seq(0, 1, length.out = n)
</code></pre>
</div>
</div>
<p>The model’s expectation is thus this straight line:</p>
<div>
<div>
<pre><code>y_true &lt;- beta * x + alpha
plot(x, y_true)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/intro-to-glms_files/figure-markdown_strict/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<p>Because we made the model up, we can say this is the true underlying<br>
relationship. Now we will add error to it and see if we can recover that<br>
relationship with a general linear model.</p>
<p>Let’s generate some error:</p>
<div>
<div>
<pre><code>sigma &lt;- 2.4
set.seed(42)
error &lt;- rnorm(n, sd = sigma)
y_obs &lt;- y_true + error
plot(x, y_obs)
lines(x, y_true)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/intro-to-glms_files/figure-markdown_strict/unnamed-chunk-4-1.png?w=456" data-recalc-dims="1"></p>
<p>Here <code>sigma</code> is our standard deviation, which measures how much the<br>
observations y vary around the true relationship. We then used <code>rnorm</code><br>
to generate <code>n</code> random normal numbers, that we just add to our predicted<br>
line <code>y_true</code> to simulate observing this relationship.</p>
<p>Congratulations, you just created a (modelled) reality a simulated an<br>
ecologist going out and measuring that reality.</p>
<p>Note the <code>set.seed()</code> command. This just ensures the random number<br>
generator produces the same set of numbers every time it is run in R and<br>
it is good practice to use it (so your code is repeatable). Here is a<br>
great explanation of seed setting and <a rel="nofollow" href="http://www.masalmon.eu/2017/04/12/seeds/" target="_blank">why 42 is so<br>
popular</a>.</p>
<p>Also, check out the errors:</p>
<div>
<div>
<pre><code>hist(error)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/images/intro-to-glms_files/figure-markdown_strict/unnamed-chunk-5-1.png?w=456" data-recalc-dims="1"></p>
<p>Looks like a normal distribution hey? That’s because we generated them<br>
from a normal distribution. That was a handy trick, because the basic<br>
linear model assumes the errors are normally distributed (but not<br>
necessarily the raw data).</p>
<p>Also note that <code>sigma</code> is constant (e.g. it doesn’t get larger as x gets<br>
larger). That is another assumption of basic linear models called<br>
‘homogeneity of variance’.</p>
<h2>Fitting a model</h2>
<p>To fit a basic linear model in R we can use the <code>lm()</code> function:</p>
<div>
<div>
<pre><code>m1 &lt;- lm(y_obs ~ x)
</code></pre>
</div>
</div>
<p>It takes a formula argument, which simply says here that <code>y_obs</code> depends<br>
on (the tilde <code>~</code>) <code>x</code>. R will do all the number crunching to estimate<br>
the parameters now.</p>
<p>To see what it came up with try:</p>
<div>
<div>
<pre><code>coef(m1)

## (Intercept)           x
##   30.163713    2.028646
</code></pre>
</div>
</div>
<p>This command tells us the estimate of the intercept (<code>(Intercept)</code>) and<br>
the slope on x under <code>x</code>. Notice they are close to, but not exactly the<br>
same as <code>alpha</code> and <code>beta</code>. So the model has done a pretty decent job of<br>
recovering our original process. The reason the values are not identical<br>
is that we simulated someone going and measuring the real process with<br>
error (that was when we added the normal random numbers).</p>
<p>We can get slightly more details about the model fit like this:</p>
<div>
<div>
<pre><code>summary(m1)

##
## Call:
## lm(formula = y_obs ~ x)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -7.2467 -1.5884  0.1942  1.5665  5.3433
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  30.1637     0.4985  60.503   &lt;2e-16 ***
## x             2.0286     0.8613   2.355   0.0205 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.511 on 98 degrees of freedom
## Multiple R-squared:  0.05357,    Adjusted R-squared:  0.04391
## F-statistic: 5.547 on 1 and 98 DF,  p-value: 0.0205
</code></pre>
</div>
</div>
<p>I’m not going to go overboard with explaining this output now, but<br>
notice a few key things. With the summary, we get standard errors for<br>
the parameter estimates (which is a measure of how much they might<br>
vary). Also notice the R-squared, which can be handy. Finally, notice<br>
that the <code>Residual standard error</code> is close to the value we used for<br>
<code>sigma</code>, which is because it is an estimate of <code>sigma</code> from our<br>
simulated data.</p>
<p>Your homework is play around with the model and sampling process. Try<br>
change <code>alpha</code>, <code>beta</code>, <code>n</code> and <code>sigma</code>, then refit the model and see<br>
what happens.</p>
<h2>Final few points</h2>
<p>So did you do the homework? If you did, well done, you just performed a<br>
simple power analysis (in the broad sense).</p>
<p>In a more formal power analysis (which is what you might have come<br>
across previously) could systematically vary <code>n</code> or <code>beta</code> and for 1000<br>
randomised data sets and then calculate the proportion out of 1000<br>
data-sets that your p-value was ‘significant’ (e.g. less than a critical<br>
threshold like the ever-popular 0.05). This number tells you how good<br>
you are at detecting ‘real’ effects.</p>
<p>Here’s a great intro to power analysis in the broad sense: <a rel="nofollow" href="https://press.princeton.edu/titles/8709.html" target="_blank">Bolker,<br>
Ecological Models and Data in<br>
R</a></p>
<p>One more point. Remember we said above about some ‘assumptions’. Well we<br>
can check those in R quite easily:</p>
<div>
<div>
<pre><code>plot(m1, 1)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i2.wp.com/www.seascapemodels.org/images/intro-to-glms_files/figure-markdown_strict/unnamed-chunk-9-1.png?w=456" data-recalc-dims="1"></p>
<p>This shows a plot of the residuals (A.K.A. errors) versus the predicted<br>
values. We are looking for ‘heteroskedasticity’ which is a fancy way of<br>
saying the errors aren’t equal across the range of predictions (remember<br>
I said <code>sigma</code> is a constant?).</p>
<p>Another good plot:</p>
<div>
<div>
<pre><code>plot(m1, 2)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i1.wp.com/www.seascapemodels.org/images/intro-to-glms_files/figure-markdown_strict/unnamed-chunk-10-1.png?w=456" data-recalc-dims="1"></p>
<p>Here we are looking for deviations of the points from the line. Points<br>
on the line mean the errors are approximately normally distributed,<br>
which was a key assumption. Points far from the line could indicate the<br>
errors are skewed left or right, too fat in the middle, or too in the<br>
middle skinny. More on that issue<br>
<a rel="nofollow" href="http://www.seascapemodels.org/rstats/2017/10/06/qqplot-non-normal-glm.html" target="_blank">here</a></p>
<h2>The end</h2>
<p>So the basics might belie the true complexity of situations we can<br>
address with General Linear Models and their relatives Generalized<br>
Linear Models. But, just to get you excited, here are a few things you<br>
can do by adding on more terms to the right hand side of the linear<br>
equation:</p>
<ol>
<li>Model multiple, interacting covariates.</li>
<li>Include factors as covariates (instead of continuous variables). Got<br>
a factor and a continuous variable? Don’t bother with the old-school<br>
ANCOVA method, just use a linear model.</li>
<li>Include a spline to model non-linear effects (that’s a GAM).</li>
<li>Account for hierarchies in your sampling, like transects sampled<br>
within sites (that’s a mixed effects model)</li>
<li>Account for spatial or temporal dependencies.</li>
<li>Model varying error variance (e.g. when the variance increases with the mean).</li>
</ol>
<p>You can also change the left-hand side, so that it no longer assumes<br>
normality (then that’s a <strong>Generalized</strong> Linear Model). Or even add<br>
chains of models together to model pathways of cause and effect (that’s<br>
a ‘path analysis’ or ‘structural equation model’)</p>
<p>If this taster has left you keen to learn more, then check out any one<br>
of the zillion online courses or books on GLMs with R, or if you can get<br>
to Brisbane, come to our next <a rel="nofollow" href="http://www.easyr.org/courses.html" target="_blank">course (which as of writing was in Feb<br>
2018, but we do them regularly)</a>.</p>
<p>Now you know the basics, practice, practice, practice and pretty soon<br>
you will be running General Linear Models behind your back while you<br>
watch your 2 year old, which is what I do for kicks.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.seascapemodels.org/rstats/2018/01/19/intro-to-glms.html"> Bluecology blog</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/W_EU6GmVs_0" height="1"></a></div>