17.6 – ANCOVA – analysis of covariance

Introduction

Analysis of covariance (ANCOVA) is intended to help with analysis of designs with categorical treatment variables on some response (dependent) variable, but a known confounding variable is also present. Thus, the researcher is also likely to know of additional ratio scale variables that covary with the response variable and, moreover, must be included in the experimental design in some way.

Take for example the well-known relationship between body size and whole-animal metabolic rate as measured by rates of carbon dioxide production or rates of oxygen consumption for aerobic organisms. We may be interested in how addition or blocking of stress hormones affects resting metabolism; we may be interested in comparing men and women for activity metabolism, and so on. We’d like to know if the regressions were the same (eg, metabolic rate covaried with body mass in the same way — that is, the slope of the relationship was the same).

This situation arises frequently in biology. For example, we might want to know if male and female birds have different mean field metabolic rates, in which case we might be tempted to use a one-way ANOVA or t-test (since one factor with two levels). However, if males and females also differ for body size, then any differences we might see in metabolic rate could be due to differences in metabolic rate are confounded by differences in the covariable body size. We already discussed one approach to correction: calculate a ratio. Thus, a logical approach to correcting or normalizing for the covariation would be to divide body mass (units of kilograms) into metabolic rate (eg, volume of oxygen, O₂, consumed), and make comparisons, say, among different species, on mass-specific trait (). However, because the regression between mass and metabolic rate is allometric, i.e., not equal to one, the ratio does not, in fact normalize for body mass.  We made this point in Chapter 6.2, and remarked that analysis of covariance ANCOVA was a solution.

ANCOVA allows you to test for mean differences in traits like metabolic rate between two or more groups, but only after first accounting for covariation due to another variable (eg, body size). However, ANCOVA makes the assumption that relationship between the covariable and the response variable is the same in the two groups. This is the same as saying that the regression slopes are the same. We discussed how to use t-test to test hypothesis of equal slopes between regression models in Chapter 17.5, but a more elegant way is to include this in your model.

Example

We return to our sample of 13 tadpoles (Rana pipiens), hatched in the laboratory. I’ve repeated the data set in this page, scroll or click here.

Our linear model was and a scatterplot of the data set is shown in Figure 1 (a repeat of Figure 3 from Chapter 17.5, but now points identified to developmental group).

Figure 1. Scatterplot of oxygen consumption by R. pipiens tadpoles vs body mass (g) by developmental group (Gosner stages I or II).

The project looked at whether metabolism as measured by oxygen consumption was consistent across two developmental stages. Metamorphosis in frogs and other amphibians represents profound reorganization of the organism as the tadpole moves from water to air. Thus, we would predict some cost as evidenced by change in metabolism associated with later stages of development. Figure 2 shows a box plot of tadpole oxygen consumption by Gosner (1960) developmental stage (Figure 2 is a repeat of Figure 4 Chapter 17.5).

Figure 2. Copy of Figure 4, Chapter 17.5; boxplot of oxygen consumption of R. pipiens tadpoles by Gosner developmental stages.

Looking at Figure 2 we see a trend consistent with our prediction; developmental stage may be associated with increased metabolism. However, older tadpoles also tend to be larger, and the plot in Figure 2 does not account for that. Thus, body mass is a confounding variable in this example. There are several options for analysis here (eg, ANCOVA), but one way to view this is to compare the slopes for the two developmental stages. While this test does not compare the means, it does ask a related question: is there evidence of change in rate of oxygen consumption relative to body size between the two developmental stages? The assumption that the slopes are equal is a necessary step for conducting the ANCOVA.

