18.1 – Multiple Linear Regression

Introduction

Last time we introduced simple linear regression:

• one independent X variable
• one dependent Y variable.

The linear relationship between Y and X was estimated by the method of Ordinary Least Squares (OLS). OLS minimizes the sum of squared distances between the observed responses, and responses predicted by the line. Simple linear regression is analogous to our one-way ANOVA — one outcome or response variable and one factor or predictor variable (Chapter 12.2).

But the world is complicated and so, our one-way ANOVA was extended to the more general case of two or more predictor (factor) variables (Chapter 14). As you might have guessed by now, we can extend simple regression to include more than one predictor variable. In fact, combining ANOVA and regression gives you the general linear model! And, you should not be surprised that statistics has extended this logic to include not only multiple predictor variables, but also multiple response variables. Multiple response variables falls into a category of statistics called multivariate statistics.

Like multi-way ANOVA, multiple regression is the extension of simple linear regression from one independent predictor variable to include two or more predictors. The benefit of this extension is obvious — our models gain realism. All else being equal, the more predictors, the better the model will be at describing and/or predicting the response. Things are not all equal, of course, and we’ll consider two complications of this basic premise, that more predictors are best; in some cases they are not.

However, before discussing the exceptions or even the complications of a multiple linear regression model, we begin by obtaining estimates of the full model, then introduce aspects of how to evaluate the model. We also introduce  and whether a reduced model may be the preferred model.

R code

Multiple regression is easy to do in Rcmdr — recall that we used the general linear model function, lm(), to analyze one-way ANOVA and simple linear regression. In R Commander, we access lm() by

Rcmdr: Statistics → Fit model → Linear model

You may, however, access linear regression through R Commander

We use the same general linear model function for cases of multi-way ANOVA and for multiple regression problems. Simply enter more than one ratio-scale predictor variable and boom!

You now have yourself a multiple regression. You would then proceed to generate the ANOVA table for hypothesis testing

Rcmdr: Models → Hypothesis testing → ANOVA tables

From the output of the regression command, estimates of the coefficients along with standard errors for the estimate and results of t-tests for each coefficient against the respective null hypotheses for each coefficient are also provided. In our discussion of simple linear regression we introduced the components: the intercept, the slope, as well as the concept of model fit, as evidenced by R², the coefficient of determination. These components exist for the multiple regression problem, too, but now we call the slopes partial regression slopes because there are more than one.

Our full multiple regression model becomes

where the coefficients β₁, β₂, … β_n are the partial regression slopes and β₀ is the Y-intercept for a model with 1 – n predictor variables. Each coefficient has a null hypothesis, each has a standard error, and therefore, each coefficient can be tested by t-test.

Now, regression, like ANOVA, is an enormous subject and we cannot do it justice in the few days we will devote to it. We can, however, walk you through a fairly typical example. I’ve posted a small data set diabetesCholStatin at the end of this page. Scroll down or click here. View the data set and complete your basic data exploration routine: make scatterplots and box plots.  We think (predict) that body size and drug dose cause variation in serum cholesterol levels in adult men. But do both predict cholesterol levels?

Selecting the best model

We have two predictor variables, and we can start to see whether none, one, or both of the predictors contribute to differences in cholesterol levels. In this case, both contribute significantly. The power of multiple regression approaches is that it provides a simultaneous test of a model which may have many explanatory variables deemed appropriate to describe a particular response. More generally, it is sometimes advisable to think more philosophically about how to select a best model.

In model selection, some would invoke Occam’s razor — given a set of explanations, the simplest should be selected — to justify seeking simpler models. There are a number of approaches (forward selection, backward selection, or stepwise selection), and the whole effort of deciding among competing models is complicated with a number of different assumptions, strengths and weaknesses. I refer you to the discussion below, which of course is just a very brief introduction to a very large subject in (bio)statistics!

Let’s get the full regression model

The statistical model is

