 |
| See Video ⮞ ☝ |
AGRON stats
November 29, 2020
Introduction
Correlation describes the direction and strength of association between variables. A pairwise analysis is bivariate; a matrix summarizes many pairs at once. These summaries are useful for exploration, but they do not by themselves identify causes or account for confounding variables.
Associations may be linear, monotonic but nonlinear, or more complex. Pearson correlation summarizes linear association; Spearman and Kendall methods use ranks and are appropriate for monotonic relationships under their respective assumptions.
- Positive association between the variables
- Negative association between the variables
- No association between the variables (fails to present any relationship)
The correlation coefficient is a standardized summary of association for a pair of variables. Its interpretation depends on the method, data range, outliers, and sampling design. A coefficient should be considered alongside a plot and its uncertainty, not as a complete description of the relationship.
For measuring linear relationships, Pearson correlation is the best measureThis tutorial demonstrates pairwise correlation tests and correlation matrices in R. We will also note what these analyses can and cannot support.
You may remove objects from the workspace if you are beginning a fresh analysis, but save any work you need first. Closing graphics devices is optional; clearing the console has no effect on the data or statistical results.
Importing data
Next step is importing the data set. The example used here is comprised of replication, two factor variables and rest of the variables are response variables. Load the package readxl by using library() function. To import the data from excel spreadsheet use read_excel() function. In argument path provide the link of the file. Type TRUE for argument col_names if the file contains first row as variable names. Your data are now available in the R console. They are stored in the object which I have chosen to call data.
Note that file paths are specified using slashes (/). In R, you cannot use backslashes (\), as you would in Microsoft Windows, unless you replace it with forward-slash (/) or double all the backslashes (\\). library(readxl)
data = read_excel(path = "E:/Correlation analysis/data_corr.xlsx",
col_names = TRUE)Visualize data
The example workbook contains replication, two treatment factors, and numeric response variables. Import it with readxl::read_excel(); set col_names = TRUE when the first row contains variable names. Confirm the column names, dimensions, data types, and missingness before analysis.
Use head() or tail() to inspect rows and View() to browse the data interactively. These functions display the data; they do not edit it.
The fix() function opens an editable data window. For reproducible work, make data corrections in the source file or with documented R transformations rather than relying on manual edits in a viewer.
Verify the structure of data
To verify the structure of the data, str() function is used.
# 'data.frame': 18 obs. of 10 variables:
# $ rep : num 1 1 1 1 1 1 2 2 2 2 ...
# $ method : chr "conv" "conv" "conv" "SRI" ...
# $ priming : chr "NP" "HP" "OP" "NP" ...
# $ em.head : num 88 92 83 85 86 82 87 89 84 86 ...
# $ head.mat : num 35 27 27 41 31 29 37 32 36 39 ...
# $ plant.height : num 113 126 118 123 120 127 116 128 119 116 ...
# $ total.tillers: num 445 498 535 497 534 562 450 498 538 452 ...
# $ th.kernal.wt : num 9.9 12 15.7 13.8 15.3 15.8 14.4 15.1 19.1 16.7 ...
# $ kernal.yield : num 3.46 3.58 3.61 3.28 3.79 3.73 3.34 3.65 4.03 3.63 ...
# $ straw.yield : num 12.9 13.3 13.9 11.2 11.3 12.3 13.5 14.2 14.2 11.7 ...Use str() to check the imported column types and confirm that the intended variables are numeric or categorical.
# 'data.frame': 18 obs. of 10 variables:
# $ rep : num 1 1 1 1 1 1 2 2 2 2 ...
# $ method : Factor w/ 2 levels "conv","SRI": 1 1 1 2 2 2 1 1 1 2 ...
# $ priming : Factor w/ 3 levels "HP","NP","OP": 2 1 3 2 1 3 2 1 3 2 ...
# $ em.head : num 88 92 83 85 86 82 87 89 84 86 ...
# $ head.mat : num 35 27 27 41 31 29 37 32 36 39 ...
# $ plant.height : num 113 126 118 123 120 127 116 128 119 116 ...
# $ total.tillers: num 445 498 535 497 534 562 450 498 538 452 ...
# $ th.kernal.wt : num 9.9 12 15.7 13.8 15.3 15.8 14.4 15.1 19.1 16.7 ...
# $ kernal.yield : num 3.46 3.58 3.61 3.28 3.79 3.73 3.34 3.65 4.03 3.63 ...
# $ straw.yield : num 12.9 13.3 13.9 11.2 11.3 12.3 13.5 14.2 14.2 11.7 ...If treatment columns have been imported as character strings, convert them to factors when they represent categories. Keep response measurements numeric and verify the conversion with str().
Correlation analysis
Using agricolae package
The original code uses attach() to access columns by name. Explicit data-frame references or a function’s data argument are safer for scripts because they reduce the risk of masked or ambiguous variable names.
We will examine the association between selected quantitative responses and then calculate a correlation matrix. Choose variables based on the research question and data scale; do not include treatment identifiers as continuous measurements without a defensible rationale.
library(agricolae)
cor.1 = correl(x = total.tillers,
y = kernal.yield,
method = "pearson",
alternative = "two.sided")
print(cor.1)# $stat
# [1] 4.419862
#
# $rho
# [1] 0.7414458
#
# $pvalue
# [1] 0.0004291801The tutorial uses agricolae for pairwise correlation tests. The function takes two numeric vectors and a method such as Pearson, Spearman, or Kendall. Select the method and alternative hypothesis before looking at the results, based on the variable scales and question.
Read the coefficient, test result, p-value, confidence interval where available, and number of complete pairs. Check a scatterplot for outliers, nonlinear patterns, and restricted ranges before interpreting the test.
# [1] 0.5497419
# [1] 0.07723301If the two-sided p-value is below the prespecified significance level, the data provide evidence against a zero correlation under the test assumptions. Report the coefficient and its uncertainty, not only the p-value. A positive estimate indicates that larger values of one variable tend to occur with larger values of the other in this sample; it does not show that one variable causes the other.
Using stats package
Squaring Pearson’s sample correlation gives r², a descriptive measure of the strength of the linear association in a simple two-variable setting. Its interpretation is not equivalent to a causal variance decomposition.
library(stats)
cor.2 = cor.test(x = total.tillers,
y = kernal.yield,
alternative = "two.sided",
method = "pearson",
conf.level = 0.95)
print(cor.2)#
# Pearson's product-moment correlation
#
# data: total.tillers and kernal.yield
# t = 4.4199, df = 16, p-value = 0.0004292
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
# 0.4199426 0.8976028
# sample estimates:
# cor
# 0.7414458Getting correlation matrix
Using stats package
Here, r² is about 0.55, meaning the fitted simple linear relationship is associated with roughly 55% of the sample variation in kernel yield around its mean. This does not mean tiller number caused that share of yield variation; other variables, design features, and measurement error may contribute.
require(stats)
cor.3 = cor(x = data[, c(4:10)], y = NULL,
use = "everything",
method = "pearson")
print(cor.3)# em.head head.mat plant.height total.tillers th.kernal.wt
# em.head 1.0000000 0.1545703 -0.40328717 -0.6078961974 -0.7123963
# head.mat 0.1545703 1.0000000 -0.33997627 -0.4923771076 -0.1018873
# plant.height -0.4032872 -0.3399763 1.00000000 0.5728506083 0.4487721
# total.tillers -0.6078962 -0.4923771 0.57285061 1.0000000000 0.5966086
# th.kernal.wt -0.7123963 -0.1018873 0.44877214 0.5966086434 1.0000000
# kernal.yield -0.6865255 -0.4265499 0.46877795 0.7414458464 0.7794977
# straw.yield 0.3121484 -0.2762724 -0.01128121 0.0004545191 -0.1204427
# kernal.yield straw.yield
# em.head -0.6865255 0.3121484037
# head.mat -0.4265499 -0.2762724365
# plant.height 0.4687779 -0.0112812050
# total.tillers 0.7414458 0.0004545191
# th.kernal.wt 0.7794977 -0.1204426847
# kernal.yield 1.0000000 0.1005602015
# straw.yield 0.1005602 1.0000000000Using agricolae package
The base R function cor.test() provides another route for testing a pairwise correlation. Specify the variables, method, alternative hypothesis, and confidence level as appropriate. Review the returned estimate, test statistic, degrees of freedom, p-value, and interval; available output depends on the method.
library(agricolae)
cor.4 = correlation(x = data[, c(4:10)], y = NULL,
method = "pearson",
alternative = "two.sided")
print(cor.4)# $correlation
# em.head head.mat plant.height total.tillers th.kernal.wt
# em.head 1.00 0.15 -0.40 -0.61 -0.71
# head.mat 0.15 1.00 -0.34 -0.49 -0.10
# plant.height -0.40 -0.34 1.00 0.57 0.45
# total.tillers -0.61 -0.49 0.57 1.00 0.60
# th.kernal.wt -0.71 -0.10 0.45 0.60 1.00
# kernal.yield -0.69 -0.43 0.47 0.74 0.78
# straw.yield 0.31 -0.28 -0.01 0.00 -0.12
# kernal.yield straw.yield
# em.head -0.69 0.31
# head.mat -0.43 -0.28
# plant.height 0.47 -0.01
# total.tillers 0.74 0.00
# th.kernal.wt 0.78 -0.12
# kernal.yield 1.00 0.10
# straw.yield 0.10 1.00
#
# $pvalue
# em.head head.mat plant.height total.tillers th.kernal.wt
# em.head 1.0000000000 0.54026966 0.09701049 0.0074434513 0.0009087942
# head.mat 0.5402696631 1.00000000 0.16747034 0.0379123294 0.6874722478
# plant.height 0.0970104894 0.16747034 1.00000000 0.0129596892 0.0617512187
# total.tillers 0.0074434513 0.03791233 0.01295969 1.0000000000 0.0089586967
# th.kernal.wt 0.0009087942 0.68747225 0.06175122 0.0089586967 1.0000000000
# kernal.yield 0.0016517901 0.07752274 0.04972248 0.0004291801 0.0001369208
# straw.yield 0.2072882086 0.26710622 0.96456396 0.9985718607 0.6340391279
# kernal.yield straw.yield
# em.head 0.0016517901 0.2072882
# head.mat 0.0775227402 0.2671062
# plant.height 0.0497224848 0.9645640
# total.tillers 0.0004291801 0.9985719
# th.kernal.wt 0.0001369208 0.6340391
# kernal.yield 1.0000000000 0.6913522
# straw.yield 0.6913521760 1.0000000
#
# $n.obs
# [1] 18A correlation matrix reports pairwise coefficients among a selected set of numeric variables. The cor() function computes the coefficients; choose the method and define how pairs with missing values are handled. Different missing-data options can produce different effective sample sizes and should be documented.
Using Hmisc package
The agricolae::correlation() function can return coefficients with corresponding test probabilities and pairwise sample sizes. Check which method and missing-data rules the function uses, and ensure the displayed entries align with the same variables and observations.
require(Hmisc)
cor.5 = rcorr(x = as.matrix(data[, c(4:10)]), y = NULL,
type ="pearson")
print(cor.5)# em.head head.mat plant.height total.tillers th.kernal.wt
# em.head 1.00 0.15 -0.40 -0.61 -0.71
# head.mat 0.15 1.00 -0.34 -0.49 -0.10
# plant.height -0.40 -0.34 1.00 0.57 0.45
# total.tillers -0.61 -0.49 0.57 1.00 0.60
# th.kernal.wt -0.71 -0.10 0.45 0.60 1.00
# kernal.yield -0.69 -0.43 0.47 0.74 0.78
# straw.yield 0.31 -0.28 -0.01 0.00 -0.12
# kernal.yield straw.yield
# em.head -0.69 0.31
# head.mat -0.43 -0.28
# plant.height 0.47 -0.01
# total.tillers 0.74 0.00
# th.kernal.wt 0.78 -0.12
# kernal.yield 1.00 0.10
# straw.yield 0.10 1.00
#
# n= 18
#
#
# P
# em.head head.mat plant.height total.tillers th.kernal.wt
# em.head 0.5403 0.0970 0.0074 0.0009
# head.mat 0.5403 0.1675 0.0379 0.6875
# plant.height 0.0970 0.1675 0.0130 0.0618
# total.tillers 0.0074 0.0379 0.0130 0.0090
# th.kernal.wt 0.0009 0.6875 0.0618 0.0090
# kernal.yield 0.0017 0.0775 0.0497 0.0004 0.0001
# straw.yield 0.2073 0.2671 0.9646 0.9986 0.6340
# kernal.yield straw.yield
# em.head 0.0017 0.2073
# head.mat 0.0775 0.2671
# plant.height 0.0497 0.9646
# total.tillers 0.0004 0.9986
# th.kernal.wt 0.0001 0.6340
# kernal.yield 0.6914
# straw.yield 0.6914Properties
Access named components of a returned R object with
$or, preferably for programming, with[[ ]]. Inspect the object structure before selecting a component so the code refers to the intended matrix.The
Hmisc::rcorr()function provides Pearson or Spearman correlation matrices, sample sizes, and p-values for numeric matrix input. Prepare the input carefully, because nonnumeric columns and missing-value handling affect the returned results.A pairwise correlation is symmetric:
cor(x, y)equalscor(y, x)for the same observations and method.Correlation is unit-free and can summarize standardized association between measurements recorded in different units. It does not make the variables scientifically comparable in every other respect.
For Pearson correlation, the coefficient ranges from −1 to +1. Values near either endpoint indicate a strong linear pattern; values near zero indicate little linear association, which does not rule out a nonlinear relationship.
Adding a constant to either variable does not change Pearson correlation. Multiplying by a positive constant leaves it unchanged; multiplying by a negative constant reverses its sign.
Limitations
The pairwise coefficient is the same whichever variable is written first:
r(x, y) = r(y, x).Independence implies zero correlation when the relevant moments exist, but zero correlation does not generally imply independence. A nonlinear relationship can have zero Pearson correlation.
Pearson correlation measures linear association. Rank correlations describe monotonic association; neither coefficient captures every possible nonlinear pattern. Plot the data before selecting a summary.
Please comment below if you have any questions.With several variables, pairwise correlations do not adjust for the others. They may therefore differ from conditional or multivariable relationships, and a matrix can contain many tests requiring attention to multiplicity.
A strong pairwise association may arise because a third variable influences both measured variables or because of selection, design, or measurement structure. Use an appropriate design and multivariable analysis to address confounding; correlation alone cannot identify its source.
Download R program —
Click_hereDownload R studio —
Click_here
Thank you so much for this excellent tutorial and explanation.
ReplyDelete