So, divide the data set (Table 1) into two groups by developmental stage (12 tadpoles could be assigned to one of two developmental stages; one was at a lower Gosner stage than the others and so is dropped from the subset.

Table 1. and of 12 R. pipiens tadpole larvae by Gosner developmental stage.

Gosner stage I

1.76	109.41
1.88	329.06
1.95	82.35
2.13	198
2.26	607.7

Gosner stage II

Body mass
2.28	362.71
2.35	556.6
2.62	612.93
2.77	514.02
2.97	961.01
3.14	892.41
3.79	976.97

The slopes and standard errors we obtained in Chapter 17.5 were

Make a plot (Figure 3).

Figure 3. Scatterplot with best-fit regression lines of by for Gosner State I (closed circle, solid line) and Gosner Stage II (open circle, dashed line) R. pipiens tadpoles.

R code for plot in Figure 3.

#Used Rcmdr scatterplot(), then modified code
scatterplot(VO2~Body.mass | Gosner, regLine=FALSE, smooth=FALSE, 
boxplots=FALSE, xlab="Body mass (g)", ylab="Oxygen consumption (ml/h)", 
main="", cex=1.4, cex.axis=1.5, cex.lab=1.5, pch=c(19,19), by.groups=TRUE,
col=c("blue","red"), grid=FALSE, 
legend=list(coords="bottomright"), data=Tadpoles)
#Get regression equations for groups, subset by Gosner
abline(lm(VO2~Body.mass, data=Stage01), lty=1, lwd=2, col="blue")
abline(lm(VO2~Body.mass, data=Stage02), lty=1, lwd=2, col="red")

Returning to the important question, are the two slopes statistically indistinguishable (Ho: b_I = b_II), where b_I is the slope for the Gosner Stage I subset and b_II is the slope for the Gosner Stage II subset? We look at the plot, and since the lines cross, we tend to see a difference. Of course, we need to consider that our perception of slope differences may simply be chance, especially because the sample size is small. Proceed to test.

R code

The ANCOVA is a new ANOVA model where the factor variables are adjusted or corrected for the effects of the continuous variable.

R code for ANCOVA example, crossed or interaction model.

tadpole.1 <- lm(VO2 ~ Body.mass*Gosner, data=example.Tadpole)
summary(tadpole.1)
Anova(tadpole.1, type="II")

Output

summary(tadpole.1)

Call:
lm(formula = VO2 ~ Body.mass * Gosner, data = example.Tadpole)

Residuals:
    Min      1Q   Median      3Q      Max 
-167.80 -117.93    13.81   94.66   214.65

Coefficients:
                       Estimate   Std. Error   t value    Pr(>|t|) 
(Intercept)             -1231.6       783.0     -1.573    0.1544 
Body.mass                 750.0       390.7      1.919    0.0912 .
Gosner[T.II]              790.4       859.2      0.920    0.3845 
Body.mass:Gosner[T.II]   -350.1       409.5     -0.855    0.4174 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 155.8 on 8 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.821, Adjusted R-squared: 0.7539 
F-statistic: 12.23 on 3 and 8 DF, p-value: 0.002336

The coefficients for the first factor (GII) and then the differences in the coefficient for the second factor. You can just add the second coefficient to the first so they’re on the same scale.

Anova Table (Type II tests)

Response: VO2
                  Sum Sq   Df    F value     Pr(>F) 
Body.mass         330046    1    13.6030   0.006146 **
Gosner              5630    1     0.2321   0.642908 
Body.mass:Gosner   17736    1     0.7310   0.417423 
Residuals 194102 8 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Suggests interaction is not significant, i.e., the slopes are not different.

We can then proceed to check to see if the intercepts are different, now that we’ve confirmed no significant difference in slope.

R code for ANCOVA as additive model

tadpole.2 <- lm(VO2 ~ Body.mass + Gosner, data=example.Tadpole)
summary(tadpole.2)
Anova(tadpole.2, type="II")

Output

> summary(tadpole.2)

Call:
lm(formula = VO2 ~ Body.mass + Gosner, data = example.Tadpole)

Residuals:
    Min       1Q    Median      3Q      Max 
-163.12  -125.53    -20.27   83.71   228.56

Coefficients:
                   Estimate   Std. Error     t value    Pr(>|t|) 
(Intercept)         -595.37       239.87      -2.482     0.03487 * 
Body.mass            431.20       115.15       3.745     0.00459 **
Gosner[T.II]          64.96       132.83       0.489     0.63648 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 153.4 on 9 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8047, Adjusted R-squared: 0.7613 
F-statistic: 18.54 on 2 and 9 DF, p-value: 0.0006432

Anova(tadpole.2, type="II")
Anova Table (Type II tests)

Response: VO2
              Sum Sq    Df    F value      Pr(>F) 
Body.mass     330046     1    14.0221   0.004593 **
Gosner          5630     1     0.2392   0.636482 
Residuals     211839     9 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Instead of the interaction, or as added, try a nested model, with body mass nested within stage.

tadpole.3 <- lm(VO2 ~ Body.mass/Gosner, data=example.Tadpole)

summary(tadpole.3)

Call:
lm(formula = VO2 ~ Body.mass/Gosner, data = example.Tadpole)

Residuals:
    Min       1Q    Median      3Q      Max 
-168.66  -131.14    -20.28   90.33   225.36

Coefficients:
                 Estimate      Std. Error    t value Pr(>|t|) 
(Intercept)           -575.10      319.51     -1.800   0.1054 
Body.mass              423.65      162.50      2.607   0.0284 *
Body.mass:Gosner[T.II]  21.95       63.73      0.344   0.7384 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 154.4 on 9 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8021, Adjusted R-squared: 0.7581 
F-statistic: 18.24 on 2 and 9 DF, p-value: 0.0006823

Gets the true coefficient (nested lm() version).

The two test different hypotheses.

lm(VO2 ~ Body.mass * Gosner) tests whether or not the regression has a nonzero slope.

lm(VO2 ~ Body.mass */Gosner) test whether or not the slopes and intercepts from different factors are statistically significant.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. An OLS approach was used for the analysis of tadpole oxygen consumption body mass. Consider the RMA approach — would that be a more appropriate regression model? Explain why or why not.
3. Consider an experiment in which you plan to administer a treatment that has a carry-over effect. For example, Compare and contrast “crossed” and “nested” designs.
4. True or False. The nested design option for the ANCOVA assumes the slopes for the two groups of tadpoles for the regression line of by are equal. Explain your choice.
5. Metabolic rates like oxygen consumption over time are well-known examples of allometric relationships. That is, many biological variables (eg, is related as aM^b, where M is body mass, slope b is scaling exponent), and best evaluated on log-log scale. Repeat the analysis above on log₁₀-transformed and for
   • crossed design (eg, tadpole.1 model)
   • added design (eg, tadpole.2 model)
   • nested design (eg, tadpole.3 model)
6. Create the plot and add the fitted lines from crossed design to the plot.

lgVO2 <- log10(VO2)

lm(log10(VO2) ~ log10(Body.mass ), data=example.Tadpole)

Quiz Chapter 17.6

ANCOVA - analysis of covariance

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Data sets

Oxygen consumption, , of Anuran tadpoles, dataset=example.Tadpole

Gosner refers to Gosner (1960), who developed a criteria for judging metamorphosis staging.

Chapter 17 contents

• Introduction
• Simple Linear Regression
• Relationship between the slope and the correlation
• Estimation of linear regression coefficients
• OLS, RMA, and smoothing functions
• Testing regression coefficients
• ANCOVA – Analysis of covariance
• Regression model fit
• Assumptions and model diagnostics for Simple Linear Regression
• References and suggested readings (Ch17 & 18)

17.4 – OLS, RMA, and smoothing functions

Introduction

OLS or ordinary least squares is the most commonly used estimation procedure for fitting a line to the data. For both simple and multiple regression, OLS works by minimizing the sum of the squared residuals. OLS is appropriate when the linear regression assumptions LINE apply. In addition, further restrictions apply to OLS including that the predictor variables are fixed and without error. OLS is appropriate when the goal of the analysis is to retrieve a predictive model. OLS describes an asymmetric association between the predictor and the response variable: the slope b_X for Y ~ b_XX will generally not be the same as the slope b_Y for X ~ b_YY.

OLS is appropriate for assessing functional relationships (i.e., inference about the coefficients) as long as the assumptions hold. In some literature, OLS is referred to as a Model I regression.

Generalized Least Squares

Generalized linear regression is an estimation procedure related to OLS but can be used either when variances are unequal or multicollinearity is present among the error terms.

Weighted Least Squares

A conceptually straightforward extension of OLS can be made to account for situation where the variances in the error terms are not equal. If the variance of Yi varies for each Xi, then a weighting function based on the reciprocal of the estimated variance may be used.

then, instead of minimizing the squared residuals as in OLS, the regression equation estimates in weighted least squares minimizes the squared residuals summed over the weights.

Weighted least squares is a form of generalized least squares. In order to estimate w_i, however, multiple values of Y for each observed X must be available.

Reduced Major Axis

There are many alternative methods available when OLS may not be justified. These approaches, collectively, may be called Model II regression methods. These methods are invariably invoked in situations in which both Y and X variables have random error associated with them. In other words, the OLS assumption that the predictor variables are measured without error is violated. Among the more common methods is one called Reduced Major Axis or RMA.

Smoothing functions

data set: atmospheric carbon dioxide (CO₂) readings Mauna Loa. Source: https://gml.noaa.gov/ccgg/trends/data.html

Fit curves without applying a known formula. This technique is called smoothing and, while there are several versions, the technique involves taking information from groups of observations — weighted averaged — and using these groups to estimate how the response variable changes with values of the independent variable. Smoothing is used to help reveal patterns, to emphasize trends by reducing noise — clearly, caution needs to be employed as smoothing necessarily hides outlier data, which can themselves be important. Smoothing techniques by name include kernel, loess, and spline. Default in the scatter plot command is loess.

CO₂ in parts per million (ppm) plotted by year from 1958 to 2014 the first CO₂ readings were recorded in April 1958 (Fig 1).

A few words of explanation for Figure 1. The green line shows the OLS line, the red line shows the loess smoothing with a smoothing parameter of 0.5 (in Rcmdr the slider value reads “5”).

Figure 1. CO₂ in parts per million (ppm) plotted by year from 1958 to 2014

R command was started with option settings available in Rcmdr context menu for scatterplot, then additional commands were added

scatterplot(CO2~Year, reg.line=lm, grid=FALSE, smooth=TRUE, spread=FALSE, boxplots=FALSE, span=0.05, lwd=2, xlab="Months since 1958", ylab="CO2 ppm", main="CO2 at Mauna Loa Observatory, April 1958 - April 2014", cex=1, cex.axis=1.2, cex.lab=1.2, pch=c(16), data=StkCO2MLO)

The next plot is for ppm CO₂ by month for the year 2013. The plot shows the annual cycle of atmospheric CO₂ in the northern hemisphere.

Figure 2. Plot of ppm CO₂ by month for the year 2013.

Again, the smoothing parameter was set to 0.5 and the loess function is plotted in red (Fig 2).

Loess is an acronym short for local regression. Loess is a weighted least squares approach which is used to fit linear or quadratic functions of the predictors at the centers of neighborhoods. The radius of each neighborhood is chosen so that the neighborhood contains a specified percentage of the data points. This percentage of data points is referred to as the smoothing parameter and this parameter may differ for different neighborhoods of points. The idea of loess, in fact, any smoothing algorithm, is to reveal pattern within a noisy sequence of observations. The smoothing parameter can be set to different values, between 0 and 1 is typical.

To get a sense of what the parameter does, Figure 3 takes the same data as in Figure 2, but with different values of the smoothing parameter (Fig 3).

Figure 3. Plot with different smoothing values (0.5 to 10.0).

plot key:

The R code used to generate Figure 3 plot was

spanList = c(0.5, 0.75, 1, 2, 10)
reg1 = lm(ppm~Month)
png(filename = "RplotCO2mo.png", width = 400, height = 400, units = "px", pointsize = 12, bg = "white")
plot(Month,ppm, cex=1.2, cex.axis=1.2, cex.lab=1.2, pch=c(16), xlab="Months", ylab="CO2 ppm", main="CO2 levels 2013")
abline(reg1,lwd=2,col="green") 
for (i in 1:length(spanList))
{
ppm.loess <- loess(ppm~Month, span=spanList[i], Dataset)  
ppm.predict <- predict(ppm.loess, Month)
lines(Month,ppm.predict,lwd=2,col=i)
}

The CO2 data constitutes a time series. Instead of loess, a simple moving average would be a more natural way to reveal trends. In principle, take a set of nearby points (odd number of points best, keeps the calculation symmetric) and calculate the average. Next, shift the points by a specified time interval (eg, 7 days), and recalculate the average for the new set of points. See Chapter 20.5 for Time series analysis.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. This is bio-statistics, so I gotta ask: What ecological/environmental process explains the shape of the relationship between ppm CO₂ and months of the year as shown in Figure 2? Hint: NOAA Global Monitoring Laboratory responsible for the CO2 data is located at Mauna Loa Observatory, Hawaiʻi (lat: 19.52291, lon: -155.61586).
3. As I wrote this question (January 2022), we were 22 months since WHO declared Covid-19 a pandemic (CDC timeline). Omicron variant is now dominant; Daily case counts State of Hawaiʻi from 1 November 2021 to 15 January 2022 reported in data set table.
   1. Make plot like Figure 2 (days instead of months)
   2. Apply different loess smoothing parameters and re-plot the data. Observe and describe the change to the trend between case reports and days.

Quiz Chapter 17.4

OLS, RMA, and smoothing functions

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Data set

Covid-19 cases reported State of Hawaiʻi from 1 November 2021 to 15 January 2022 (data extracted from Wikipedia)

Chapter 17 contents

• Introduction
• Simple Linear Regression
• Relationship between the slope and the correlation
• Estimation of linear regression coefficients
• OLS, RMA, and smoothing functions
• Testing regression coefficients
• ANCOVA – analysis of covariance
• Regression model fit
• Assumptions and model diagnostics for Simple Linear Regression
• References and suggested readings (Ch17 & 18)

4.5 – Scatter plots

Introduction

Scatter plots, also called scatter diagrams, scatterplots, or XY plots, display associations between two quantitative, ratio-scaled variables. Each point in the graph is identified by two values: its X value and its Y value. The horizontal axis is used to display the dispersion of the X variable, while the vertical axis displays the dispersion of the Y variable.

The graphs we just looked at with Tufte’s examples Anscombe’s quartet data were scatter plots (Chapter 4 – How to report statistics).

Here’s another example of a scatter plot, data from Francis Galton, data contained in the R package HistData.

Figure 1. Scatterplot of mid-parent (vertical axis) and their adult children’s (horizontal axis) height, in inches. data from Galton’s 1885 paper, “Regression towards mediocrity in hereditary stature.” The red line is the linear regression fitted line, or “trend” line, which is interpreted in this case as the heritability of height.

The commands I used to make the plot in Figure 1 were

library(HistData)
data(GaltonFamilies, package="HistData")
attach(GaltonFamilies)

plot(childHeight~midparentHeight, xlab="Mid-parent height", ylab="Mean adult children's height", xlim=c(55, 80), ylim=c(55,80), cex=0.8, cex.axis=1.2, cex.lab=1.3, pch=c(19), data=GaltonFamilies)
abline(lm(childHeight~midparentHeight), col="red", lwd=2)

I forced the plot function to use the same range of values, set by providing values for xlim and ylim; the default values of the plot command picks a range of data that fits each variable independently. Thus, the default X axis values ranged from 64 to 76 and the Y variable values ranged from 55 to 80. This has the effect of shifting the data, reducing the amount of white space, which a naïve reading of Tufte would suggest is a good idea, but at the expense of allowing the reader to see what would be the main point of the graph: that the children are, on average, shorter than the parents, mean height = 67 vs. 69 inches, respectively. Therefore, Galton’s title begins with the word “regression,” as in the definition of regression as a “return to a former … state” (Oxford Dictionary).

For completeness, cex sets the size of the points (default = 1), and therefore cex.axis and cex.lab apply size changes to the axes and labels, respectively; pch refers to the graph elements or plotting characters, further discussed below (see Fig 8); lm() is a call to the linear model function; col refers to color.

Figure 2 shows the same plot, but without attention to the axis scales, and, more in keeping with Tufte’s principle of maximize data, minimize white space.

Figure 2. Same plot as Figure 1, but with default settings for axis scales.

Take a moment to compare the graphs in Figure 1 and 2. Setting the scales equal allows you to see that the mid-parent heights were less variable, between 65 and 75 inches, than the mean children height, which ranged from 55 to 80 inches.

And another example, Figure 3. This plot is from the ggplot2() function and was generated from within R Commander’s KMggplot2 plug-in.

Figure 3. Finishing times in minutes of 1278 runners by age and gender at the 2013 Jamba Juice Banana 5K in Honolulu, Hawaiʻi. Loess smoothing functions by groups of female (red) and male (blue) runners are plotted along with 95% confidence intervals.

Figure 3 is a busy plot. Because there were so many data points, it is challenging to view any discernible pattern, unlike Figure 1 and 2 plots, which featured less data. Use of the loess smoothing function, a transformation of the data to reduce data “noise” to reveal a continuous function, helps reveal patterns in the data:

1. across most ages, men completed 5K faster than did females and
2. there was an inverse, nonlinear association between runner’s age and time to complete the 5K race.

Take a look at the X-axis. Some runners ages were reported as less than 5 years old (trace the points down to the axis to confirm), and yet many of these youngsters were completing the 5K race in less than 30 minutes. That’s 6-minute mile pace! What might be some explanations for how pre-schoolers could be running so fast?

Design criteria

As in all plotting, maximize information to background. Keep white space minimal and avoid distorting relationships. Some things to consider:

1. keep axes same length
2. do not connect the dots UNLESS you have a continuous function
3. do not draw a trend line UNLESS you are implying causation

Scatter plots in R

We have many options in R to generate scatter plots. We have already demonstrated use of plot() to make scatter plots. Here we introduce how to generate the plot in R Commander.

Rcmdr: Graphs → Scatterplot…

Rcmdr uses the scatterplot function from the car package. In recent versions of R Commander the available options for the scatterplot command are divided into two menu tabs, Data and Options, shown in Figure 4 and Figure 5.

Figure 4. First menu popup in R Commander Scatterplot command, Rcmdr ver. 2.2-3.

Select X and Y variables, choose Plot by groups if multiple grounds are included, eg, male, female, then click Options tab to complete.

Figure 5. Second menu popup in R Commander scatterplot command., Rcmdr ver. 2.2-3

Set graph options including axes labels and size of the points.

Here is another example using the default settings in scatterplot() function in the car package, now the default scatter plot command via R Commander (Fig. 4), along with the same graph, but modified to improve the look and usefulness of the graph (Fig. 6). The data set was Puromycin in the package datasets.

Figure 6. Default scatterplot, package car, from R Commander, version 2.2-4.

Grid lines in graphs should be avoided unless you intend to draw attention to values of particular data points. I prefer to position the figure legend within the frame of the graph, eg, the open are at the bottom right of the graph. Modified graph shown in Figure 7.

Figure 7. Modified scatterplot, same data from Figure 6

R commands used to make the scatter plot in Figure 7 were

scatterplot(rate~conc|state, col=c("blue", "red"), cex=1.5, pch=c(19,19), 
bty="n", reg=FALSE, grid=FALSE, legend.coords="bottomright")

A comment about graph elements in R

In some ways R is too rich in options for making graphs. There are the plot functions in the base package, there’s lattice and ggplot2 which provide many options for graphics, and more. The advice is to start slowly, for example taking advantage of the Figure 8 displays R’s plotting characters and the number you would invoke to retrieve that plotting character.

Figure 8. R plotting characters pch = 1 – 25 along with examples of color.

colors()

  [1] "white"                "aliceblue"            "antiquewhite" 

[655] "yellow3"              "yellow4"              "yellowgreen"  

There’s a lot more to R plotting. For example, you are not limited to just 25 possible characters. R can print any of the ASCII characters 32:127 or from the extended ASCII code 128:255. See Wikipedia to see the listing of ASCII characters.

In addition to pch, you can change the size of the plotting character with cex, color of the character with col for lines or outlines, and bg to change background color. lwd is used to set the character’s border thickness.

Here’s the R code used to generate the graph in Figure 8. Remember, any line beginning with # is a comment line, not an R command.

#create a vector with 26 numbers, from 0 to 25
stuff <-c(0:25)
plot(stuff, pch=c(32:58), cex = 2.5, col = c(1:26), 'xlab' = "pch number", 'ylab' = "Plotting character number")

Is it “scatter plot” or “scatterplot”?

Spelling matters, of course, and yet there are many words for which the correct spelling seems to be like “beauty,” it is “in the eye of the beholder” (Molly Bawn, 1878, by Margaret Hungerford). Scatter plot is one of these — is it one word or two, or is it something else entirely?

Scatter plot is one of these terms: you’ll find it spelled as “scatterplot” or as “scatter plot,” in the dictionary (eg, Oxford English dictionary), with no guidance to choose between them. And I’m not just talking about the differences between British and American English spelling traditions. 

Thus, in these situations as an author, you can turn to which of the spellings is in common use. I first looked at some of the statistics books on my shelves. I selected 14 (bio)statistics textbooks and checked the index and if present, chapters on graphics for term usage.

Table 1. Frequency of use of different terms for scatter plot in 14 (bio)statistics books currently on Mike’s shelves.

Not much help, basically, it is a tie between “scatter plot” and “scatterplot.”

Next, I searched six journals for the interval 1990 – 2016 for use of these terms. Results presented in Table 6 along with journal impact factor for 2014 and number of issues (Table 2).

Table 2. Impact factor and number of issues 1990 – 2016 for six science journals.

My methods? I used the journal’s online search functions for the various usage for scatter plot and the results are shown in Figure 9.

Figure 9. Usage of terms for X Y plots in research articles normalized to number of issues in six journals between 1990 and 2016.

The journals have different numbers of articles; I partially corrected for this by calculating the ratio number of articles with one of the terms divided by the number of issues for the interval 1990 – 2016. It would have been better to count all of the articles, but even I found that to be an excessive effort given the point I’m trying to make here.

Not much help there, although we can see a trend favoring “scatter plot” over any of the other options.

And finally, to completely work over the issue I present results from use of Google’s Ngram Viewer. Ngram Viewer allows you to search words in all of the texts that Google’s folks have scanned into digital form. I searched on the terms in texts between 1950 and 2015, and results are displayed in Figure 10 and Figure 11.

Figure 10. Results from Ngram Viewer for American English, “scatterplot” (blue), “scatter plot” (red), “scatter diagram” (green), “scattergram” (orange), and “XY plot” (purple).

And the same plot, but this time for British sources

Figure 11. Results from Ngram Viewer for British English. See Figure 10 for key.

Conclusion? It looks like “scatterplot” (blue line) is the preferred usage, but it is close. Except for “scattergram” and “XY plot,” which, apparently, are rarely used. After all of this, it looks like you’re free to make your choice between “scatterplot” or “scatter plot.” I will continue to use “scatter plot.”

Bland-Altman plot

Also known as Tukey mean-difference plot, the Bland-Altman plot is used to describe agreement between two ratio scale variables (Bland and Altman 1986, Giavarina 2015), for example agreement between two different methods used to measure the same samples.

Consider use of imageJ by two different observers to record number of pixels of a unit measure (1 cm) on a series of digital images — there’s subjectivity in drawing the lines (where to start, where to end) — do they agree?  Data set below. blandr package, blandr.draw function.

Figure 12. Bland-Altman plot of 1 cm unit measure in pixel number by imageJ from digital images by two independent observers. Purple central region is 95% CI. Lower and upper dashed horizontal lines represent bounds of “acceptable” agreement.

blandr.draw(Obs1, Obs2)

The plot makes it easy to identify questionable points, for example, the one point in upper right quadrant looks suspect.

Visualizing large data sets

Since the early 2000s, size of data sets in biology have increased dramatically, leaving many statistical graphics solutions ill-equipped to tell stories. The data scientists involved with the genomics and health care communities have introduced a number of useful visualizations. We introduce just two, the volcano plot (Fig 13) and the Manhattan-style plot (Fig 14).

Volcano plot

The volcano plot (Fig 13) is a common way to show events that differ between two groups of subjects (eg, p-values), and is common often provides to help visualize gene expression studies of an exposed group vs a control group (eg, fold changes). For the volcano plot, fold changes are reflected on horizontal axis, indicating how much the gene expression level has increased or decreased. The vertical axis typically represents the negative logarithm of the p-value (-logP), which indicates the statistical significance of the change.

Below is R code — with assistance from Gemini AI (Fall 2026) — for a set of five replicate simulated gene expression fold change data sets on stress test (NaCl, H₂O₂) on yeast grown on YPDA plates. Gene expression profiles were simulated for a set of  100 genes (oxidative pathway, DNA repair pathway); as noted, such studies typically include many more genes. 

Figure 13. Volcano plot, gene expression fold change on simulated stress test (NaCl, H₂O₂) yeast data sets, 100 genes (oxidative pathway, DNA repair pathway).

Inspecting the plots (Fig 13), we see a symmetric distribution: The overall distribution of genes is cleanly centered around a Log2 Fold Change (LFC) of zero. This mirrors a realistic transcriptomic profile where the vast majority of the background genome (shown in grey) remains stable, displaying only technical and biological noise. There is a distinct, asymmetrical shift toward positive fold changes among specific gene pathways (targeted activation). Genes belonging to the Oxidative Stress pathway (red) and DNA Repair pathway (blue) consistently surpass both the statistical threshold () and the biological threshold (). Additionally, we see evidence of synergistic effects: The simultaneous application of both stressors (“Both”) drives these pathways to their highest expression levels, pushing them into the top-right quadrant of the plots as highly significant, upregulated features.

R code — with assistance from Gemini AI (Fall 2026) — for yeast gene expression and volcano plots; graphs by base package only.

# ==============================================================================
# 1. SIMULATE 5 REASONABLE DATASETS WITH BIOLOGICAL & TECHNICAL NOISE
# ==============================================================================
set.seed(123) # New seed for realistic distribution

n_genes <- 100
n_datasets <- 5
groups <- c("Control", "Salt", "H2O2", "Both")
n_groups <- length(groups)

# Real-world experiment: each dataset has 3 biological replicates per group
n_reps_per_group <- 3

# Define gene pathways (15 oxidative stress, 15 DNA repair, 70 background/unchanged)
gene_ids <- paste0("Gene_", 1:n_genes)
pathways <- c(rep("Oxidative_Stress", 15), rep("DNA_Repair", 15), rep("Background", 70))

# Initialize a list to hold the 5 realistic datasets
datasets <- list()

for (d in 1:n_datasets) {

# 1. Establish a realistic base expression distribution across 100 genes
# Real data has high-expressing and low-expressing genes
base_expression <- rnorm(n_genes, mean = 7, sd = 2)

# Create a blank sample matrix: 100 genes x 12 samples (4 groups x 3 replicates)
sample_names <- paste0(rep(groups, each = n_reps_per_group), "_R", 1:n_reps_per_group)
expr_matrix <- matrix(0, nrow = n_genes, ncol = length(sample_names))
rownames(expr_matrix) <- gene_ids
colnames(expr_matrix) <- sample_names

# 2. Add realistic baseline biological variation to every sample
for (col_idx in 1:ncol(expr_matrix)) {
expr_matrix[, col_idx] <- base_expression + rnorm(n_genes, mean = 0, sd = 0.4)
}

# 3. Inject subtle, realistic treatment effects (not over-exaggerated)
for (i in 1:n_genes) {
# Every gene gets a tiny bit of random directional drift across groups (real-world noise)
drift <- rnorm(n_groups, mean = 0, sd = 0.15)

if (pathways[i] == "Oxidative_Stress") {
# Subtle upregulation in H2O2 and Both (mean fold change ~ 1.0 to 1.3)
drift[3] <- drift[3] + rnorm(1, mean = 1.1, sd = 0.2) # H2O2
drift[4] <- drift[4] + rnorm(1, mean = 1.3, sd = 0.2) # Both
} else if (pathways[i] == "DNA_Repair") {
# Subtle upregulation in Salt and Both (mean fold change ~ 0.8 to 1.2)
drift[2] <- drift[2] + rnorm(1, mean = 0.9, sd = 0.2) # Salt
drift[4] <- drift[4] + rnorm(1, mean = 1.2, sd = 0.2) # Both
}

# Apply the group drifts to the respective sample columns
for (g_idx in 1:n_groups) {
start_col <- (g_idx - 1) * n_reps_per_group + 1
end_col <- g_idx * n_reps_per_group
expr_matrix[i, start_col:end_col] <- expr_matrix[i, start_col:end_col] + drift[g_idx]
}
}

datasets[[d]] <- expr_matrix
}

# ==============================================================================
# 2. STATISTICAL ANALYSIS (Calculated directly from replicates)
# ==============================================================================
# We will compare the double stressor 'Both' vs 'Control' using real t-tests.
all_lfc <- matrix(0, nrow = n_genes, ncol = n_datasets)
all_pvals <- matrix(0, nrow = n_genes, ncol = n_datasets)
rownames(all_lfc) <- rownames(all_pvals) <- gene_ids

for (d in 1:n_datasets) {
mat <- datasets[[d]]

# Identify the columns for each group
control_cols <- grep("^Control_", colnames(mat))
both_cols <- grep("^Both_", colnames(mat))

for (i in 1:n_genes) {
control_vals <- mat[i, control_cols]
both_vals <- mat[i, both_cols]

# Real Log2 Fold Change (Difference of means since data is on log scale)
all_lfc[i, d] <- mean(both_vals) - mean(control_vals)

# Real Student's t-test
ttest_res <- t.test(both_vals, control_vals, var.equal = TRUE)
all_pvals[i, d] <- ttest_res$p.value
}
}

# ==============================================================================
# 3. GENERATE REALISTIC VOLCANO PLOTS
# ==============================================================================
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 2.5, 1.5))

