<jupyter_start><jupyter_text>Linear Models: Regression  ContentsIntroductionChapter aimsExploring the dataCorrelationsTransformations and allometric scalingPerforming the Regression analysisExercisePlotting the fitted Regression modelModel diagnosticsReporting the model  IntroductionLienar regression analysis, that is, fitting a linear model of the regression variety is very often an good choice if both, your response (dependent) and your predictor (independent) variables are continuous.  Chapter aimsIn this chapter, you will learn:* More R functions for plotting data and models to visualize your data* To calculate correlation coefficients to explore your data* To fit a linear regression model to data and test its significance* To use diagnostic plots to assess whether the lienar model is appropriate for your dataAs with the previous chapters, we'll start with creating a new blank script for you to fill in during the practical (the script file `regress.R` is also available from TheMulQuaBio). We'll also be using the genome size data again, so:$\star$ Open R and `setwd` to your `code` directory.$\star$ Create a new blank script called `Regression.R` and add some introductory comments.$\star$ Add code to your script to load the genome size data into R and check it:<jupyter_code>genome <- read.csv('../data/GenomeSize.csv')

head(genome)<jupyter_output><empty_output><jupyter_text>Exploring the dataIn previous chapters we used `plot` to create a scatterplot between two variables. If you have a set of variables to explore, writing code for each plot is tiresome, so R provides a the function `pairs`, which creates a grid of scatterplots between eachpair of variables. All it needs is a dataset.$\star$ Add `pairs(genome, col=genome$Suborder)` into your script and run the code:<jupyter_code>pairs(genome, col=genome$Suborder)<jupyter_output><empty_output><jupyter_text>The result is messy! There are far too many variables in `genome` for this to be useful. We need to cut down the data to fewer variables. [Previously](12-ExpDesign.ipynbExplore-further-by-scatter-plotting-two-variables) we used indices to select colours; here, we can use indices to select columns from the data frame. This again uses square brackets (`x[]`), but a data frame has two dimensions, rows and columns, so you need to provide an index for each dimension, separated by commas. If an index is left blank, then all of that dimension (i.e. all rows or columns) are selected. Try the following to re-acquaint yourself to access data frame content using indices:<jupyter_code># create a small data frame:
dat <- data.frame(A = c("a", "b", "c", "d", "e"), B = c(1, 2, 3, 4, 5))
dat[1, ] # select row 1 (all columns selected)
dat[, 2] # select column 2 (all rows selected)
dat[2, 1] # select row 2, column 1<jupyter_output><empty_output><jupyter_text>Now let's get started with the actual analysis. We will look at five key variables: genome size, body weight, total length, forewing length and forewing area. If you look at the output of `str(genome)`, you'll see that these are in columns 4, 7, 8, 12 and 14. We can record the indices of these columns and use this to select the data in the pairs plot:<jupyter_code>morpho_vars <- c(4, 7, 8, 12, 14) # store the indices
pairs(genome[, morpho_vars], col = genome$Suborder)<jupyter_output><empty_output><jupyter_text>$\star$ Add the code above to your script and run it.In the figure above, each scatterplot is shown twice, with the variables swapping between the $x$ and $y$ axes. You can see immediately that the relationships between the four morphological measurements and genome size are fairly scattered but that the plots comparing morphology show much clearer relationships. CorrelationsOne way of summarising how close strong the relationship between these variables are is to calculate a correlation coefficient. Pearson correlations look at the difference of each point from the mean of each variable (and since it uses means, it is a parametric statistic).It is calculated using of the differences from the mean on each axis. The key calculation is — for each point – to get the product of the differences on each axis and add them up. If the points are mostly top left ($-x$, $y$) or bottom right ($x$, $-y$) then these products are mostly negative ($-xy$); if the points are mostly top right ($x$, $y$) or bottom left ($-x$, $-y$) then the products are mostly positive ($xy$).---                                 Illustration of what the Pearson correlation coefficient means.                        ---The plots above show three clear cases where all the values of $xy$ are negative or positive or where both are present and sum to zero. The Pearson correlation coefficient simply scales these sums of $xy$ to be between -1 (perfectly negatively correlated) and 1 (perfectly positively correlated) via zero (no correlation).We will use two functions to look at correlations. The first is `cor`, which can calculate correlations between pairs of variables, so is a good partner for `pairs` plots. The second is `cor.test`, which can only compare a single pair of variables, but uses a $t$ test to assess whether the correlation is significant.$\star$ Try the following (and include it in your R script file):<jupyter_code>cor(genome[, morpho], use = "pairwise")<jupyter_output><empty_output><jupyter_text>This is the correlation matrix. Then:<jupyter_code>cor.test(genome$GenomeSize, genome$TotalLength, use = "pairwise")<jupyter_output><empty_output><jupyter_text>The `use='pairwise'` tells R to omit observations with missing data and use complete pairs of observations. The first function confirms our impressions from the graphs: thecorrelations between genome size and morphology are positive but comparatively weak and the correlations between morphological measurements are positive and very strong (i.e. close to 1). The correlation test tells us that genome size and body length are positively correlated (r=0.34, $t$ = 3.5507, df = 96, $p$ = 0.0006). *Remember this example when reporting correlations in your reports &ndash; all this information needs to be included* Transformations and allometric scalingThere is one problem with the correlations above: *the correlation coefficient calculation assumes a straight line relationship*. Some of the scatterplots above are fairly straight but there are some strongly curved relationships. This is due to [allometric scaling](https://en.wikipedia.org/wiki/Allometry), where one body measure changes (or grows) disproportionately with respect to another. Here, two of the variables are in linear units (total and forewing length), one is in squared units (forewing area) and one in cubic units (body weight, which is approximately volume). That these measures are in different units itself guarantees that they will scale allometrically with respect to each other. The relationships between these variables can be described using a power law: $$y = ax^b$$Fortunately, if we log transform this equation, we get $\log(y) = \log(a) + b \log(x)$. This is the equation of a straight line ($y=a+bx$), so we should be able to make these plots straighter by logging both axes. We saw [previously](13-t_F_tests.ipynb$t$-tests-revisited) that we can create a new logged variable in the data frame like this:<jupyter_code>genome$logGS <- log(genome$GenomeSize)<jupyter_output><empty_output><jupyter_text>$\star$ Using this command as a template, create a new logged version of the fivevariables listed above:<jupyter_code>genome$logGS <- log(genome$GenomeSize)
genome$logBW <- log(genome$BodyWeight)
genome$logTL <- log(genome$TotalLength)
genome$logFL <- log(genome$ForewingLength)
genome$logFA <- log(genome$ForewingArea)<jupyter_output><empty_output><jupyter_text>$\star$ Then, using `str`, work out which column numbers the logged variables are and create a new variable called `logmorpho` containing these numbers:<jupyter_code>str(genome)
logmorpho <- c(17,18,19,20,21)<jupyter_output><empty_output><jupyter_text>We can now use the `pairs` and `cor` test as beforefor the columns in `logmorpho`:<jupyter_code>pairs(genome[, logmorpho], col=genome$Suborder)
cor(genome[, logmorpho], use='pairwise')<jupyter_output><empty_output><jupyter_text>The scatterplots show that logging the data has very successfully addressed curvature (non-linearity) due to allometric scaling between variables in the data. Performing the Regression analysisWe'll now look at fitting the first linear model of this course, to explore whether log genome size explains log body weight. The first thing to do is to plot the data (as you did [previously](12-ExpDesign.ipynbExplore-further-by-scatter-plotting-two-variables)):<jupyter_code>myColours <- c('red', 'blue') # Choose two colours
mySymbols <- c(1,3) # And two different markers
colnames(genome)
plot(logBW ~ GenomeSize , data = genome, #Now plot again
col = myColours[Suborder], pch = mySymbols[Suborder],
xlab='Total length (mm)', ylab='Genome size (pg)')

legend("topleft", legend=levels(genome$Suborder), #Add legend at top left corner
       col= myColours, pch = mySymbols)<jupyter_output><empty_output><jupyter_text>It is clear that the two suborders have very different relationships: to begin with we will look at dragonflies (Anisoptera). We will calculate two linear models:* **The Null Model**:  This is the simplest linear model: nothing is going on and the response variable just has variation around the mean: $y = \beta_1$. This is written as an R formula as `y ~ 1`.* **The Linear (Regression) Model**: This models a straight line relationship between the response variable and a continuous explanatory variable: $y= \beta_1 + \beta_{2}x$.The code below fits these two models.<jupyter_code>nullModelDragon <- lm(logBW ~ 1, data = genome, subset = Suborder ==  "Anisoptera")
genomeSizeModelDragon <- lm(logBW ~ logGS, data = genome, subset =  Suborder == "Anisoptera")<jupyter_output><empty_output><jupyter_text>* Note the long names for the models. Short names are easier to type but calling R objects names like `mod1`, `mod2`, ` xxx` swiftly get confusing!*$\star$ Add these model fitting commands into your script and run them.Now we want to look at the output of the model. Remember from the lecture that a model has *coefficients* (the $\beta$ values in the equation of the model) and *terms* which are theexplanatory variables in the model. We'll look at the *coefficients* first:<jupyter_code>summary(genomeSizeModelDragon)<jupyter_output><empty_output><jupyter_text>There is a lot of information there: the model description ('`Call`'), a summary of the residuals, a table of coefficients and then information on residual standard error, r squared and an $F$ test.All of these will become clearer during this course — for the moment, concentrate on the coefficients table.There are two rows in the coefficient table, one for each coefficient in $y=\beta_1 + \beta_2x$ — these are the intercept and the slope of the line. The rest the details on each row are a $t$ test of whether the slope and intercept are significantly different from zero.Now we will look at the *terms* of the model using the `anova` function. We will have a proper look at ANOVA (Analysis of Variance) [later](15-anova.ipynb). Meanwhile, recall from your lecture that ANOVA tests how much variation in the response variable is explained by each explanatory variable. We only have one variable and so there is only one row in the output:<jupyter_code>anova(genomeSizeModelDragon)<jupyter_output><empty_output><jupyter_text>This table is comparing the variation in log body weight explained by log genome size to the total variation in log body weight. We are interested in how much smaller the residuals are for the genome size model than the null model. Graphically, how much shorter are the red residuals than the blue residuals:---                                             Comparison of the null and the regression models.                        ---We can get the sums of the squares of these residuals from the twomodels using the function `resid`, and then square them and add them up:<jupyter_code>sum(resid(nullModelDragon) ^ 2)
sum(resid(genomeSizeModelDragon) ^ 2)<jupyter_output><empty_output><jupyter_text>So we have five columns in the ANOVA table:* **`Df`**: This shows the degrees of freedom. Each fitted parameter/coefficient takes up a degree of freedom from the total sample size, and the left over are the residuals degree of freedom. In this case, genome size adds a slope (compare the null model $y=\beta_1$ and this model $y=\beta_1 + \beta_2x$ — there is one more $\beta$).* **`Sum Sq`**: This shows sums of squares. The bottom line is the residual sum of squares for the model and the one above is the variation explained by genome size. Using the two values from above, the sums of square residuals for the null model are 36.67. In the genome size model, the sum of square residuals are 28.14 and so $36.67-28.14=8.53$ units of variance have been explained by this model.* **`Mean Sq`**: These are just the Sum Sq (Sum of Squares) values divided by the degreesof freedom. The idea behind this is simple: if we explain lots of variation with one coefficient, that is good (the null model), and if we explain a small amount of variation with a loss of degree of freedom (by adding and then estimating more parameters), then that is bad.* **`F value`**: This is the ratio of the Mean Sq for the variable and the residual MeanSq. This is used to test whether the explained variation is large or small.* **`Pr(>F)`**: This is the $p$ value — the probability of the x-variable (the fitted model) explaining this much variance by chance.In this case, it is clear that genome size explains a significant variation in body weight.$\star$ Include the `summary` and `anova` commands for `genomeSizeModelDragon` in your script, run them and check you are happy with the output. Exercise Using the above code as a template, create a new model called `genomeSizeModelDamsel` that fits log body weight as a function of log genome size for damselflies. Write and run code to get the `summary` and `anova` tables for this new model. For example, the first step would be:<jupyter_code>genomeSizeModelDamsel <- lm(logBW ~ logGS, data=genome,subset=Suborder=='Zygoptera')<jupyter_output><empty_output><jupyter_text>Plotting the fitted Regression modelNow we can plot the data and add lines to show the models. For simple regression models, we can use the function `abline(modelName)` to add a line based on the model.$\star$ You already know how to create and customise scatterplots from previous chapters. Create a plot of log body weight as a function of log genome size, picking your favourite colours for the points.Use `abline` to add a line for each model and use the ` col` option in the function to colour each line to match thepoints. For example: `abline(genomeSizeModelDragon, col='red')`.<jupyter_code>myCol <- c('red','blue')

plot(logBW ~ logGS, data=genome, col=myCol[Suborder], xlab='log Genome Size (pg)', ylab='log Body Weight (g)')

abline(genomeSizeModelDragon, col='red')<jupyter_output><empty_output><jupyter_text>Your final figure should look something like this:                                             Linear regression models fitted to the body weight vs. genome size to the Dragonfly (red) and Damselfly (blue) subsets of the data.                         Model diagnosticsNow that we have our models, we need to check that they are appropriate for the data. For this, we will inspect "diagnostic plots". Producing diagnostic plots is easy in R — if you `plot` the model object, then R produces a set of diagnostic plots!$\star$ Try the following code (and include in the R script file):<jupyter_code>par(mfrow = c(2, 2), mar = c(5, 5, 1.5, 1.5))
plot(genomeSizeModelDragon)<jupyter_output><empty_output><jupyter_text>These are the Diagnostics for the `lm` fit to the Dragonfly data subset.<jupyter_code>par(mfrow = c(2, 2), mar = c(5, 5, 1.5, 1.5))
plot(genomeSizeModelDamsel)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07437348423605497, 'fraction_numerical': 0.007089111373670792, 'mean_word_length': 5.009715994020927, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 78, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '592308', 'license': 'MIT', 'n_tokens_mistral': 4469, 'n_tokens_neox': 4132, 'n_words': 2199, 'path': 'notebooks/14-regress.ipynb', 'repo_name': 'mathemage/TheMulQuaBio', 'size': '16081'}