What is linear regression?

• explains the linear relationship between dependent (Y) and independent (X) variables
• "line of best fit" to data
• a "simple" linear regression is equivalent to a correlation
• can be considered supervised machine learning technique (see here)
  • the model learns from known data, and can be used to predict
  • sometimes we are looking for trends and don't care about predictions

linear regression takes the linear equation of:

\[ Y = mX + b \] But we're going to use it more like:

\[ \hat{y_{i}} = \hat{\beta_{0}} + \hat{\beta_{1}}x_{i} + \epsilon_{i}\]

where:
\(\hat{y_{i}} =\) predicted response for an experimental unit \(i\)
\(x_{i} =\) predictor (or independent variable) of experimental unit \(i\)
\(\hat{\beta_{0}} =\) is expected value when \(x_{i} = 0\) (intercept)
\(\hat{\beta_{1}} =\) slope
\(\epsilon =\) some error, because models are never perfect!

\(\beta_{0}\) and \(\beta_{1}\) are unknown. We will estimate these coefficients from known data. To do this, we need to estimate a line of best fit than minimizes error... this is where linear regression comes in!

How to manipulate data

Sometimes we need to access specific data from the data set because we aren't interested in the entire thing. R let's us manipulate data a number of different ways but we'll go over some of the most important ones.

• data[rows,columns]
• e.g. data[1:20, c(2,5,6)] would give us rows 1-20 and columns 2, 5, 6, of "data"
• data$column_name
• e.g. data$species would give us the column names "species" in "data"

Play around with your penguins data set in here with these functions. I'll give you a head start

### Use this space to complete your exercise.
penguins[1:4,c(2:4)]

penguins$island 

## your turn! :)

We will focus on two penguin measurements for today:

Artwork by @allison_horst.

Simple visual check

