Correlation matrix in R with p-values and intervals

A correlation matrix is a table of the correlation between every pair of your variables. In R, cor() gives you the coefficients in one line, but a table of coefficients alone says nothing about uncertainty and nothing about how many pairs you looked at. This post builds the matrix with a p-value and a confidence interval for every pair, adds an adjustment for testing many pairs at once, and shows how to read the result.

The data are the swiss data that ship with R. One pair from them is rerun through the deterministic engine so you can check it against base R. Everything else is printed R code that you can run yourself, and the post says which is which.

The data and the question

swiss holds a standardised fertility measure and five socioeconomic indicators for 47 French-speaking Swiss provinces around 1888, one row per province, all six variables in percent (R documentation for swiss, after Mosteller and Tukey, 1977). The documentation describes the period as the start of Switzerland’s demographic transition, when fertility was beginning to fall from the high level typical of less developed countries.

The six variables are Fertility, Agriculture (share of males working in agriculture), Examination (share of army draftees with the highest examination mark), Education (share of draftees with schooling beyond primary), Catholic (share of the population that is Catholic as opposed to Protestant) and Infant.Mortality (share of live births who die before age one).

Six variables make 15 distinct pairs, and that number is the whole reason this post exists. Each pair is its own hypothesis test.

What the engine ran

The engine’s Pearson method takes one pair of variables, not a whole matrix. So the deterministic engine ran one cell of the matrix: Fertility against Education, on the 47 provinces.

• Pearson’s r = -0.6638, 95% CI -0.7987 to -0.4653
• t = -5.954, p < 0.001, n = 47
• Engine assumption check: sample size adequate for the method

The engine’s output for this method includes no normality check, so the shape of the data is left to you. The next section does that work in base R. Everything from here on, apart from this one cell, is base R code and its printed output, not engine output.

The matrix in base R

First look at the data, because a correlation coefficient hides the shape of the scatter. Anscombe (1973) built four data sets with the same correlation and very different pictures, which is the standing reason to draw before you compute.

pairs(swiss, panel = panel.smooth)   # 15 scatterplots; look before you test
r <- round(cor(swiss), 2)            # Pearson, the default method
r

The coefficient matrix, rounded to two decimals:

	Fertility	Agriculture	Examination	Education	Catholic	Infant.Mortality
Fertility	1	0.35	-0.65	-0.66	0.46	0.42
Agriculture	0.35	1	-0.69	-0.64	0.40	-0.06
Examination	-0.65	-0.69	1	0.70	-0.57	-0.11
Education	-0.66	-0.64	0.70	1	-0.15	-0.10
Catholic	0.46	0.40	-0.57	-0.15	1	0.18
Infant.Mortality	0.42	-0.06	-0.11	-0.10	0.18	1

The matrix is symmetric, and the diagonal is always 1, so only the 15 cells on one side of the diagonal carry information. Fertility against Education, the engine’s pair, is -0.66, matching the engine’s -0.6638 once rounded.

By default cor() returns NA for any pair that involves a missing value (R documentation for cor). The swiss data have no missing values, so that does not bite here. In your own data you must choose. use = "complete.obs" drops every row with any missing value, so all cells share one sample. use = "pairwise.complete.obs" keeps all available pairs for each cell, but the documentation warns that it can produce a matrix that is not positive semi-definite, and the cells then rest on different sample sizes.

A p-value and an interval for every pair

cor() gives no p-values or intervals. cor.test() does, one pair at a time, and looping over all 15 pairs collects them. The interval for a Pearson correlation is an asymptotic one based on Fisher’s z transformation (R documentation for cor.test; Fisher, 1915).

vars  <- names(swiss)
pairs <- combn(vars, 2, simplify = FALSE)

res <- do.call(rbind, lapply(pairs, function(p) {
  t <- cor.test(swiss[[p[1]]], swiss[[p[2]]])
  data.frame(a = p[1], b = p[2], r = unname(t$estimate),
             lo = t$conf.int[1], hi = t$conf.int[2], p = t$p.value)
}))

res$p_holm <- p.adjust(res$p, method = "holm")
res$p_bonferroni <- p.adjust(res$p, method = "bonferroni")