lfc_threshold <- 1.0 
p_threshold <- 0.05

for (d in 1:n_datasets) {
lfc <- all_lfc[, d]
neg_log10_p <- -log10(all_pvals[, d])

col_vector <- rep("grey75", n_genes)

is_sig_up <- (lfc >= lfc_threshold) & (all_pvals[, d] <= p_threshold)
is_sig_down <- (lfc <= -lfc_threshold) & (all_pvals[, d] <= p_threshold)

col_vector[is_sig_up & pathways == "Oxidative_Stress"] <- "firebrick2"
col_vector[is_sig_up & pathways == "DNA_Repair"] <- "dodgerblue2"
col_vector[(is_sig_up | is_sig_down) & col_vector == "grey75"] <- "grey30"

plot(lfc, neg_log10_p,
col = col_vector,
pch = 19,
cex = 0.9,
xlab = expression(Log ~ Fold ~ Change),
ylab = expression(-Log ~ p - value),
main = paste("Dataset Replicate", d),
xlim = c(-2.5, 2.5), 
ylim = c(0, 5), 
font.main = 1)

abline(v = c(-lfc_threshold, lfc_threshold), col = "grey50", lty = 3)
abline(h = -log10(p_threshold), col = "grey50", lty = 3)

# Set the legend.
if (d == 1) {
legend("topleft", 
legend = c("Oxidative Stress (Up)", "DNA Repair (Up)", "Other Significant", "Not Significant"),
col = c("firebrick2", "dodgerblue2", "grey30", "grey75"), 
pch = 19, 
cex = 0.6,
box.lty = 0,
bg = "transparent")
}
}