As written in R format, our model is ChLDL ~ BMI + Dose + Statin.

We can now proceed in R Commander to fit the model.

Rmdr: Statistics → Fit models → Linear model

How the model is inputted into linear model menu is shown in Figure 1.

Figure 1. Screenshot of Rcmdr linear model menu with our model elements in place.

The output

summary(LinearModel.1)
Call:
lm(formula = ChLDL ~ BMI + Dose + Statin, data = cholStatins)

Residuals:
    Min      1Q  Median     3Q    Max
-3.7756 -0.5147 -0.0449 0.5038 4.3821

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)           1.016715   1.178430   0.863 0.39041
BMI                   0.058078   0.047012   1.235 0.21970
Dose                 -0.014197   0.004829  -2.940 0.00411 **
Statin[Statin2]       0.514526   0.262127   1.963 0.05255 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.31 on 96 degrees of freedom
Multiple R-squared: 0.1231, Adjusted R-squared: 0.09565
F-statistic: 4.49 on 3 and 96 DF, p-value: 0.005407

Question. What are the estimates of the model coefficients (rounded)?

• = intercept = 1.017
• = slope for variable BMI = 0.058
• = slope for variable Dose = -0.014
• = slope for variable Statin = -0.515

Question. Which of the three coefficients were statistically different from their null hypothesis (p-value < 5%)?

Answer: Only coefficient was judged statistically significant at the Type I error level of 5% (p = 0.0041).  Of the four null hypotheses we have for the coefficients:

• = 0
• = 0
• = 0; we only reject the null hypothesis for Dose coefficient.
• = 0)

The important take-home concept is about the lack of a direct relationship between the magnitude of the estimate of the coefficient and the likelihood that it will be statistically significant! In absolute value terms , but was not close to statistical significance (p = 0.220).

We generate a 2D scatterplot and include the regression lines (by group=Statin) to convey the relationship between at least one of the predictors (Fig 2).

Figure 2. Scatter plot of predicted LDL against dose of a statin drug. Regression lines represent the different statin drugs (Statin1, Statin2).

Question. Based on the graph, can you explain why there will be no statistical differences between levels of the statin drug type, Statin1 (shown open circles) vs. Statin2 (shown closed red circles).

Because we have two predictors (BMI and Statin Dose), you may also elect to use a 3D-scatterplot. Here’s one possible result (Fig 3).

Figure 3. 3D plot of BMI and dose of Statin drugs on change in LDL levels (green Statin2, blue Statin1).

R code for Figure 3.

Graph made in Rcmdr: Graphs → 3D Graph → 3D scatterplot …

scatter3d(ChLDL~BMI+Dose|Statin, data=diabetesCholStatin, fit="linear", 
residuals=TRUE, parallel=FALSE, bg="white", axis.scales=TRUE, grid=TRUE, 
ellipsoid=FALSE)

Figure 4. An example of a possible interactive 3D plot; the file embedded in this page is not interactive, just an animation.

Diagnostic plots

While visualization concerns are important, let’s return to the statistics. All evaluations of regression equations should involve an inspection of the residuals. Inspection of the residuals allows you to decide if the regression fits the data; if the fit is adequate, you then proceed to evaluate the statistical significance of the coefficients.

The default diagnostic plots (Fig 5) R provides are available from

Rcmdr: Models → Graphs → Basic diagnostics plots

Four plots are returned

Figure 5. R’s default regression diagnostic plots.

Each of these diagnostic plots in Figure 4 gives you clues about the model fit.

1. Plot of residuals vs. fitted helps you identify patterns in the residuals
2. Normal Q-Q plot helps you to see if the residuals are approximately normally distributed
3. Scale-location plot provides a view of the spread of the residuals
4. The residuals vs. leverage plot allows you to identify influential data points.

We introduced these plots in Chapter 17.8 when we discussed fit of simple linear model to data. My conclusion? No obvious trend in residuals, so linear regression is a fit to the data; data not normally distributed Q-Q plot shows S-shape.