The result, with the intervals unadjusted (each is a 95% interval for that pair on its own):

Pair	r	95% CI	p	Holm p
Fertility, Agriculture	0.35	0.07 to 0.58	0.015	0.090
Fertility, Examination	-0.65	-0.79 to -0.44	< 0.001	< 0.001
Fertility, Education	-0.66	-0.80 to -0.47	< 0.001	< 0.001
Fertility, Catholic	0.46	0.20 to 0.66	0.001	0.009
Fertility, Infant.Mortality	0.42	0.15 to 0.63	0.004	0.029
Agriculture, Examination	-0.69	-0.81 to -0.50	< 0.001	< 0.001
Agriculture, Education	-0.64	-0.78 to -0.43	< 0.001	< 0.001
Agriculture, Catholic	0.40	0.13 to 0.62	0.005	0.036
Agriculture, Infant.Mortality	-0.06	-0.34 to 0.23	0.684	1.000
Examination, Education	0.70	0.51 to 0.82	< 0.001	< 0.001
Examination, Catholic	-0.57	-0.74 to -0.34	< 0.001	< 0.001
Examination, Infant.Mortality	-0.11	-0.39 to 0.18	0.445	1.000
Education, Catholic	-0.15	-0.42 to 0.14	0.302	1.000
Education, Infant.Mortality	-0.10	-0.38 to 0.19	0.507	1.000
Catholic, Infant.Mortality	0.18	-0.12 to 0.44	0.238	1.000

Report the interval, not only the p-value. Fertility against Examination is -0.65 with an interval from -0.79 to -0.44: the data are compatible with a moderate to strong negative association in these provinces, and the interval shows how precisely 47 provinces pin it down. Agriculture with Infant.Mortality has an interval from -0.34 to 0.23, which is compatible with anything from a moderate negative to a weak positive correlation.

Why a grid of stars is not enough

Many tutorials print the matrix with an asterisk beside every cell whose p-value is under 0.05. That treats 15 tests as if each were the only one. If all 15 true correlations were zero and the tests were independent, the chance of at least one p-value under 0.05 would be 1 - 0.95^15, about 54 percent. The 15 tests here are not independent, because the pairs share variables, so that figure is a guide to the size of the problem and not an exact rate.

In the table above, 10 of the 15 raw p-values are below 0.05. The adjustment changes the count. Our multiple comparisons guide covers the families of methods in more depth; the short version for a matrix is:

• Bonferroni multiplies each p-value by the number of tests, here 15. It controls the chance of any false positive, but it is the most conservative option. Seven pairs stay below 0.05.
• Holm (1979) is a step-down version that controls the same error rate. The R documentation says it is valid under arbitrary assumptions and dominates the unmodified Bonferroni method. Nine pairs stay below 0.05.
• Benjamini-Hochberg (1995), method = "BH", controls the expected proportion of false discoveries among the pairs you flag, a weaker guarantee than controlling any false positive. Benjamini and Hochberg proved it for independent tests. The pairs in a correlation matrix share variables, so the tests are not independent; if you want FDR control under any dependence, the R documentation offers method = "BY", which is more conservative.

Holm is a reasonable default for an exploratory matrix because it needs no assumption about how the tests relate. The two pairs where Holm and Bonferroni disagree are visible in the table. Fertility with Infant.Mortality has a raw p of 0.004 and a Holm p of 0.029, which clears 0.05, while its Bonferroni p is 0.054 and does not. Agriculture with Catholic is the same story: Holm 0.036, Bonferroni 0.078. Fertility with Agriculture shows the cost of any adjustment: raw p 0.015, Holm p 0.090.

The intervals in the table are not adjusted. Each 95% interval is built to cover its own true value with 95% probability, and all 15 together cover with a lower probability. If you need simultaneous coverage, say so and adjust the intervals too; psych::corr.test() reports both raw and adjusted intervals (psych documentation).

Reading the table

A correlation matrix is easy to over-read. Four points matter here.

The unit is a province. These are percentages for provinces, not measurements on people. A correlation between provincial rates does not tell you how individuals behave, and the gap between the two can be large (Robinson, 1950). Fertility and Education moving together across provinces is a statement about provinces. Neighbouring provinces may also resemble each other, which makes the 47 rows less independent than 47 separate draws.