# Reset layout parameters
par(mfrow = c(1, 1))

End R code.

About the simulation

1. Symmetric, Balanced Center: Most of the 100 genes (the 70 background genes) are strictly driven by standard biological and technical noise. This distributes them evenly around 0 on the X-axis, meaning you will naturally see negative fold changes alongside positive ones.
2. Realistic Effect Scales: True upregulations are dialed down to max out around to Fold Changes, mirroring standard biological responses on YPDA plates.
3. True T-Tests: Instead of faking p-values, each simulated dataset now contains 3 biological samples per group. The p-values are calculated dynamically using t.test(), generating a mathematically sound distribution pattern on the Y-axis.

Manhattan plot

A Manhattan plot (Fig 14) is a type of scatter plot used to display the statistical significance of large sets of data points simultaneously. Some feature of the experiment, for example location in the genome or sequential test identifiers, is displayed on the horizontal axis. Like the volcano plot (Fig 13), the vertical axis displays the statistical significance (eg, -log₁₀(p-value)).

Figure 14. Manhattan plot, gene expression fold change on simulated stress test (NaCl, H2O2) yeast data sets, 100 genes (oxidative pathway, DNA repair pathway).

Looking over the plots (Fig 14) we see evidence of pathway clustering: By arranging the 100 genes sequentially across three chromosomes, we see distinct localized peaks rather than a random distribution of significant points. We also see chromosome hotspots: Chromosome I (Oxidative Stress genes) and Chromosome II (DNA Repair genes) display pronounced vertical clusters that reliably break through the red significance line across all five experimental replicates. This simulates a real-world genomic architecture where functionally related gene clusters or operon-like structures respond in unison to environmental shifts.

