Objectives

The objectives of this tutorial are:

1. To introduce general linear models
   • As statistical models for describing a population
   • As a tool for estimating quantities of interest and making decisions from data
2. The concept of response and explanatory variables
3. Explain how general linear model can be divided into deterministic and random parts
4. Develop statistical models for a hypothesis (H₁) and a null-hypothesis (H₀)
5. Introduce the concept of model parameters
6. Introduce R’slm()function for fitting a general linear model to data
7. Demonstrate the use oflm()by fitting models to cortisol data in the WOLF.CSV file

Statistical modelling

Review this online lesson introducing statistical models before continuing.

Online Lesson:Statisticial Modelling

This lesson is athttps://www.ucd.ie/ecomodel/OnlineLessons/lesson4_statisticalmodelling_Website.html

A story about wolves

Every winter in Canada hunters shoot and trap wolves for their pelts. This hunting reduces the population size of wolves, but it has many other knock on consequences. For example, hunting also removes individuals from social groups which disrupts the wolf’s social hierarchy.

Hypotheses

We’re going to write our question as atestable statement, which is called a hypothesis. It is very important that this hypothesis is testable so that we can use our data to quantitatively test the hypothesis.

Start with a hypothesis about the effect of hunting on cortisol concentrations, which we’ll call thealternative hypothesis (H~1):

• H₁: Wolves in population 2 (heavily hunted) have average cortisol levels that are higher than wolves in population 1 (almost not hunted).

We’ll also write a hypothesis that is a statement about no effect, which we’ll call thenull hypothesis (H₀)

• H₀: Wolves in population 2 have the same average cortisol levels as wolves in population 1

Importing the data

First we import our data (our sample which gives partial infomration about our system). We make sure that the data type of thePopulationvariable is set as qualitative using the functionas.factor().

wolf = read.table('WOLF.CSV', header=T, sep=',')  # Import wolf data set
wolf$Population = as.factor(wolf$Population)      # Set Population as a qualitative variable

We will focus on populations 1 and 2 (light and heavy hunting respectively), so let’s subset our data to contain only these two populations.

wolf12 = droplevels(subset(wolf, Population!=3))  # Extract data from population 1 & 2

Notice that we have used thedroplevels()function to ‘update’ the levels of the factors in thewolf12data frame. Thisdroplevels()function will force R to forget about the level for population 3.

Exploratory data analysis

After importing the data it’s important to start making a picture of the data by visualising the data using numerical and graphical summaries.

summary(wolf12)        # Numerical summary of the data

   Individual         Sex            Population    Colour         
 Min.   :  1.00   Length:148         1: 45      Length:148        
 1st Qu.: 37.75   Class :character   2:103      Class :character  
 Median : 74.50   Mode  :character              Mode  :character  
 Mean   : 74.50                                                   
 3rd Qu.:111.25                                                   
 Max.   :148.00                                                   
                                                                  
     Cpgmg           Tpgmg            Ppgmg      
 Min.   : 4.75   Min.   : 3.250   Min.   :12.76  
 1st Qu.:12.16   1st Qu.: 4.378   1st Qu.:19.50  
 Median :15.38   Median : 5.030   Median :25.00  
 Mean   :16.61   Mean   : 5.617   Mean   :25.89  
 3rd Qu.:19.98   3rd Qu.: 6.067   3rd Qu.:30.01  
 Max.   :40.43   Max.   :15.130   Max.   :53.28  
                                  NA's   :79     

# Calculate mean Cortisol in each population
aggregate(Cpgmg~Population, data=wolf12, FUN=mean)  

  Population    Cpgmg
1          1 15.56222
2          2 17.07495

Graphical data exploration

Exploring the data with graphics is also useful

library(ggplot2)       # Load ggplot2

First we’ll look at thedistributionof the cortisol concentrations

# Plot the distribution of the cortisol data
ggplot(data=wolf12,
       aes(x=Cpgmg)) +
  geom_histogram(bins=20) +
  labs(x='Cortisol Concentration [pg/mg]') + 
  theme_bw()

The distribution looks slightly right-skewed.

We can plot the raw data

# Plot the raw cortisol data
ggplot(data=wolf12,
       aes(x=Population,
           y=Cpgmg,
           colour=Population)) +
  geom_point(position='jitter') +
  labs(x='Population',
       y='Cortisol Concentration [pg/mg]') + 
  theme_bw()

