A common task with matrix data is to fit the same test or model to
every row or every column: which genes differ between conditions, which
survey questions are answered differently by different groups.
tidymatrix does this with compute_across() and a set of
wrappers around standard tests.
The active component decides the direction. With columns active, the function is applied to each column, and the row metadata supplies the variables to test against. With rows active it is the other way round.
library(tidymatrix)
library(dplyr, warn.conflicts = FALSE)Data
We use the big5 personality survey (see
?big5), with careless respondents removed and reverse-keyed
items re-scored, so that a higher value always means more of the trait
(see Matrix operations).
tm <- tidymatrix(big5_responses, big5_respondents, big5_items) |>
activate(rows) |>
filter(completion_min > 3.5) |>
transform_matrix(\(x, flip) ifelse(flip, 6L - x, x), flip = reversed)Comparing two groups
t-test
Do women and men answer the items differently? The survey also has a few non-binary respondents, too few for a separate group, so we compare women and men only:
tm_fm <- tm |>
activate(rows) |>
filter(gender %in% c("Female", "Male"))
ttest <- tm_fm |>
activate(columns) |>
compute_ttest(group_col = "gender", control = "Male", treatment = "Female")
ttest |>
select(item_id, trait, p.value, log2fc, p.adj) |>
arrange(p.value)
#> # A tibble: 30 × 5
#> item_id trait p.value log2fc p.adj
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 N3 Neuroticism 0.00000154 0.272 0.0000303
#> 2 N5 Neuroticism 0.00000202 0.276 0.0000303
#> 3 N1 Neuroticism 0.0000535 0.229 0.000535
#> 4 N4 Neuroticism 0.0000717 0.218 0.000538
#> 5 N6 Neuroticism 0.000428 0.168 0.00257
#> 6 A3 Agreeableness 0.00103 0.170 0.00514
#> 7 N2 Neuroticism 0.00424 0.165 0.0182
#> 8 A6 Agreeableness 0.00950 0.123 0.0356
#> 9 A4 Agreeableness 0.0145 0.106 0.0485
#> 10 A5 Agreeableness 0.0240 0.134 0.0694
#> # ℹ 20 more rowsThe result has one row per item: the item metadata, followed by the
p-value, the log2 ratio of the group means
(treatment / control) and the FDR-adjusted p-value. If you
leave out control and treatment, they are
inferred — in the order of the factor levels, or alphabetically for
other columns — and a message reports which comparison was made. Setting
them explicitly documents the direction of the fold change and also lets
you compare two groups out of several without filtering first.
Which traits do the significant items belong to?
ttest |>
filter(p.adj < 0.05) |>
count(trait)
#> # A tibble: 2 × 2
#> trait n
#> <chr> <int>
#> 1 Agreeableness 3
#> 2 Neuroticism 6Women score higher on agreeableness and neuroticism items — the same pattern as in real personality data.
Arguments not used by compute_ttest() go to
t.test(), e.g. var.equal = TRUE. Use
adjust to choose another p.adjust() method, or
adjust = "none" to skip the adjustment.
Wilcoxon test
Likert answers are ordinal, so a rank-based test is a reasonable
alternative. compute_wilcox() has the same interface and
adds the difference in medians:
tm_fm |>
activate(columns) |>
compute_wilcox(group_col = "gender", control = "Male", treatment = "Female") |>
filter(p.adj < 0.05) |>
select(item_id, trait, median_diff, p.value, p.adj)
#> # A tibble: 8 × 5
#> item_id trait median_diff p.value p.adj
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 A3 Agreeableness 0 0.00168 0.00840
#> 2 A6 Agreeableness 1 0.00916 0.0343
#> 3 N1 Neuroticism 0 0.0000525 0.000486
#> 4 N2 Neuroticism 0 0.00348 0.0149
#> 5 N3 Neuroticism 0 0.00000160 0.0000481
#> 6 N4 Neuroticism 0 0.0000648 0.000486
#> 7 N5 Neuroticism 0 0.00000466 0.0000698
#> 8 N6 Neuroticism 1 0.000404 0.00242Continuous variables
Correlation
compute_correlation() correlates each item with a
numeric metadata variable (method can be
"pearson", "spearman" or
"kendall"; note that the rank-based methods warn about ties
with Likert data). Personality changes slowly over adulthood:
age_cor <- tm |>
activate(columns) |>
compute_correlation(var = "age")
age_cor |>
filter(p.adj < 0.05) |>
select(item_id, trait, correlation, p.adj) |>
arrange(correlation)
#> # A tibble: 16 × 4
#> item_id trait correlation p.adj
#> <chr> <chr> <dbl> <dbl>
#> 1 N4 Neuroticism -0.236 0.0000436
#> 2 N2 Neuroticism -0.203 0.000249
#> 3 N5 Neuroticism -0.198 0.000313
#> 4 N1 Neuroticism -0.196 0.000328
#> 5 N3 Neuroticism -0.193 0.000407
#> 6 E5 Extraversion -0.118 0.0377
#> 7 A6 Agreeableness 0.124 0.0294
#> 8 A4 Agreeableness 0.125 0.0292
#> 9 A2 Agreeableness 0.139 0.0137
#> 10 A3 Agreeableness 0.164 0.00295
#> 11 C4 Conscientiousness 0.174 0.00157
#> 12 C5 Conscientiousness 0.203 0.000249
#> 13 C2 Conscientiousness 0.211 0.000162
#> 14 C3 Conscientiousness 0.218 0.000112
#> 15 C6 Conscientiousness 0.228 0.0000581
#> 16 C1 Conscientiousness 0.235 0.0000436Neuroticism items correlate negatively with age, conscientiousness and agreeableness items positively.
Simple regression
compute_lm_simple() fits value ~ predictor
and reports the slope, intercept and R²:
tm |>
activate(columns) |>
compute_lm_simple(predictor = "life_satisfaction") |>
select(item_id, trait, slope, r.squared, p.adj) |>
arrange(slope)
#> # A tibble: 30 × 5
#> item_id trait slope r.squared p.adj
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 N4 Neuroticism -0.286 0.239 3.46e-23
#> 2 N3 Neuroticism -0.272 0.212 1.46e-20
#> 3 N5 Neuroticism -0.258 0.199 2.43e-19
#> 4 N1 Neuroticism -0.254 0.183 6.47e-18
#> 5 N2 Neuroticism -0.246 0.185 5.17e-18
#> 6 N6 Neuroticism -0.194 0.113 5.85e-11
#> 7 O5 Openness 0.0259 0.00202 3.77e- 1
#> 8 O3 Openness 0.0341 0.00325 2.72e- 1
#> 9 O2 Openness 0.0363 0.00414 2.21e- 1
#> 10 O6 Openness 0.0458 0.00609 1.39e- 1
#> # ℹ 20 more rowsMultiple regression
compute_lm() takes the right-hand side of a model
formula, so you can adjust for other variables. With coef
set, it reports that coefficient; without it, the overall fit of the
model.
Is the gender difference in neuroticism still there after adjusting
for age and education? We make Male the reference level so
that the coefficient is called genderFemale:
tm_fm <- tm_fm |>
activate(rows) |>
mutate(gender = factor(gender, levels = c("Male", "Female")))
tm_fm |>
activate(columns) |>
filter(trait == "Neuroticism") |>
compute_lm(~ age + gender + education, coef = "genderFemale") |>
select(item_id, item_text, estimate, se, p.adj)
#> # A tibble: 6 × 5
#> item_id item_text estimate se p.adj
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 N1 I get stressed out easily. 0.461 0.119 0.000223
#> 2 N2 I worry about things. 0.323 0.114 0.00502
#> 3 N3 My mood changes often. 0.564 0.117 0.0000126
#> 4 N4 I get irritated easily. 0.444 0.116 0.000223
#> 5 N5 I stay calm under pressure. 0.533 0.114 0.0000127
#> 6 N6 I seldom feel blue. 0.385 0.118 0.00146
tm_fm |>
activate(columns) |>
compute_lm(~ age + gender + education) |>
select(item_id, trait, r.squared, p.value) |>
arrange(desc(r.squared)) |>
head()
#> # A tibble: 6 × 4
#> item_id trait r.squared p.value
#> <chr> <chr> <dbl> <dbl>
#> 1 N5 Neuroticism 0.103 0.000000310
#> 2 N4 Neuroticism 0.102 0.000000322
#> 3 N3 Neuroticism 0.102 0.000000368
#> 4 O5 Openness 0.0919 0.00000237
#> 5 O4 Openness 0.0913 0.00000264
#> 6 N1 Neuroticism 0.0912 0.00000268More than two groups
compute_anova() and its rank-based counterpart
compute_kruskal() compare several groups. Openness
increases with education:
tm |>
activate(columns) |>
compute_anova(group_col = "education") |>
filter(p.adj < 0.05) |>
select(item_id, trait, f.statistic, p.adj)
#> # A tibble: 6 × 4
#> item_id trait f.statistic p.adj
#> <chr> <chr> <dbl> <dbl>
#> 1 O1 Openness 3.92 0.0197
#> 2 O2 Openness 4.61 0.00721
#> 3 O3 Openness 5.96 0.00117
#> 4 O4 Openness 9.06 0.0000159
#> 5 O5 Openness 7.22 0.000195
#> 6 O6 Openness 4.62 0.00721
tm |>
activate(columns) |>
compute_kruskal(group_col = "occupation") |>
filter(p.adj < 0.05) |>
select(item_id, trait, statistic, p.adj)
#> # A tibble: 15 × 4
#> item_id trait statistic p.adj
#> <chr> <chr> <dbl> <dbl>
#> 1 A2 Agreeableness 16.5 0.0239
#> 2 A4 Agreeableness 15.2 0.0368
#> 3 C1 Conscientiousness 26.5 0.000659
#> 4 C2 Conscientiousness 30.1 0.000381
#> 5 C3 Conscientiousness 39.3 0.0000191
#> 6 C4 Conscientiousness 26.3 0.000659
#> 7 C5 Conscientiousness 27.5 0.000528
#> 8 C6 Conscientiousness 27.8 0.000528
#> 9 N1 Neuroticism 19.2 0.00891
#> 10 N2 Neuroticism 22.8 0.00232
#> 11 N3 Neuroticism 24.3 0.00136
#> 12 N4 Neuroticism 31.1 0.000362
#> 13 N5 Neuroticism 27.4 0.000528
#> 14 O4 Openness 28.3 0.000528
#> 15 O5 Openness 22.2 0.00277Your own function: compute_across()
All the wrappers above are built on compute_across(). It
takes a function of two arguments — the values of one row or column, and
the full metadata of the other dimension — that returns a named list or
vector. Each element becomes a column of the result.
For example, an effect size is often more useful than a p-value. Cohen’s d for the gender difference:
cohens_d <- function(values, respondents) {
f <- values[respondents$gender == "Female"]
m <- values[respondents$gender == "Male"]
pooled_sd <- sqrt(
((length(f) - 1) * var(f) + (length(m) - 1) * var(m)) /
(length(f) + length(m) - 2)
)
list(
mean_female = mean(f),
mean_male = mean(m),
d = (mean(f) - mean(m)) / pooled_sd
)
}
effect_sizes <- tm_fm |>
activate(columns) |>
compute_across(cohens_d)
effect_sizes |>
group_by(trait) |>
summarise(mean_d = mean(d), min_d = min(d), max_d = max(d))
#> # A tibble: 5 × 4
#> trait mean_d min_d max_d
#> <chr> <dbl> <dbl> <dbl>
#> 1 Agreeableness 0.247 0.153 0.340
#> 2 Conscientiousness 0.0701 -0.0316 0.148
#> 3 Extraversion 0.0862 -0.0300 0.220
#> 4 Neuroticism 0.416 0.294 0.500
#> 5 Openness -0.0916 -0.157 -0.0690Extra arguments to compute_across() are passed on to the
function, and errors in individual rows or columns become warnings
rather than stopping the whole computation.
Row-wise statistics
With rows active the function gets one respondent’s answers together with the item metadata. A classic check for careless responding is the longstring index: the longest run of identical answers in the order the questions were asked. On the raw data, before removing anyone:
longstring <- function(answers, items) {
runs <- rle(answers[order(items$position)])
list(longstring = max(runs$lengths))
}
tidymatrix(big5_responses, big5_respondents, big5_items) |>
activate(rows) |>
compute_across(longstring) |>
arrange(desc(longstring)) |>
select(respondent_id, longstring, completion_min) |>
head(8)
#> # A tibble: 8 × 3
#> respondent_id longstring completion_min
#> <chr> <int> <dbl>
#> 1 R352 27 2.6
#> 2 R074 24 2
#> 3 R170 20 2.5
#> 4 R219 20 2.2
#> 5 R307 20 3.2
#> 6 R392 16 1.5
#> 7 R049 9 15
#> 8 R176 8 10The respondents who gave the same answer over and over are the ones who finished in two or three minutes.
The function can return several values at once. Here, one score per trait:
tm |>
activate(rows) |>
compute_across(\(answers, items) tapply(answers, items$trait, mean)) |>
select(respondent_id, Agreeableness:Openness)
#> # A tibble: 388 × 6
#> respondent_id Agreeableness Conscientiousness Extraversion Neuroticism
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 R001 3.33 1.5 3.33 3
#> 2 R002 2.83 2.67 2.5 4.5
#> 3 R003 3.33 2.17 2.67 3.17
#> 4 R004 2.83 2.5 2.67 4.17
#> 5 R005 4.83 3.67 1.33 3.33
#> 6 R006 4 2.17 2.5 5
#> 7 R007 4.83 3.17 3 4.17
#> 8 R008 3 2.33 2.67 5
#> 9 R009 4.5 3.83 2.67 3
#> 10 R010 3.67 4.67 4.17 2.33
#> # ℹ 378 more rows
#> # ℹ 1 more variable: Openness <dbl>Storing results in the metadata
By default the result is returned as a separate tibble. With
add_to_data = TRUE it is added to the active metadata
instead, and the tidymatrix is returned. Use prefix to keep
the columns of different analyses apart — compute_across()
refuses to overwrite existing columns:
tm_stats <- tm_fm |>
activate(columns) |>
compute_ttest(
group_col = "gender", control = "Male", treatment = "Female",
add_to_data = TRUE, prefix = "gender"
) |>
compute_correlation(var = "age", add_to_data = TRUE, prefix = "age") |>
compute_across(cohens_d, add_to_data = TRUE, prefix = "gender")
tm_stats |>
activate(columns) |>
select(item_id, trait, gender_p.adj, gender_d, age_correlation, age_p.adj)
#> # A tidymatrix: 382 x 30 matrix
#> # Active: columns
#> #
#> # Row data: 382 rows x 8 columns
#> # Column data: 30 rows x 6 columns
#> #
#> # Active data (columns):
#> item_id trait gender_p.adj gender_d age_correlation age_p.adj
#> 1 E1 Extraversion 0.62744096 0.07339704 -0.05473334 0.32994274
#> 2 E2 Extraversion 0.27020931 0.16225759 -0.11856460 0.03835155
#> 3 E3 Extraversion 0.76012764 0.04717771 -0.08161302 0.13907874
#> 4 E4 Extraversion 0.76012764 0.04440616 -0.10341683 0.06196758
#> 5 E5 Extraversion 0.77169042 -0.02999294 -0.11444498 0.04216239
#> 6 E6 Extraversion 0.08418223 0.21970081 -0.10444446 0.06196758Now the statistics can be used like any other metadata, for example to keep only the items that differ between women and men:
tm_stats |>
activate(columns) |>
filter(gender_p.adj < 0.05, abs(gender_d) > 0.2)
#> # A tidymatrix: 382 x 9 matrix
#> # Active: columns
#> #
#> # Row data: 382 rows x 8 columns
#> # Column data: 9 rows x 14 columns
#> #
#> # Active data (columns):
#> item_id trait reversed
#> 1 A3 Agreeableness FALSE
#> 2 A4 Agreeableness FALSE
#> 3 A6 Agreeableness TRUE
#> 4 N1 Neuroticism FALSE
#> 5 N2 Neuroticism FALSE
#> 6 N3 Neuroticism FALSE
#> item_text position gender_p.value
#> 1 I trust what people tell me. 12 1.028150e-03
#> 2 I am quick to forgive. 17 1.454847e-02
#> 3 I am not really interested in others' problems. 27 9.500961e-03
#> 4 I get stressed out easily. 4 5.349460e-05
#> 5 I worry about things. 9 4.240518e-03
#> 6 My mood changes often. 14 1.539203e-06
#> gender_log2fc gender_p.adj age_correlation age_p.value age_p.adj
#> 1 0.1699653 5.140752e-03 0.1729122 6.886390e-04 0.0018781064
#> 2 0.1055770 4.849489e-02 0.1189164 2.008060e-02 0.0383515535
#> 3 0.1227857 3.562861e-02 0.1209797 1.800708e-02 0.0383515535
#> 4 0.2292184 5.349460e-04 -0.2056356 5.137184e-05 0.0002201650
#> 5 0.1647665 1.817365e-02 -0.2103225 3.417688e-05 0.0001733995
#> 6 0.2724246 3.025276e-05 -0.1965773 1.100582e-04 0.0003301747
#> gender_mean_female gender_mean_male gender_d
#> 1 3.466667 3.081395 0.3402807
#> 2 3.809524 3.540698 0.2572107
#> 3 3.747619 3.441860 0.2697998
#> 4 3.366667 2.872093 0.4252665
#> 5 3.076190 2.744186 0.2943243
#> 6 3.328571 2.755814 0.5002835See also
- Plotting and exporting for a volcano plot of these results.
- PCA and clustering for unsupervised analyses.