Conversely, Chromosome III (Background genes, ie, baseline Control) remained entirely flat and well below the significance threshold, validating that the observed stress responses are highly specific and not an artifact of global transcriptional runaway.

R code — with assistance from Gemini AI (Fall 2026) —  for Manhattan plot.

# ==============================================================================
# 1. ASSIGN GEOMETRIC CHROMOSOME POSITIONS TO SIMULATED GENES
# ==============================================================================
# We will use the exact same 5 datasets generated in the previous step.
# To plot a Manhattan plot, every gene needs a Chromosome and a Base Pair (BP) position.

# Assign chromosomes: 
# Chromosome 1: Genes 1-40 (Oxidative Stress pathway)
# Chromosome 2: Genes 41-80 (DNA Repair pathway)
# Chromosome 3: Genes 81-100 (Background/Control pathway)
chroms <- c(rep(1, 40), rep(2, 40), rep(3, 20))

# Assign incremental Base Pair positions along each chromosome
bp_positions <- c(seq(1000, 40000, length.out = 40), 
seq(1000, 40000, length.out = 40), 
seq(1000, 20000, length.out = 20))

# To make a continuous X-axis across multiple chromosomes, we calculate cumulative positions
chr_lengths <- c(45000, 45000, 25000)
chr_offsets <- c(0, cumsum(chr_lengths)[1:2])
coord_x <- bp_positions + chr_offsets[chroms]