Interpreting the diagnostic plots for this problem

The “Normal Q-Q” plot allows us to view our residuals against a normal distribution (the dotted line). Our residuals do no show an ideal distribution: low for the first quartile, about on the line for intermediate values, then high for the 3rd and 4th quartile residuals. If the data were bivariate normal we would see the data fall along a straight line. The “S-shape” suggests log-transformation of the response and or one or more of the predictor variables.

There also seems to be a pattern in residuals vs the predicted (fitted) values. There is a trend of increasing residuals as cholesterol levels increase, which is particularly evident in the “scale-location” plot. Residuals tended to be positive at low and high doses, but negative at intermediate doses. This suggests that the relationship between predictors and cholesterol levels may not be linear, and it demonstrates what statisticians refer to as a monotonic spread of residuals.

The last diagnostic plot looks for individual points that influence, change, or “leverage” the regression — in other words, if a point is removed, does the general pattern change? If so, then the point had “leverage” and thus we need to decide whether or not to include the datum. diagnostic plots Cook’s distance is a measure of the influence of a point in regression. Points with large Cook’s distance values warrant additional checking.

The multicollinearity problem

Statistical model building is a balancing act by the statistician. While simpler models may be easier to interpret and, perhaps, to use, it is a basic truism that the more predictor variables the model includes, the more realistic the statistical model. However, each additional parameter that is added to the statistical model must be independent of all other parameters already in the model. To the extent that this assumption is violated, the problem is termed multicollinearity. If predictor variables are highly correlated, then they are essentially just linear combinations and do not provide independent evidence. For example, one would naturally not include two core body temperature variables in a statistical model on basal metabolic rate, one in degrees Fahrenheit and the other in degrees Celsius, because  it is a simple linear conversion between the two units. This would be an example of structural collinearity: the collinearity is because of misspecification of the model variables. In contrast, collinearity among predictor variables may because the data are themselves correlated. For example, if multiple measures of body size are included: weight, height, length of arm, etc., then we would expect these to be correlated, ie, data multicollinearity.

Collinearity in statistical models may have a number of undesirable effects on a multiple regression model. These include

• estimates of coefficients not stable: with collinearity, values of coefficients depend on other variables in the model; if collinear predictors, then the assumption of independent predictor variables is violated.
• precision of the estimates decreases (standard error of estimates increase).
• statistical power decreases.
• p-values for individual coefficients not trust worthy.

Tolerance and Variance Inflation Factor

Absence of multicollinearity is important assumption of multiple regression. A partial test is to calculate product moment correlations among predictor variables. For example, when we calculate the correlation between BMI and Dose for our model, we get r = 0.101 (P = 0.3186), and therefore would tentatively conclude that there was little correlation between our predictor variables.

A number of diagnostic statistics have been developed to test for multicollinearity. Tolerance for a particular independent variable (Xi) is defined as 1 minus the proportion of variance it shares with the other independent variables in the regression analysis (1 − R²_i) (O’Brien 2007). Tolerance reports the proportion of total variance explained by adding the X_i–^th predictor variable that is unrelated to the other variables in the model. A small value for tolerance indicates multicollinearity — and that the predictor variable is nearly a perfect combination (linear) of the variables already in the model and therefore should be omitted from the model. Because tolerance is defined in relation to the coefficient of determination you can interpret a tolerance score as the unique variance accounted for by a predictor variable.

A second, related diagnostic of multicollinearity is called the Variance Inflation Factor, VIF. VIF is the inverse of tolerance.

VIF shows how much of the variance of a regression coefficient is increased because of collinearity with the other predictor variables in the model. VIF is easy to interpret: a tolerance of 0.01 has a VIF of 100; a tolerance of 0.1 has a VIF of 10; a tolerance of 0.5 has a VIF of 2, and so on. Thus, small values of tolerance and large values of VIF are taken as evidence of multicollinearity.