plot(penguins$bill_depth_mm, penguins$bill_length_mm,
     main="Bill length vs bill depth", xlab = "Bill depth (mm)", ylab = "
     Bill lengh (mm)")

Hmmm... not sure yet. But! Remember! There are three types of penguins in the data. Maybe we'll see some patterns later on...

Summary output

Our summary output gives us a lot of information:

1. Info about distribution of residuals (errors)
2. The estimates of our coefficients
3. Standard error of coefficient estimates
   • square root of the variance
   • \(\sqrt{\sigma}^{2}\)
4. t-values (coefficient estimates/standard error)
5. \(R^{2}\) = (want closer to 1)
6. p-values
   • testing if your coefficient = 0

Let's refer back to the linear model equation.

\[ \hat{y_{i}} = \hat{\beta_{0}} + \hat{\beta_{1}}x_{i} + \epsilon_{i}\]

Here, \(\hat{y_{i}}\) is bill length, \(x_{i}\) is bill depth. Now we've estimated:

• \(\hat{\beta_{0}}\), the intercept = 55.0674
• \(\hat{\beta_{1}}\) the bill depth coefficient = -0.6498

Note: \(\hat{\beta_{0}}\), the intercept, is the intercept of the linear model when \(x=0\), so we don't see it on this graph.

We can re-write out equation as:

\[ \hat{y_{i}} = 55.0674 - 0.6598x_{i} + \epsilon_{i}\]

Multiple linear model without interations, lm2

lm2 will find different slopes for each species in our data, while maintaining the same intercept. This is a linear model without interactions. The confusing part is when we do these multiple regression models, there is always a "reference" variable. In this case, our reference will always be Adelie penguins. This reference was chosen automatically because it starts with 'A'.

We've estimated more coefficients this time, so our equation is getting a bit longer. Now it looks like:

\[ \hat{y_{i}} = 13.2164 + 1.3948x_{i} + 9.939\text{Chinstrap} + 13.4033 \text{Gentoo} + \epsilon_{i}\]

1. Notice the positive slope! This is different from last time.
2. What are we supposed to do with "Chinstrap" and "Gentoo"?

"Chinstrap" and "Gentoo" are in there to make our equation "easier" on us. If the data we are analysing are from a Chinstrap penguin, then Chinstrap will become a "1" and Gentoo a "0" and vice versa. If the data are from an Adelie penguin, both will be 0.

This is what the Adelie equation looks like:

\[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + 9.939(0) + 13.4033(0) + \epsilon_{i}\] \[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + \epsilon_{i} \] Chinstrap:

\[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + 9.939(1) + 13.4033(0) + \epsilon_{i}\] \[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + 9.939 \epsilon_{i}\] \[ \hat{y_{i}} = 23.1554 + 1.3940x_{i} +\epsilon_{i}\]

Gentoo:

\[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + 9.939(0) + 13.4033(1) + \epsilon_{i}\] \[ \hat{y_{i}} = 13.2164 + 1.3940x_{i} + 13.4033(1) + \epsilon_{i}\] \[ \hat{y_{i}} = 26.6197 + 1.3940x_{i} + \epsilon_{i}\]

What if we used a different penguin as the reference? Let's try it out

penguins$species <- relevel(as.factor(penguins$species), ref="Chinstrap")
lm2_relevel <- lm(formula = bill_length_mm ~ bill_depth_mm + species, data = penguins)
summary(lm2_relevel)

Notice how the slope doesn't change? So far, we've only told R to look for different slopes for each species.

Multiple linear model with interations, lm4

In this regression, lm4 is giving us a model that includes: - different intercepts per species - different bill length coefficients per species

Again, Adelie is the reference penguin, so everything is compared to Adelie. But this time we have MORE coefficients since we're looking at how the slope is affected.

Full equation

\[ \hat{y_{i}} = 23.0681 + 0.857x_{i} + (-9.6402\text{Chinstrap}) + (-5.8386 \text{Gentoo}) + 1.0651x_{i}\text{Chinstrap} + 1.1637x_{i}\text{Gentoo} + \epsilon_{i}\]

Adelie's equation: \[ \hat{y_{i}} = 23.0681 + 0.857x_{i} -9.6402(0) -5.8386(0) + 1.0651x_{i}(0) + 1.1637x_{i}(0) + \epsilon_{i}\] \[ \hat{y_{i}} = 23.0681 + 0.857x_{i} + \epsilon_{i}\] Chinstrap's equation: \[ \hat{y_{i}} = 23.0681 + 0.857x_{i} -9.6402(1)-5.8386(0) + 1.0651x_{i}(1) + 1.1637x_{i}(0) + \epsilon_{i}\] \[ \hat{y_{i}} = 23.0681 + 0.857x_{i} -9.6402 + 1.0651x_{i} + \epsilon_{i}\] \[ \hat{y_{i}} = 13.4279 + 0.857x_{i} + 1.0651x_{i} + \epsilon_{i}\] \[ \hat{y_{i}} = 13.4279 + 1.9221x_{i} + \epsilon_{i}\]

What would Gentoo's equation look like?

Finding the mean bodyweight of penguins

We have a pretty big dataset of info on penguins. We don't really need all of it, do we? Right now we're only interested in body weight info of Chinstrap and Gentoo penguins.

We're going to have to subset our data so that we can focus on what we're interested in right now. We are going to be using the function filter(). To do this, we're also going to have to learn about "logical operators".

You will use logical operators often in R programming. Here are some examples:

• ==: equals
• !=: does not equal
• &: AND
• |: OR
• >: greater than
• >: less than
• >=: greater or equal to
• <=: less than or equal to

hote: We're using functions filter() and select(). These functions are from the tidyverse and use special "syntax". Most tidyverse functions always consider the first argument as the data, so we don't need explicitly write data = penguins. Fun!

Today we want to filter for Chinstrap OR (|) Gentoo penguins. It might seem like we're interested in Chinstrap AND (&) Gentoo penguins, but in programming logic, that is impossible. Let's try both and see what happens..

## loading ALL tidyverse packages at once because we'll be using many of them, normally required
#library(tidyverse)
## filtering for Chinstrap AND Gentoo penguins from the species column of penguins
penguin_chin_gent <- filter(penguins, species == "Chinstrap" & species == "Gentoo")
head(penguin_chin_gent)

## filtering for Chinstrap OR Gentoo penguins from the species column of penguins
penguin_chin_gent <- filter(penguins, species == "Chinstrap" | species == "Gentoo")
# print the first couple rows
head(penguin_chin_gent)

Now, we're only interested in the bodyweight of the penguins. We should make an even smaller dataset that only includes the body weight information of Gentoo and Chinstrap penguins.

penguin_chin_gent <- filter(penguins, species == "Chinstrap" | species == "Gentoo")

## what are the column names of our dataset again? Let's check
colnames(penguin_chin_gent)

## OK, so body weight is callled "body_mass_g". Let's select body_mass_g and species so that we
## can calculate the mean body weight of each species.
penguin_chin_gent_bw <- select(penguin_chin_gent, species, body_mass_g)
head(penguin_chin_gent_bw)
tail(penguin_chin_gent_bw)

## what does this function do?
count(penguin_chin_gent_bw, species)

Now that we have our dataset, how about we plot some data???

penguin_chin_gent <- filter(penguins, species == "Chinstrap" | species == "Gentoo")
penguin_chin_gent_bw <- select(penguin_chin_gent, species, body_mass_g)

Plotting bill length vs depth with ggplot

Let's use a scatter plot to explore our data. Modify this code by updating the ???? to make a scatter plot of penguin bill length by bill depth that has points coloured in by the penguin species.

ggplot(data=penguins,
           y=????, 
           colour=????)
       ) + 
  geom_????()

ggplot(data=penguins, 
       aes(x=bill_length_mm, 
           y=bill_depth_mm, 
           colour=species)
       ) + 
  geom_point()

Adding our linear model to our plot

Looking at our plot above, there may be some patterns that can be explained by species. We also found this when we completed our linear models. Cool! So, let's add our linear models to our scatter plot. And while we're at it, let's add custom colours to our plots.