# Calculate where the chromosome tick marks should sit on the X-axis (the midpoints)
chr_midpoints <- chr_offsets + (chr_lengths / 2)

# ==============================================================================
# 2. GENERATE MANHATTAN PLOTS FOR ALL 5 REPLICATES
# ==============================================================================
# Set up a clean 2x3 plotting grid to display all 5 replicates together
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 2.5, 1.5))

# Standard scientific significance line (-log10(0.05) is ~1.3)
p_threshold <- 0.05
suggestive_line <- -log10(p_threshold)

for (d in 1:n_datasets) {
neg_log10_p <- -log10(all_pvals[, d])

# Alternating color scheme per chromosome to differentiate boundaries
# Chromosome 1 = Slate Blue, Chromosome 2 = Deep Blue, Chromosome 3 = Dark Grey
chr_colors <- c("lightslateblue", "steelblue4", "grey40")
col_vector <- chr_colors[chroms]

# Create the Manhattan plot grid
plot(coord_x, neg_log10_p,
col = col_vector,
pch = 19,
cex = 1.0,
xaxt = "n", # Suppress the standard numeric X-axis to draw custom Chromosome ticks
xlab = "Chromosome",
ylab = expression(-Log[10] ~ p - value),
main = paste("Dataset Replicate", d),
xlim = c(0, sum(chr_lengths)),
ylim = c(0, 5),
font.main = 1)

# Add the standard genome-wide significance horizontal threshold line
abline(h = suggestive_line, col = "firebrick3", lty = 2, lwd = 1.2)

# Draw the custom structural X-axis labeling the specific chromosomes
axis(1, at = chr_midpoints, labels = c("Chr I", "Chr II", "Chr III"))

# Add an explanatory legend to the first plot panel to map pathways to locations
if (d == 1) {
legend("topright", 
legend = c("Chr I (Oxidative Stress)", "Chr II (DNA Repair)", "Chr III (Background)"),
col = c("lightslateblue", "steelblue4", "grey40"), 
pch = 19, 
cex = 0.65,
box.lty = 0,
bg = "transparent")
}
}