Rcmdr: Models → Numerical diagnostics → Variation-inflation factors

vif(RegModel.2)
      BMI     Dose 
1.010256 1.010256

A rule of thumb is that if VIF is greater than 5 then there is multicollinearity; with VIF values close to one we would conclude, like our results from the partial correlation estimate above, that there is little evidence for a problem of collinearity between the two predictor variables. They can therefore remain in the model.

Solutions for multicollinearity

If there is substantial multicollinearity then you cannot simply trust the estimates of the coefficients. Assuming that there hasn’t been some kind of coding error on your part, then you may need to find a solution. One solution is to drop one of the predictor variables and redo the regression model. Another option is to run what is called a Principal Components Regression (PCR).  One takes the predictor variables and runs a Principal Component Analysis to reduce the number of variables, then the regression is run on the PCA components. By definition, the PCA components are independent of each other. Another option is to use ridge regression approach.

Like any diagnostic rule, however, one should not blindly apply a rule of thumb. A VIF of 10 or more may indicate multicollinearity, but it does not necessarily lead to the conclusion that the linear regression model requires that the researcher reduce the number of predictor variables or analyze the problem using a different statistical method to address multicollinearity as the sole criteria of a poor statistical model. Rather, the researcher needs to address all of the other issues about model and parameter estimate stability, including sample size. Unless the collinearity is extreme (like a correlation of 1.0 between predictor variables!), larger sample sizes alone will work in favor of better model stability (by lowering the sample error) (O’Brien 2007).

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. Please explain why the magnitude of the slope is not the key to statistical significance of a slope? Hint: look at the equation of the t-test for statistical significance of the slope.
3. Consider the following scenario. A researcher measures his subjects multiple times for blood pressure over several weeks, then plots all of the values over time. In all, the data set consists of thousands of readings. He then proceeds to develop a model to explain blood pressure changes over time. What kind of collinearity is present in his data set? Explain your choice.
4. We noted that Dose could be viewed as categorical variable. Convert Dose to factor variable (fDose) and redo the linear model. Compare the summary output and discuss the additional coefficients.
   • Use Rcmdr: Data → Manage variables in active data set →  Convert numeric Variables to Factors to create a new factor variable fDose. It’s ok to use the numbers as factor levels.
5. We flagged the change in LDL as likely to be not normally distributed. Create a log10-transformed variable for ChLDL and perform the multiple regression again.
   1. Write the new statistical model
   2. Obtain the regression coefficients — are they statistically significant?
   3. Run basic diagnostic plots and evaluate for

Quiz Chapter 18.1

Multiple Linear Regression

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

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

Your Answers

Data set

Chapter 18 contents

• Introduction
• Multiple linear regression
• Nonlinear regression
• Logistic regression
• Generalized Linear Squares
• Selecting the best model
• Compare two linear models
• References and suggested readings (Ch17 & 18)

20.8 – Diversity indexes

Introduction

“Biodiversity,” coined in the mid 1980s (see discussion in Sarkar 2021), is an important and a well used term — over 100 million hits in Google search as of Fall 2025. From ecology we learn that all things being equal, high numbers of species living in the same area rather than low or few kinds living in the area is associated with high productivity and is likely a better predictor of ecosystem persistence over time. And in many ways, we generally know the difference between high and low diversity when we see it (college campus vs arboretum vs forest reserve on Oahu, Fig 1).

A

B

C

Figure 1. (A) Chaminade University, a portion of lower campus; (B) A portion of Lyon Arboretum; (C) Portion of Roundtop Forest Reserve. Google satellite images, approximately the same altitude and sized areas.