or summarise the data with a boxplot. Importantly, The boxplot not only visualises the central tendency of the data from the two populations (the median) but also the spread (in terms of the 25% and 75% quantiles as well as the whiskers).

# Draw a boxplot of the cortisol data from populations 1 and 2
ggplot(data=wolf12,
       aes(x=Population,
           y=Cpgmg)) +
  geom_boxplot() +
  labs(x='Population',
       y='Cortisol Concentration [pg/mg]') + 
  theme_bw()

Median cortisol levels in population 2 are slightly higher than population 1, but there is a lot of spread about the central tendency.

A violin plot is an alternative to a boxplot

# Draw a violin plot of the cortisol data from populations 1 and 2
ggplot(data=wolf12,
       aes(x=Population,
           y=Cpgmg)) +
  geom_violin() +
  labs(x='Population',
       y='Cortisol Concentration [pg/mg]') + 
  theme_bw()

We may want to look at how several variables are able to explain cortisol values (the response variable). For example, differences between males and females and differences between the two populations.

# Draw a violin plot of the cortisol data from populations 1 and 2
# for males and females
ggplot(data=wolf12,
       aes(x=Population,
           y=Cpgmg,
           fill=Sex)) +
  geom_violin() +
  labs(x='Population',
       y='Cortisol Concentration [pg/mg]') + 
  theme_bw()

Balanced and unbalanced data

Use thetable()function to look at how much data we have from each population.

# Look at the balance of the data across
# levels of Population
table(wolf12$Population)

  1   2 
 45 103 

The sample size for population 1 is lower than population 2. This is called anunbalanced designand is very common when working with data collected from the field. Abalanced designis when sample sizes are the same for all levels of a factor. Classical ANOVA analysis required the data to be balanced. So these data would require special treatment if analysed with a classical ANOVA. General linear models can be used with unbalanced data, because they use maximum likelihood to fit the model. However, general linear model analysis is easiest for balanced data sets, and they should not be used for data that are extremely unbalanced.

Analyse the data

Modelling the population

Our data are asample. The data contain partial information about the topic of ourhypotheses(i.e. thepopulation). The population is the set of cortisol concentrations from all wolves from the two locations, at any time of the year, and maybe even across more of ta wolf’s body than just the hair.

If we are going to investigate our hypotheses then we need a way to describe the population. We will use a statistical model to describe the distribution of data we may expect across the whole population.

A statistical model

The figure below sketches out a possible statistical model. It sketches the distribution of cortisol for populations 1 and 2. The model has three important aspects:

1. Theresponse variableis cortisol concentration (on the x-axis), and it is a quantitative variable
2. The average value of cortisol in each population is described by a number (represented as μ1​and μ2​, whose value is as yet not known). In other words, the average value in each population has a unique value than can bedetermined.Given complete infomration about the population we could determine the value precisely.
3. The variation of cortisol values in each population is being described by anormal distribution. The variation within a population is being described as beingrandombecause we are not trying to specify unique values of cortisol for every individual wolf. Instead we are describing all the possible outcomes of cortisol values by a normal distribution. There is an important differnce to the last point:Given complete infomration about the population we could not determine the precise value of a wolf’s cortisol value, but we can assign a probability to observing a cortisol concentration.

This model in the above sketch could be written in equation form as

Cpgmg=μi​+Normal(0,σ)

where μi​is the mean cortisol for population $i$($i$is an index that can be 1 or 2) and $\sigma$is the standard deviation in cortisol values (i.e. the spread in cortisol values).

Note that:

• the spread is the same in both population 1 and 2.
• $\text{Normal}(0,\sigma )$is shorthand for a normal distribution with meanof zero and standard deviation $\sigma$.

General linear models

The sketch above is a representation of a general linear model. General linear models are a hugely important family of statistical models.

A general linear model always uses a normal distribution to describe variation in the population that we will not try to associate with another variable (e.g. variation in cortisol concentrations between wolves in a population). The normal distribution in a general linear model is assumed to have a mean of zero and a constant standard deviation that must be estimated from data.

Special cases of general linear models are the methods of: linear regression, t-tests, analysis of variance (ANOVA) and analysis of covariance (ANCOVA).

General linear models are a statistical tool that allow us to:

1. Analyse a very wide range of data with one statistical methodology
2. Always describe randomness using a normal distribution
3. Estimatemeans and other parameters of the population (e.g. the difference in mean cortisol between populations)
4. Estimate the uncertaintyin population parameters from point 3.
5. Quantify the strength of evidenceagainst the null-hypothesis

