Learning Objectives

Review summary statistic functions

Learn basic plotting to view variables

Use the hist() function to view distributions of variables

Learn to transform variables to meet assumptions of normality

Construct linear models

Open Project

Log in to Google Drive Desktop, open your Environmental Data Analysis project in RStudio, and create a new script named “A3_LinearModels”.

Download Data

R has a number of built-in data packages that developers use for examples. If you ever need to test code or practice in R, you can view a list of available datasets by calling the data() function.

data()

Use the data() function to download the airquality data set. Once loaded, it will appear as an object in your Environment tab, under Data.

data(airquality)

You can use the str() function to inspect the data frame:

str(airquality)

## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

The data set contains air quality measurements in New York from May to September 1973 over a period of 5 months.

Variable	Type	Unit	Description	Detail
Ozone	numeric	parts per billion	mean ozone concentration	Roosevelt island
Solar.R	numeric		Solar radiaion	central park
Wind	numeric	miles per hour	average wind speed	LaGuardia airport
Temp	numeric	Fahrenheit	maximum daily temperature	LaGuardia airport
Month	numeric		Month of observation	5 values
Day	numeric		day of month	153 days total

You will also see that the Month variable is an integer. We will create a new column named Month_Name that expresses Month as an ordered factor.

# Create a new column Month_Name in the airquality data frame that is a factor
airquality$Month_Name <- as.factor(airquality$Month)
# Determine the levels of your new factor column
levels(airquality$Month_Name)

## [1] "5" "6" "7" "8" "9"

# Rename the levels
levels(airquality$Month_Name) <- c("May", "June", "July", "August", "September")

Summary Statistics

Last week, we learned the mean() function, which calculates the mean of a vector, and the length() function, which returns the number of elements in a vector. You can combine these summary functions with subsetting brackets to quickly calculate means, standard deviations, or sample size for various groups in your data.

# Calculate mean of Ozone variable
mean(airquality$Ozone, na.rm = T)

## [1] 42.12931

# Count the number of observations in Ozone variable
length(airquality$Ozone)

## [1] 153

# If you wish to omit rows with missing data
length(na.omit(airquality$Ozone))

## [1] 116

# Calculate the standard deviation of the Ozone variable
sd(airquality$Ozone, na.rm = T)

## [1] 32.98788

# Calculate the mean of the Ozone variable for observations made in May
may <- airquality[airquality$Month_Name == "May",]
mean(may$Ozone, na.rm = T)

## [1] 23.61538

You can also use the summary() function to quickly calculate summary statistics for every numeric variable in your data frame.

summary(airquality)

##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##  NA's   :37       NA's   :7                                       
##      Month            Day           Month_Name
##  Min.   :5.000   Min.   : 1.0   May      :31  
##  1st Qu.:6.000   1st Qu.: 8.0   June     :30  
##  Median :7.000   Median :16.0   July     :31  
##  Mean   :6.993   Mean   :15.8   August   :31  
##  3rd Qu.:8.000   3rd Qu.:23.0   September:30  
##  Max.   :9.000   Max.   :31.0                 
##

This function will return the range, quartiles, median, and mean for each variable. It will also indicate if the column contains any missing values, or NAs.

Challenge 1

Calculate the mean Ozone concentration for the month of August.

Toggle answer

# Create a data frame that contains only August observations
aug <- airquality[airquality$Month_Name == "August",]
# Calculate mean of the Ozone column in the new data frame
mean(aug$Ozone, na.rm = T)

## [1] 59.96154

Summary Plots

Later we will learn to use the package ggplot2, part of the tidyverse, to create attractive graphs in R. For now, it is useful to know how to use the base plotting system in R. It allows you to visually inspect variables quickly and will work for any version of R, even when add-on packages fail to load properly.

Histograms

You can use the hist() function to display the frequency distribution of a variable. Use hist() to show the distribution of the Ozone variable in the airquality data frame.

hist(airquality$Ozone)

All base plots may be customized by including additional arguments. Use ?hist() to open the help file for the hist() function and view potential arguments. For example, you can remove the plot title, relabel the x axis, and change the number of bins by adding the following arguments:

hist(airquality$Ozone, 
     main = NULL,
     xlab = "Tropospheric Ozone (ppb)",
     nclass = 12)

This histogram shows us that Ozone has a right-skewed distribution. If we can transform the data such that they follow a normal distribution, then we can use parametric tests to analyze the data. Let’s try a log10 transformation.

hist(log10(airquality$Ozone))

Challenge 2

Create a histogram for the column of temperature data.
What kind of distribution do the temperature data most closely follow?
What is the mean temperature over the entire study period?
What is the mean temperature for the month of August?
What is the mean temperature for the month of June?

Toggle answer

# Create a histogram for temperature data
hist(airquality$Temp)