We all agree biodiversity captures the variety of life on Earth, but quantifying — and asking questions about biodiversity differences among ecosystems — requires quantification. Over the years a variety of metrics have been proposed and considerable knowledge about their form and applicability is required — I point interested readers to the many excellent resources on this topic, including Jost 2006, Tuomisto 2010, and Sîrbu et al 2025, to name just a few. A thorough and free resource about R and diversity indexes can be found at Diversity Analysis, part of Dr David Zelený’s (National Taiwan University, Taipei, Taiwan), excellent site on Analysis of Community Ecology Data in R.

We introduced index concepts in Chapter 6.2 – Ratios and probabilities. At the time, we did not emphasize an important distinction between parameters and indexes. As we introduced before, parameters are considered fixed, but unknown quantities about a population. Diversity indexes are not properties of a community (statistical population), they are estimates from samples. Therefore, there can be no theoretical distribution to which a diversity index can be compared, hence, not a parameter. This has consequences, which we briefly introduce below.

A diversity index is a measure of how many different kinds (eg, species) are present in a dataset.  These indexes are more than a count of the different types; they also account for how common (or rare) a kinds is. Diversity indexes are examples of applied multivariate statistics: two or more predictor variables and two or more response variables.

Richness, R, quantifies how many types or species. Thus, at any specific site one can report number of species found

The general equation of diversity is

There are many varieties of diversity indexes (Bollarapu et al 2024), but three are well known and two are commonly used.

Shannon’s diversity index

H, Shannon’s Index, Shannon Entropy, accounts for abundance and evenness of all species present in an area. Evenness refers to how close the numbers of each species are in an area.

(1)  

No upper limit, but useful for comparing diversity among sites.

Google Scholar, 83,800 hits (October 2025).

Simpson diversity index

Simpson index (ecology) is the number of species or species richness (S) and evenness — proportional abundance — of populations; independently developed, called Herfindahl–Hirschman index (HHI) in economics.

(2)  

Simpson usually expressed as , ranges from 0 to 1.

Google Scholar, 74,900 hits (October 2025).

Other indices of note

Pielou’s evenness index

(3)  

where is the total number of species and is the Shannon’s diversity index.

Google Scholar, 23,600 hits (October 2025).

Fisher’s diversity index, or alpha (α). The formula is

(4)  

where is the total number of species, is the total number of individuals, and is the Fisher’s alpha index.

Google Scholar, 178,000 hits (October 2025).

To find you need to solve the equation for using the known values of and . This requires an iterative approach, such as the Newton-Raphson method, and you start iterations with a reasonable guess for starting value.

Rao’s diversity, or more formally, Rao’s Quadratic Entropy index (Ricotta & Szeidl 2006). In addition to the proportions (relative abundance) of each species, the index adds a measure of distance between species.

(5)  

where is the distance between species i and species j, is relative abundance of species i, with the relative abundance of species j.

The measure of pairwise distances among the species in the community, which can provide insight into functional diversity.

Google Scholar, 17,400 hits (October 2025).

Comparing index properties

Considering each index, we ask: emphasize richness? Relative abundance or evenness? Both Shannon and Simpson address richness.

R packages and other software

Download vegan package. The function diversity() returns Shannon, Simpson, and Fisher diversity indices and species richness.

Download and install the BiodiversityR package, a GUI for biodiversity, suitability and community ecology analysis. Package provides GUI to many vegan functions, utilizes R Commander. Requires

Example

We need a community matrix, which is a table of data that shows the number of individuals of each species in a community or at a site.

Table 1. A made up example matrix.

R code

library(vegan)
community_data <- data.frame(
Site1 = c(10, 5, 2, 0, 8),
Site2 = c(3, 12, 7, 1, 0),
Site3 = c(0, 4, 15, 6, 2)
)
# matrix transpose
community_matrix <- t(community_data)
H <- diversity(community_matrix, index="simpson")

R output

# returns 1 - D
    Site1     Site2     Site3 
0.6912000 0.6162571 0.6145405

# species richness

S <- specnumber(community_matrix)

R output

Site1 Site2 Site3 
    4     4     4

Calculate Pielou’s evenness index