The work-horse function in R for fitting general linear models to data is thelm()function.

We will use the general linear model framework to estimate mean cortisol values in populations 1 & 2, and then test the hypothesis that these means are the same. We will then compare this approach to a t-test approach and finally to a permutation test approach.

Hypotheses revisited

Each hypothesis (the null-hypothesis, H₀and the alternative hypothesis, H₁) can be expressed in terms of our statistical model.

Fitting a model to data

Use thelm()function to fit a general linear model to data.

To use thelm()function we specify the model using aformula(see above) and specify the data frame containing the data.

Our formulas are

• H₀model:Cpgmg ~ 1
• H₁model:Cpgmg ~ 1 + Population

Theresponse(Cpgmg) is on the left hand side of the tilde (~) and thedeterministic partof the model (theexplanatory variables) is on the right hand side of the tilde. Therandom partis not represented in the formula because it is automatically taken to be a normal distribution. The functionlm()always assumes the random part of a model is a normal distribution.

Here is the R code to fit the H₀model (null model) and the H₁model (alternative model)

# Fit a linear model for the null-model 
# (i.e. Population is not part of the model)
m0 = lm(Cpgmg ~ 1, data=wolf12)



# Fit a linear model for the alternative model 
# (i.e. Population affects cortisol)
m1 = lm(Cpgmg ~ 1 + Population, data=wolf12)

summary(m0)        # Summarise the null model

Call:
lm(formula = Cpgmg ~ 1, data = wolf12)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.865  -4.455  -1.240   3.360  23.815 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  16.6150     0.5051    32.9   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.145 on 147 degrees of freedom

# Model coefficient estimates and uncertainties

coef(m0)      # Coefficients for model m0

(Intercept) 
     16.615 

confint(m0)   # 95% confidence intervals on the coefficients for model m0

               2.5 %   97.5 %
(Intercept) 15.61685 17.61315

coef(m1)      # Coefficients for model m1

(Intercept) Population2 
  15.562222    1.512729 

confint(m1)   # 95% confidence intervals on the coefficients for model m1

                 2.5 %    97.5 %
(Intercept) 13.7575212 17.366923
Population2 -0.6505745  3.676033

Residuals

A residual is the difference between a data point and the prediction from a model. The residual is therefore associated with the random part of the model. It is the part of the response that cannot be explained by the model.

Most of the assumptions in our statistical model concern the residuals. So we will look at residuals in more detail when we discuss if the fitted model obeys the assumptions (model validation is the subject of a later tutorial).

The residuals of a fitted linear model can be obtained with theresiduals()function

m1_resid = residuals(m1)   # Extract the residuals for model m1

The distribution of residuals has a mean of zero,

mean(m1_resid)             # Mean of the residuals from model m1

[1] 5.521109e-16

The calculated mean is extremely close to zero. It is not exactly zero because of the limited accuracy of the computer (remember1.0e-16is computer shorthand for 1x10^-16= 0.0000000000000001).

The distribution of residuals can also be visualised by plotting a histogram

# Plot the distribution of residuals from model m1
ggplot(data=as.data.frame(m1_resid),
       aes(x=m1_resid)) +
  geom_histogram(bins=30) +
  theme_bw()

sd(m1_resid)              # Standard deviation of the residuals

[1] 6.104729

m1_resid2 = wolf12$logCpgmg - predict(m1)  # Calculate residuals by hand

Summary of the topics covered

• The statistical model for general linear models revolves around the Normal distribution
• General linear models always use a normal distribution as the random part
• The variable of interest is called the response variable
• Statistical models of data always have a deterministic and a random part
• Statistical models have parameters
• The values of a model’s parameters can be estimated by fitting a model to data
• Thelm()function fits general linear models to data

Further Reading

All these books can be found in UCD’s library

• (opens in a new window)Cohen, J. (1990). Things I have learned (so far). American Psychologist, 45(12), 1304–1312.
• (opens in a new window)Cohen, J. (1994). The earth is round (p < .05). American Psychologist, 49(12), 997–1003.
• Andrew P. Beckerman and Owen L. Petchey, 2012Getting Started with R: An introduction for biologists(Oxford University Press, Oxford) [Chapter 5]
• Michael J. Crawley, 2015Statistics : an introduction using R(John Wiley & Sons, Chichester) [Chapters 7, 8]
• Tenko Raykov and George A Marcoulides, 2013Basic statistics: an introduction with R(Rowman and Littlefield, Plymouth)