# Reset layout parameters back to default single screen layout
par(mfrow = c(1, 1))

About the R code.

The Genomic Peaks: Because these simulated Oxidative Stress genes sit together on Chromosome I and DNA Repair genes sit on Chromosome II,  we see pronounced, beautiful vertical clusters extending upwards over the red significance line.

The Noise Flatline: Chromosome III hosts the background genes, creating a realistically flat landscape where points rarely cross the significance threshold.

Alternating Chromosome Tracks: Using the standard coordinate offset method, the data points stay organized left-to-right without overlapping.

Questions

1. Using our Comet assay data set (Table 1, Chapter 4.2), create scatter plots to show associations between tail length, tail percent, and olive moment. 
2. Explore different settings including size of points, amount of white area, and scale of the axes. Evaluate how these changes change the “story” told by the graph.

Quiz Chapter 4.5

Scatter plots

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Data sets

Number of pixels of 1 cm length from unit measure on ten digital images. Recorded by two independent observers, Obs1 and Obs2.

 Obs1, Obs2
171.026, 171.105
136.528, 138.521
148.084, 144.222
142.014, 140.057
150.213, 153.118
187.011, 195.092
168.760, 168.668
154.302 ,160.381
209.022, 209.876
240.067, 161.805

Gene expression

ccc

Chapter 4 contents

• Graphs and tables (How to report statistics)
• Bar (column) charts
• Histograms
• Box plot
• Mosaic plots
• Scatter plots
• Add a second Y axis
• Q-Q plot
• Ternary plots
• Heat maps
• Graph software
• References and suggested readings