J <- H/log(S)

R output

    Site1     Site2     Site3 
0.4985954 0.4445355 0.4432972

Compare two indexes.

Not uncommon, diversity indexes may be derived for a number of different locations and thus, logical to wish to compare diversity between sites. For example, diversity in disturbed compared to undisturbed habitats (eg, Hu et al 2022). While desirable, I will simply warn that such comparisons are not at all straight-forward and point interested readers to many resources to explore the statistical issues at hand (eg, Hulbert 1971, Jost 2006). To illustrate,  While it may be tempting to apply an independent sample t-test, you can’t. Shannon’s index is a single summary value, not a distribution; you only have one H value per site, not multiple samples (replicates) to calculate sample variance. Hutcheson’s t-test (Hutcheson 1970) can be used instead.

(6)  

and the variance of Shanon’s Index is given by

(7)  

Hutcheson’s t-test calculates the variance of the Shannon index itself, allowing you to compare two H values statistically. After calculating Shannon’s index (eg, vegan::diversity()), use R package ecolTest, function Hutcheson_t_test().

And, degrees of freedom? Yes, not simply the as is the case for the two-sample independent t-test. Degrees of freedom, , follow from the equation for the variance.

(8)  

With in hand, the critical value is found in the standard t-distribution table (Chapter 6.10 , Appendix).

Alternatively, in collecting data for diversity indexes, a better approach is to divide the area into quadrants and calculate diversity separately for each quadrant, hence, replication. See also Hill numbers, used to calculate effective number of species from traditional diversity indexes (Joost 2006). The effective number of species is the number of equally abundant species necessary to produce the observed value of diversity.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. Be able to define and distinguish between alpha diversity and beta diversity.

Quiz Chapter 20.8

Diversity indexes

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

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

Your Answers

References and suggested readings

Bollarapu, M. J., Kuchibhotla, S., Kvsn, R., & Patel, H. (2024). Dynamic perspectives on biodiversity quantification: Beyond conventional metrics. PeerJ, 12, e17924.

Hu, J., Amor, D. R., Barbier, M., Bunin, G., & Gore, J. (2022). Emergent phases of ecological diversity and dynamics mapped in microcosms. Science, 378(6615), 85–89.

Hutcheson, K. (1970). A test for comparing diversities based on Shannon formula. Journal of theoretical Biology, 29, 151-154.

Jost, L. (2006). Entropy and diversity. Oikos, 113(2), 363–375.

Morris, E. K., Caruso, T., Buscot, F., Fischer, M., Hancock, C., Maier, T. S., Meiners, T., Müller, C., Obermaier, E., Prati, D., Socher, S. A., Sonnemann, I., Wäschke, N., Wubet, T., Wurst, S., & Rillig, M. C. (2014). Choosing and using diversity indices: Insights for ecological applications from the German Biodiversity Exploratories. Ecology and Evolution, 4(18), 3514–3524.

Ricotta, C., & Szeidl, L. (2006). Towards a unifying approach to diversity measures: Bridging the gap between the Shannon entropy and Rao’s quadratic index. Theoretical Population Biology, 70(3), 237–243.

Roswell, M., Dushoff, J., & Winfree, R. (2021). A conceptual guide to measuring species diversity. Oikos, 130(3), 321–338.

Sarkar, S. (2021). Origin of the Term Biodiversity. BioScience, 71(9), 893.

Sîrbu, I., Benedek, A.-M., & Sîrbu, M. (2025). Rethinking composite quantification by capturing biological and ecological diversity across multiple dimensions. Scientific Reports, 15(1), 27822.

Tuomisto, H. (2010). A consistent terminology for quantifying species diversity? Yes, it does exist. Oecologia, 164(4), 853–860.

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD

20.6 – Dimensional analysis, clustering methods

Introduction