Pearson measures a straight line. An r near zero does not mean no relationship, only no straight-line one, and a single extreme point can create or erase a correlation. The scatterplots from pairs() are the check. Catholic is a good case: its median is 15.14 percent and its mean 41.14, and 35 of the 47 provinces are either below 10 percent or above 90 percent. The variable is close to two groups, not a smooth continuum, so a single r for it describes a contrast between two kinds of province.

Pearson’s p-value and interval lean on normality. The cor.test documentation ties the Pearson test to samples that follow normal distributions, and the cor documentation says the rank methods are more robust when the data may not come from a bivariate normal distribution. Simulation work found that with clearly non-normal data the usual Pearson test can give inflated false positives and lose power (Bishara and Hittner, 2012). For a variable like Catholic, run the rank version too:

cor(swiss, method = "spearman")
cor.test(swiss$Fertility, swiss$Catholic, method = "spearman", exact = FALSE)

For Fertility and Catholic, Spearman’s rho is 0.41 (p = 0.004) against Pearson’s 0.46, so the direction and rough size survive the change of method. The cor.test documentation notes that it gives a confidence interval only for the Pearson method.

A matrix does not separate causes. Education and Examination both fall as Agriculture rises, and they rise together, and a matrix cannot untangle which variable drives which. If the question is about one outcome and several predictors at once, the next step is a regression, covered in the multiple regression guide.

The shortcut with psych

If you want all of this from one call, psych::corr.test() returns the coefficients, the sample size for each cell, raw and adjusted p-values and intervals. Its documented defaults are use = "pairwise", method = "pearson" and adjust = "holm", and for a symmetric matrix the adjusted p-values sit above the diagonal and the raw ones below it (psych documentation).

library(psych)
ct <- corr.test(swiss)          # Holm-adjusted by default
print(ct, short = FALSE)        # includes the confidence intervals

Know what the defaults do before you trust them. The call uses pairwise deletion, so with missing data the cells can rest on different samples. State the adjustment you used, because a reader who sees p-values above and below the diagonal needs to know which is which.

Reporting a correlation matrix

A short methods sentence is enough: “We computed Pearson correlations between all pairs of the six variables (n = 47 provinces). Two-sided p-values were adjusted for 15 comparisons with Holm’s method; intervals are 95% intervals for each pair and are not adjusted.” Then give a table like the one above, with r, the interval and the adjusted p-value.

In the text, lead with the strongest findings and give their intervals: Fertility was lower where Education was higher across these provinces (r = -0.66, 95% CI -0.80 to -0.47). Say that the data are descriptive and at province level. Keep the exact p-values for the table, and report an adjusted value as adjusted. Whatever you decide about the adjustment, make the decision before you look at which cells are significant.

Rerunning the example

The script in the engine package reproduces the Fertility and Education correlation and the attached sensitivity power calculation, from the bundled data. It does not build the matrix, the 15 tests, the Holm and Bonferroni adjustments or the Spearman comparison; those are the base R code above, run on R 4.6.1, the same version as the engine image, and are not engine output. The reproducible analysis guide explains what such a package contains, and the Pearson and Spearman reference pages list the engine’s inputs and outputs. For choosing between the two, the test chooser walks through the decision.

This example is about the method, not about Switzerland or fertility: the 1888 data are public, small and chosen only because they ship with R.

# Results

- Test: Pearson correlation
- n: 47
- statistic: -5.954
- p-value: < 0.001
- effect size: Pearson's r = -0.6638
- 95% CI: [-0.7987, -0.4653]
- seed: 20261023

## Assumption checks
- sample size: passed=TRUE (adequate n)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=47: correlation r = 0.3957
- to detect a conventional medium effect: 85 total to detect r=0.30
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Pearson correlation  |  engine 0.31.1  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261023)
.data_path <- "data/correlation-matrix-in-r.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("Fertility", "Education")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- cor.test(data[["Fertility"]], data[["Education"]], method = "pearson")
print(result)

# Sensitivity: minimum detectable effect at n=47 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.r.test(n = 47, sig.level = 0.05, power = 0.8)
print(power)