Let's go back to our ggplot2 cheatsheet. Under the Two Variables section, check out geom_smooth(). This geom can calculate a linear model and plot it for us, kind of like the abline() we used previously. In this case, we can plot it all at once. Let's try it out!

ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species)
       ) +
  geom_point() +
  geom_smooth(method = "lm") # method = "lm" tells ggplot we want to use the linear model function lm()

Looks good so far, except for a few things.

1. geom_smooth adds an error estimate from the linear model by default (those grey ribbons around the lines). It's not really helpful for us here, so we should remove that ribbon for now.
2. We should be more creative with the colours we pick.
3. Just in case anyone looking at our graph has a hard time telling colours apart, let's also change the points in our plot to shapes that match species.
4. We should update our axes labels to be more informative.
5. We keep getting these warnings that tell us we have NAs and that ggplot is removing them... we probably shouldn't ignore these warnings and should investigate.

Let's attack these points one at a time.

1. Removing error estimate from geom_smooth

To remove these unnecessary ribbons, we need to tell geom_smooth() that we don't want them.

(no_ribbon <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species)
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE)) #se = FALSE tells ggplot

2. Updating colours

One thing I love about making plots is making them look fun and pretty. Problem is...I'm not very artistically inclined. I use a lot of different websites for plot color theme inspiration. One of my favourites is coolors.co. It randomly generates complementary colour palettes and I use this all the time. Coolors gives us colours in HEX codes that we can use, or we can type in colours manually (here's a link to all the colours R understands)

Let's make our plot pretty both ways: using HEX codes and the names of colours. We're going to use scale_colour_manual() to tell ggplot that we want to manually change the colours.

(new_colours <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species)
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("royalblue", "orange", "limegreen")))

I like the blue and orange...but I'm not so sure that the green goes. I'm going to use coolors.co to help me out!

Colour palette by coolors.

I like this palette that coolors suggested. That combination of letters and numbers is the colour's HEX code! We have so many more colour options using these codes. We can input these codes instead of the colour name using the exact same ggplot function. Feel free to edit these HEX codes as you'd like!

(new_colours <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species)
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD"))) # hex codes start with a "#"

3. Adding shapes

It's useful to not only change the colours of the points, but also have shapes reflect the data as well. This makes our figures more accessible! We have to add an extra option to aes(). In this case, we add shape = species.

(new_colours <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species,
       shape = sex, 
       group = species) 
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE)) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD")) # hex codes start with a "#"

However, what if we're interested in adding more information to our plot? What about if we wanted to look at how the penguin's sex is reflected in these bill measurements? We would tell ggplot to change either the colour or the shape from species to sex. Try it out in the above plot. What happens?

Yikes.I bet something didn't look so right when we added that extra grouping. That's because our original ggplot(data = penguins,aes(x = bill_length_mm, y = bill_depth_mm,colour = species,shape = sex) call told ggplot to use of that information for the rest of the plot. Our linear model is now trying to add information for both sex and species, making it a little fuzzy. We can force ggplot to group any statistical functions by species while allowing us to colour/shape the rest of the data how we want. To do this, we add group = species to our aes() function. (This may be useful for your mini-homework)

(new_colours <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = sex,
       shape = species, 
       group = species) 
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE)) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD")) # hex codes start with a "#"

Now we can see how ignoring those warnings is coming back to haunt us... What do you think is happening here? Let's discuss

4. Changing axes labels

We're already learned how to do this hint xlab() might ring a bell.... I'll leave it up to you to rename our axes to appropriate names.

(axes <- ggplot(data = penguins,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species,
       shape = sex, 
       group = species) 
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD"))) # hex codes start with a "#"

5. Dealing with warnings

Back to those warnings... ggplot is telling us that our dataset is incomplete and contains NA. Researchers use NA when recording data if they were unable to take the measurement, or the measurement is missing for some reason. Keeping these NA is not always helpful to us. Thankfully we have a function that can help us remove any rows (observations) with NA. This function is drop_na().

penguins_removeNA <- drop_na(penguins)
(axes <- ggplot(data = penguins_removeNA,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species,
       shape = sex, 
       group = species) 
       ) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD"))) # hex codes start with a "#"

Yay! We've removed our warnings!

6. Our final plot!

Here's an example of a final plot that went a bit further than what we did in our steps above.

What has been added to this plot? If you change some options around, what happens? You might want to use some of these options for your "homework".

(final <- ggplot(data = penguins_removeNA,
       aes(x = bill_length_mm,
       y = bill_depth_mm,
       colour = species,
       shape = species, 
       group = species)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se=FALSE) + #se = FALSE tells ggplot
  scale_colour_manual(values=c("#1E3888", "#FFAD69", "#47A8BD")) + # hex codes start with a "#"
 labs(x = "Bill length (mm)",
       y = "Bill depth (mm)",
       color = "Penguin species",
       shape = "Penguin species")+
   theme_bw())