This chapter is just a sampling of the many techniques for “reducing the dimensions” of a data set that contains many variables, aka high-dimensional space, without eliminating the information contained among the many variables. The logic is similar to our discussions about model building — we may seek simple models at least in part so that we can see the big picture, not just each tree in an unhealthy forest.

Approaches can be categorized as feature extraction — a processing step to identify key characteristics contained in the larger data set into fewer, now transformed variables — or feature selection — the selection of a subset of the variables from among the full list which, presumably, contain redundant information. Principal Component Analysis is an example of feature extraction approaches. From this perspective, stepwise regression, discussed in Ch18.5, is an example of a feature selection approach. Clustering methods, however, groups objects based on existing features but does not create new features.

Dimensional analysis and clustering are complementary tools for exploring complex data (Ding and He 2004). Clustering methods use features — or their reduced representations — to group similar objects based on shared patterns. Dimensional analysis and clustering methods, presented together, show how we can first simplify biological complexity and then uncover natural groupings within it, whether we are analyzing candies by color and size or cancerous cells from normal cells by gene expression.

Principal component analysis

High dimensional data are reduced by linear transformation to sets of new vectors or coordinates. Each subsequent vector is orthogonal to the one proceeding it. By orthogonal we mean that variables are uncorrelated, completely independent of each other.

For example, consider body weight and height. Tall people tend to be heavier than shorter people, suggestive of the well-known correlation between the two — body weight and height are not orthogonal. For model building, one can’t simply include both as predictors (multicollinearity). Instead, we use PCA to transform body weight and height into a set of orthogonal variables called Principal Components. The first principal component, PC1, contains the largest amount of variation shared along the first linear adjustment, the second component, PC2, contains the next largest amount of variation, and so on. The result? Because PC2 is orthogonal to PC1, it contains completely different information — the weight that is not explained by height, or residuals.

We demonstrate utility of PCA with Bumpus data from MorphoFun/psa, variable names changed.

First, we plot (Fig 1) and add whether or not the individuals survived (TRUE) or not (FALSE). Figure 1 shows substantial overlap — there’s no distinct separation between the two groups — but survival group (blue) appears to be more concentrated in the lower-left to middle area of the plot. The non survival group (red) extends to higher masses and lengths, with several data points falling outside the main cluster and its corresponding ellipse.

Figure 1. Scatterplot English swallow mass (g) by total length (mm) by survival following winter storm

R code for Figure 1 graph.

car::scatterplot(Weight~Total_length | Survival, regLine=FALSE, 
      smooth=FALSE, boxplots=FALSE,
      ellipse=list(levels=c(.9)), 
      by.groups=TRUE, grid=FALSE, 
      pch=c(19,19), cex=1.5, col=c("red","blue"), 
      xlab="Total length, mm", ylab="Mass, g", data=Bumpus)

Bumpus measured several traits, we want to use all of the data. However, highly correlated (Fig 2) and therefore multicollinear. Manly (1986) and others (see Vieira 2012) have applied PCA to this data set.

Figure 2. Scatterplot matrix of Bumpus English sparrow traits. Traits were (left-right): Alar extent (mm), length (tip of beak to tip of tail), length of head (mm), length of femur (in.), length of humerus (in.), length of sternum (in.), skull width (in.), length of tibio-taurus (in.), and weight (g)

R code for graph

car::scatterplotMatrix(~Alar_extent+Beak_head_Length+Femur+Humerus+Keel_Length+Skull_width+Tibiotarsus+Total_length+Weight, regLine=FALSE, smooth=FALSE, diagonal=list(method="density"), data=Bumpus)

R packages

factoextra

psa package from MorphoFun/psa/

Rcmdr: Statistics > Dimensional analysis > Principal component analysis …