# Calculate mean temperature over entire study period
mean(airquality$Temp, na.rm = T)

## [1] 77.88235

# Calculate mean temperature for the month of August
# (Hint: We already created a data frame that contains only August data)
mean(aug$Temp, na.rm = T)

## [1] 83.96774

# Create data frame that contains only data collected in June
jun <- airquality[airquality$Month_Name == "June",]
# Calculate mean of Temp column in jun data frame
mean(jun$Temp, na.rm = T)

## [1] 79.1

Boxplots

You can also compare the distribution of a variable across multiple groups using the function boxplot().
The form of the function is boxplot(DepVar ~ GroupVar, data = dataobject).
Let’s use this function to create a boxplot showing the variation in Ozone concentrations by month:

# Create boxplot
boxplot(Ozone~Month_Name, data = airquality)

Now compare this boxplot to one that plots your log10-transformed Ozone variable.
Note that additional arguments have been added to the boxplot function to rename axes.

# Create boxplot for log10-transformed ozone variable
boxplot(log10(Ozone)~Month_Name,
        data = airquality,
        xlab = NULL,
        ylab = "log10[Tropospheric Ozone (ppb)]")

Challenge 3

Compare the boxplot of your untransformed Ozone variable to your log10-transformed Ozone variable.
Did this transformation improve the distribution of ozone observations within groups? (Hint: Are the distributions more symmetrical?)
Did this transformation improve the homogeneity of variances among groups? (Hint: Do error bars look more similar among groups?)

Comparing Groups using Analysis of Variance (ANOVA)

Once you have established that your data meet the assumptions of normality, or you have transformed your data so that they meet the assumptions of normality, you can analyze your data using parametric tests like an Analysis of Variance (ANOVA), t-test, or Analysis of Covariance (ANCOVA). Here we would like to compare Ozone concentrations among five months, so we will use the simplest ANOVA function, aov().

The form of the function is aov(DepVar ~ IndVar, data = data)
If you wish to add additional variables and their interactions: aov(DepVar ~ IndVar1 * IndVar2, data = data)
If you wish to add additional variables without interactions: aov(DepVar ~ IndVar1 + IndVar2, data = data)

ANOVA model

We can create a new model object and specify the ANOVA formula.
Then we use the summary() function to call the results of the ANOVA:

# ANOVA model
aov1 <- aov(log10(Ozone)~Month_Name, data = airquality)
# Call results of ANOVA
summary(aov1)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Month_Name    4  4.033   1.008   9.163 1.96e-06 ***
## Residuals   111 12.214   0.110                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 37 observations deleted due to missingness

Post-hoc tests

If we want to know which months differ from each other, we can use a post-hoc test such as Tukey’s Honestly Significant Difference (Tukey HSD).

# Run Tukey post-hoc comparisons for ANOVA model
TukeyHSD(aov1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log10(Ozone) ~ Month_Name, data = airquality)
## 
## $Month_Name
##                          diff         lwr         upr     p adj
## June-May          0.184400319 -0.17136565  0.54016629 0.6050267
## July-May          0.465458802  0.21032580  0.72059180 0.0000164
## August-May        0.448743151  0.19361015  0.70387615 0.0000352
## September-May     0.176635722 -0.07181141  0.42508285 0.2867042
## July-June         0.281058483 -0.07470749  0.63682445 0.1908453
## August-June       0.264342831 -0.09142314  0.62010880 0.2447291
## September-June   -0.007764598 -0.35876681  0.34323762 0.9999967
## August-July      -0.016715652 -0.27184865  0.23841735 0.9997526
## September-July   -0.288823081 -0.53727021 -0.04037595 0.0140810
## September-August -0.272107429 -0.52055456 -0.02366030 0.0243361

Challenge 4

Interpret the results of your ANOVA and posthoc tests:
What is the F-value for your ANOVA?
What is the p-value for your ANOVA?
Does ozone concentration differ significantly among months?
If yes, which pairs of months are significantly different?

Linear Regression

Once you have established that your data meet the assumptions of normality, or you have transformed your data so that they meet the assumptions of normality, you can analyze your data using a linear regression model. For example, we know that tropospheric ozone responds to temperature, so I may want to test for a linear relationship between these two variables.

Scatterplots

I typically like to plot my variables first. You can create a scatterplot using base graphics using the plot() function.

# Plot temperature and ozone
plot(x = airquality$Temp, 
     y = airquality$Ozone,
     xlab = "Temperature (F)",
     ylab = "Tropospheric Ozone (ppb)")

Challenge 5

Create a similar scatterplot for your log10-transformed ozone data.
Which lines of the plot code do you need to change to plot your transformed data?
How does the transformation change the shape of the relationship?

Toggle answer

# Plot temperature and log10-transformed ozone
plot(x = airquality$Temp, 
     y = log10(airquality$Ozone),
     xlab = "Temperature (F)",
     ylab = "log10[Tropospheric Ozone (ppb)]")

Running Linear Models

Let’s say we wish to run a linear regression to test for a relationship between temperature and ozone concentration. We can use the lm() function to specify our model.
The form of the function is lm(DepVar~IndVar, data = data)
We can create a new model object and specify the linear model.
Then we use the summary() function to call the results of the linear regression:

# Linear model
mod1 <- lm(Ozone~Temp, data = airquality)
# Call results of linear model
summary(mod1)

## 
## Call:
## lm(formula = Ozone ~ Temp, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.729 -17.409  -0.587  11.306 118.271 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -146.9955    18.2872  -8.038 9.37e-13 ***
## Temp           2.4287     0.2331  10.418  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.71 on 114 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.4877, Adjusted R-squared:  0.4832 
## F-statistic: 108.5 on 1 and 114 DF,  p-value: < 2.2e-16

Including Transformed Variables

But we know that our Ozone variable is not normally distributed. We can run the same model for our log10-transformed variable.

# Linear model for log10-transformed ozone
mod2 <- lm(log10(Ozone)~Temp, data = airquality)
# Call results of linear model
summary(mod2)

## 
## Call:
## lm(formula = log10(Ozone) ~ Temp, data = airquality)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93143 -0.14373  0.01286  0.15855  0.64893 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.798220   0.195866  -4.075 8.53e-05 ***
## Temp         0.029316   0.002497  11.741  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.254 on 114 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.5473, Adjusted R-squared:  0.5434 
## F-statistic: 137.8 on 1 and 114 DF,  p-value: < 2.2e-16

Challenge 6

Interpret the results of the two linear regressions:
What is the R² value for the first regression (including the untransformed ozone variable)?
What is the slope of the first regression (including the untransformed ozone variable)?
What is the R² value for the second regression (including the log10-transformed ozone variable)?
What is the slope of the second regression (including the log10-transformed ozone variable)?
Is there a significant relationship between temperature and ozone concentration?

Adding Trend lines to Scatterplots

Once you have saved your linear model as an object (mod2), you can refer to this object to add a trend line to your scatterplot. You set up your scatterplot as before, and use the function abline() to add a trend line to the base plot.

# Create base plot
plot(x = airquality$Temp, 
     y = airquality$Ozone,
     xlab = "Temperature (F)",
     ylab = "Tropospheric Ozone (ppb)")
# Add trend line, using linear regression mod2
# Argument lwd used to increase width of trend line
abline(mod1, lwd = 2)

Challenge 7

Create a similar scatterplot for your log10-transformed ozone data.
Be sure to add appropriate axis labels, as well as the trendline for the appropriate model.
Consult the “R Base Graphics Cheatsheet” (available on Brightspace) to change one other feature of your graph. Export your final graph by clicking “Export” in the Plots panel, and saving your plot as a .jpg.

Toggle answer

# Create base plot
plot(x = airquality$Temp, 
     y = log10(airquality$Ozone),
     xlab = "Temperature (F)",
     ylab = "log10[Tropospheric Ozone (ppb)]")
# Add trend line, using linear regression mod2
# Argument lwd used to increase width of trend line
abline(mod2, lwd = 2)

Checking Distribution of Model Residuals

Though transforming variables to meet the assumptions of linear models is generally good practice, linear models are robust to most violations of these assumptions. The most important thing is that residual error in your model is normally distributed. That means that your model does equally well explaining variance across the range of observations that you collected. If a model performs better at one end of the relationship than the other, that could present a problem for any predictions you make based on the model. Let’s use the hist() function to check the distribution of residual error in our model of temperature and ozone:

# Plot model residuals
hist(mod1$residuals)

…so that doesn’t look great. This distribution is right-skewed. Let’s inspect the distribution of residual error for our model of temperature and log10-transformed ozone:

# Plot model residuals
hist(mod2$residuals)

This looks better. The model containing our log10-transformed ozone variable results in a more normal distribution of model residuals. When running a linear model, always remember to check that residuals are normally distributed. For small datasets, visual inspection using a histogram works well. For much larger datasets, you may want use statistical techniques like a Lillifors test to determine whether your distribution differs significantly from a normal distribution.

Assignment 3: Data Distributions & Linear Models

Mary Alldred

last updated 2024-01-15

Learning Objectives

Open Project

Download Data

Summary Statistics

Challenge 1

Summary Plots

Histograms

Challenge 2

Boxplots

Challenge 3

Comparing Groups using Analysis of Variance (ANOVA)

ANOVA model

Post-hoc tests

Challenge 4

Linear Regression

Scatterplots

Challenge 5

Running Linear Models

Including Transformed Variables

Challenge 6

Adding Trend lines to Scatterplots

Challenge 7

Checking Distribution of Model Residuals