.PC <-
princomp(~Alar_extent+Beak_head_Length+Femur+Humerus+Keel_Length+Skull_width+Tibiotarsus+Total_length+Weight, cor=TRUE, data=Bumpus)
cat("\nComponent loadings:\n")
print(unclass(loadings(.PC)))
cat("\nComponent variances:\n")
print(.PC$sd^2)
cat("\n")
print(summary(.PC))
screeplot(.PC)
Bumpus <<- within(Bumpus, {
PC2 <- .PC$scores[,2]
PC1 <- .PC$scores[,1]
})
})

Importance of components

                              Comp.1     Comp.2
Standard deviation         2.3046882  0.9988978
Proportion of Variance     0.5901764  0.1108663
Cumulative Proportion      0.5901764  0.7010427

Interpretation

Comp.1 has a significantly higher standard deviation — it captures much more data variation than Comp.2.

Comp.1 explains 59% of the total variation, Comp.2 by itself describes an additional 11% of the variation in the data set. Between the two, 70% of the variation in the data set is explained by just these two components — “high data compression.”

Cluster analysis

Cluster analysis or clustering is a multivariate analysis technique that includes a number of different algorithms for grouping objects in such a way that objects in the same group (called a cluster) are more similar to each other than they are to objects in other groups. A number of approaches have been taken, but loosely can be grouped into distance clustering methods (see Chapter 16.6 – Similarity and Distance) and linkage clustering methods: Distance methods involve calculating the distance (or similarity) between two points and whereas linkage methods involve calculating distances among the clusters. Single linkage involves calculating the distance among all pairwise comparisons between two clusters, then

Results from cluster analyses are often displayed as dendrograms. Clustering methods include a number of different algorithms hierarchical clustering: single-linkage clustering; complete linkage clustering; average linkage clustering (UPGMA) centroid based clustering: k-means clustering.

Cluster analysis is common to molecular biology and phylogeny construction and more generally is an approach in use for exploratory data mining. Unsupervised machine learning (see 20.14 – Binary classification) used to classify, for example, methylation status of normal and diseased tissues from arrays (Clifford et al 2011).

K-means clustering

-means clustering is an unsupervised learning algorithm that automatically partitions data into -distinct, nonoverlapping groups (clusters) based on feature similarity. It finds hidden patterns by iteratively assigning data points to the nearest centroid (mean) and updating the centroids until clusters stabilize.

Ahh, the clusters. The -means problem is to find cluster centers that minimize the sum of squared distances from each data point being clustered to its cluster center (the center that is closest to it). Finding the optimal number of clusters is computationally hard, hard, because it requires finding the absolute best partition of points into clusters, which involves searching an exponential number of possible configurations.

Number of clusters

Iterations

Figure 3. Bi-plot of clusters, Skittles mini bags.

Hierarchical clustering

Ward’s method

Complete linkage

McQuitty’s method

Centroid linkage

A common way to depict the results of a cluster analysis is to construct a dendogram.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.

Quiz Chapter 20.6

Dimensional analysis

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

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

Your Answers

References and further reading

Bumpus, H. C. (1898). Eleventh lecture. The elimination of the unfit as illustrated by the introduced sparrow, Passer domesticus. (A fourth contribution to the study of variation.). Biology Lectures: Woods Hole Marine Biological Laboratory, 209–255.

Clifford, H., Wessely, F., Pendurthi, S., & Emes, R. D. (2011). Comparison of clustering methods for investigation of genome-wide methylation array data. Frontiers in genetics, 2, 88.

Ferreira, L., & Hitchcock, D. B. (2009). A comparison of hierarchical methods for clustering functional data. Communications in Statistics-Simulation and Computation, 38(9), 1925-1949.

Fraley C, Raftery AE. (2002) Model-based clustering, discriminant analysis, and density estimation. J Am Stat Assoc. 2002; 97(458):611–31.

Ding, C., & He, X. (2004). K-means Clustering via Principal Component Analysis. Proceedings of the 21 St International Conference on Machine Learning.

Vieira, V. M. N. C. S. (2012). Permutation tests to estimate significances on Principal Components Analysis. Computational Ecology and Software, 2.

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD