spicy is an R package for descriptive statistics and data analysis,
designed for data science and survey research workflows. It covers
variable inspection, frequency tables, cross-tabulations with
chi-squared tests and effect sizes, and publication-ready summary
tables, offering functionality similar to Stata or SPSS but within a
tidyverse-friendly R environment. This vignette walks through the core
workflow using the bundled sochealth dataset, a simulated
social-health survey with 1200 respondents and 29 variables.
varlist() (or its shortcut vl()) gives a
compact overview of every variable in a data frame: name, label,
representative values, class, number of distinct values, valid
observations, and missing values. In RStudio or Positron, calling
varlist(mydata) with the default tbl = FALSE
opens an interactive viewer - this is the most common usage in practice.
Here we use tbl = TRUE to produce static output for the
vignette:
varlist(sochealth, tbl = TRUE)
#> # A tibble: 29 × 7
#> Variable Label Values Class N_distinct N_valid NAs
#> <chr> <chr> <chr> <chr> <int> <int> <int>
#> 1 sex Sex Femal… fact… 2 1200 0
#> 2 age Age (years) 25, 2… nume… 51 1200 0
#> 3 age_group Age group 25-34… orde… 4 1200 0
#> 4 education Highest education le… Lower… orde… 3 1200 0
#> 5 social_class Subjective social cl… Lower… orde… 5 1200 0
#> 6 region Region of residence Centr… fact… 6 1200 0
#> 7 employment_status Employment status Emplo… fact… 4 1200 0
#> 8 income_group Household income gro… Low, … orde… 4 1182 18
#> 9 income Monthly household in… 1000,… nume… 1052 1200 0
#> 10 smoking Current smoker No, Y… fact… 2 1175 25
#> # ℹ 19 more rowsYou can also select specific columns with tidyselect syntax:
varlist(sochealth, starts_with("bmi"), income, weight, tbl = TRUE)
#> # A tibble: 4 × 7
#> Variable Label Values Class N_distinct N_valid NAs
#> <chr> <chr> <chr> <chr> <int> <int> <int>
#> 1 bmi Body mass index 16, 1… nume… 177 1188 12
#> 2 bmi_category BMI category Norma… orde… 3 1188 12
#> 3 income Monthly household income (… 1000,… nume… 1052 1200 0
#> 4 weight Survey design weight 0.294… nume… 794 1200 0code_book() turns the same inspection into a document:
one row per variable, the unweighted count of every category, and a
header with the title, authors, and notes. It prints the list of
variables and returns the codebook, which
output = "<path>.xlsx" writes to Excel and
output = "<path>.pdf" compiles to a PDF (this needs
the quarto package and Quarto 1.7 or later).
code_book(sochealth, sex, age, smoking, bmi)
#> Codebook
#>
#> Date: 2026-10-24
#> Observations: 1200
#> Variables: 4
#> Note: Counts and percentages are unweighted: they describe the data file and are not estimates for a population.
#>
#> Pos. │ Variable Label Type Valid Missing
#> ────────┼────────────────────────────────────────────────────────────────────────────
#> 1 │ sex Sex categorical (nominal) 1200 0
#> 2 │ age Age (years) numeric 1200 0
#> 10 │ smoking Current smoker categorical (nominal) 1175 25
#> 20 │ bmi Body mass index numeric 1188 12freq() produces frequency tables with counts,
percentages, valid percentages (shown by default whenever missing values
are present), and optionally cumulative percentages.
freq(sochealth, education)
#> Frequency table: education
#>
#> Category │ Values Freq. Percent
#> ────────────┼───────────────────────────────────────
#> Valid │ Lower secondary 261 21.8
#> │ Upper secondary 539 44.9
#> │ Tertiary 400 33.3
#> ────────────┼───────────────────────────────────────
#> Total │ 1200 100.0
#>
#> Label: Highest education level
#> Class: ordered, factor
#> Data: sochealthWeighted frequencies use the weights argument. With
rescale = TRUE, the total weighted N matches the unweighted
N:
freq(sochealth, education, weights = weight, rescale = TRUE)
#> Frequency table: education
#>
#> Category │ Values Freq. Percent
#> ────────────┼───────────────────────────────────────
#> Valid │ Lower secondary 259 21.6
#> │ Upper secondary 546 45.5
#> │ Tertiary 395 32.9
#> ────────────┼───────────────────────────────────────
#> Total │ 1200 100.0
#>
#> Label: Highest education level
#> Class: ordered, factor
#> Data: sochealth
#> Weight: weight (rescaled)cross_tab() crosses two categorical variables. By
default it shows counts, a chi-squared test, and Cramer’s V:
cross_tab(sochealth, smoking, education)
#> Crosstable: smoking x education (N)
#>
#> Values │ Lower secondary Upper secondary Tertiary │ Total
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> No │ 179 415 332 │ 926
#> Yes │ 78 112 59 │ 249
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> Total │ 257 527 391 │ 1175
#>
#> Chi-2(2) = 21.6, p < .001
#> Cramer's V = 0.14
#> Missing values removed: smoking (25).Add percentages with percent:
cross_tab(sochealth, smoking, education, percent = "column")
#> Crosstable: smoking x education (Column %)
#>
#> Values │ Lower secondary Upper secondary Tertiary │ Total
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> No │ 69.6 78.7 84.9 │ 78.8
#> Yes │ 30.4 21.3 15.1 │ 21.2
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> Total │ 100.0 100.0 100.0 │ 100.0
#> N │ 257 527 391 │ 1175
#>
#> Chi-2(2) = 21.6, p < .001
#> Cramer's V = 0.14
#> Missing values removed: smoking (25).Group by a third variable with by:
cross_tab(sochealth, smoking, education, by = sex)
#> Crosstable: smoking x education (N) | sex = Female
#>
#> Values │ Lower secondary Upper secondary Tertiary │ Total
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> No │ 95 220 160 │ 475
#> Yes │ 38 62 31 │ 131
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> Total │ 133 282 191 │ 606
#>
#> Chi-2(2) = 7.1, p = .029
#> Cramer's V = 0.11
#> Missing values removed: smoking (14).
#>
#> Crosstable: smoking x education (N) | sex = Male
#>
#> Values │ Lower secondary Upper secondary Tertiary │ Total
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> No │ 84 195 172 │ 451
#> Yes │ 40 50 28 │ 118
#> ──────────┼──────────────────────────────────────────────────┼─────────
#> Total │ 124 245 200 │ 569
#>
#> Chi-2(2) = 15.6, p < .001
#> Cramer's V = 0.17
#> Missing values removed: smoking (11).When both variables are ordered factors, cross_tab()
automatically selects an ordinal measure (Kendall’s Tau-b) instead of
Cramer’s V:
cross_tab(sochealth, self_rated_health, education)
#> Crosstable: self_rated_health x education (N)
#>
#> Values │ Lower secondary Upper secondary Tertiary │ Total
#> ─────────────┼──────────────────────────────────────────────────┼─────────
#> Poor │ 28 28 5 │ 61
#> Fair │ 86 118 62 │ 266
#> Good │ 102 263 193 │ 558
#> Very good │ 44 118 133 │ 295
#> ─────────────┼──────────────────────────────────────────────────┼─────────
#> Total │ 260 527 393 │ 1180
#>
#> Chi-2(6) = 73.2, p < .001
#> Kendall's Tau-b = 0.20
#> Missing values removed: self_rated_health (20).For a quick overview of all available association statistics, pass a
contingency table to assoc_measures():
tbl <- xtabs(~ smoking + education, data = sochealth)
assoc_measures(tbl)
#> Measure Estimate SE CI lower CI upper p
#> Cramer's V 0.136 – 0.079 0.191 <.001
#> Contingency Coefficient 0.134 – – – <.001
#> Lambda symmetric 0.000 0.000 0.000 0.000 –
#> Lambda R|C 0.000 0.000 0.000 0.000 –
#> Lambda C|R 0.000 0.000 0.000 0.000 –
#> Goodman-Kruskal's Tau R|C 0.018 0.008 0.003 0.034 .023
#> Goodman-Kruskal's Tau C|R 0.008 0.003 0.001 0.014 .022
#> Uncertainty Coefficient symmetric 0.011 0.005 0.002 0.021 .021
#> Uncertainty Coefficient R|C 0.018 0.008 0.003 0.032 .021
#> Uncertainty Coefficient C|R 0.009 0.004 0.001 0.016 .021
#> Goodman-Kruskal Gamma -0.268 0.056 -0.378 -0.158 <.001
#> Kendall's Tau-b -0.126 0.027 -0.180 -0.073 <.001
#> Stuart's Tau-c -0.117 0.026 -0.167 -0.067 <.001
#> Somers' D R|C -0.091 0.020 -0.131 -0.052 <.001
#> Somers' D C|R -0.175 0.038 -0.249 -0.101 <.001Individual functions such as cramer_v(),
gamma_gk(), or kendall_tau_b() return a scalar
by default. Pass detail = TRUE for the confidence interval
and p-value:
table_categorical() covers grouped or one-way summary
tables for categorical variables:
table_categorical(
sochealth,
select = c(smoking, physical_activity, dentist_12m),
by = education
)
#> Categorical table by education
#>
#> Variable │ Lower secondary n Lower secondary % Upper secondary n
#> ───────────────────┼─────────────────────────────────────────────────────────
#> smoking │
#> No │ 179 68.6 415
#> Yes │ 78 29.9 112
#> (Missing) │ 4 1.5 12
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> physical_activity │
#> No │ 177 67.8 310
#> Yes │ 84 32.2 229
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> dentist_12m │
#> No │ 113 43.3 174
#> Yes │ 148 56.7 365
#>
#> Variable │ Upper secondary % Tertiary n Tertiary % Total n
#> ───────────────────┼────────────────────────────────────────────────────
#> smoking │
#> No │ 77.0 332 83.0 926
#> Yes │ 20.8 59 14.8 249
#> (Missing) │ 2.2 9 2.2 25
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> physical_activity │
#> No │ 57.5 163 40.8 650
#> Yes │ 42.5 237 59.2 550
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> dentist_12m │
#> No │ 32.3 67 16.8 354
#> Yes │ 67.7 333 83.2 846
#>
#> Variable │ Total % p V
#> ───────────────────┼─────────────────────
#> smoking │ <.001 .14
#> No │ 77.2
#> Yes │ 20.8
#> (Missing) │ 2.1
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> physical_activity │ <.001 .21
#> No │ 54.2
#> Yes │ 45.8
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> dentist_12m │ <.001 .22
#> No │ 29.5
#> Yes │ 70.5
#>
#> Note. P-values: Pearson chi-squared test. V = Cramer's V.table_continuous() summarizes continuous variables,
either overall or by a categorical by variable; when
by is supplied it also reports a group-comparison test (a
Welch test by default):
table_continuous(
sochealth,
select = c(bmi, life_sat_health),
by = education
)
#> Descriptive statistics by education
#>
#> Variable │ Group M SD Min Max 95% CI LL
#> ─────────────────┼───────────────────────────────────────────────────────
#> bmi │ Lower secondary 28.09 3.47 18.20 38.90 27.66
#> │ Upper secondary 26.02 3.43 16.00 37.10 25.73
#> │ Tertiary 24.39 3.52 16.00 33.00 24.04
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> life_sat_health │ Lower secondary 2.71 1.20 1.00 5.00 2.57
#> │ Upper secondary 3.53 1.19 1.00 5.00 3.43
#> │ Tertiary 4.11 1.04 1.00 5.00 4.01
#>
#> Variable │ Group 95% CI UL n p
#> ─────────────────┼────────────────────────────────────────
#> bmi │ Lower secondary 28.51 260 <.001
#> │ Upper secondary 26.31 534
#> │ Tertiary 24.74 394
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> life_sat_health │ Lower secondary 2.86 259 <.001
#> │ Upper secondary 3.63 534
#> │ Tertiary 4.21 399
#>
#> Note. Missing values removed: bmi (12), life_sat_health (8). P-values: Welch one-way ANOVA. CI = confidence interval; LL = lower limit; UL = upper limit.table_continuous_lm() covers the same reporting
territory when you want to stay in a linear-model framework, for example
with robust or cluster-robust standard errors, case weights, or additive
covariate adjustment:
table_continuous_lm(
sochealth,
select = c(wellbeing_score, bmi),
by = sex,
vcov = "HC3"
)
#> Continuous outcomes by sex
#>
#> M
#> ─────────────
#> Variable │ Female Male Δ 95% CI LL 95% CI UL p R² n
#> ─────────────────┼──────────────────────────────────────────────────────────────
#> wellbeing_score │ 67.16 71.05 3.89 2.12 5.65 <.001 0.02 1200
#> bmi │ 25.69 26.20 0.51 0.09 0.93 .018 0.00 1188
#>
#> Note. Missing values removed: bmi (12). P-values: t-test from the linear model. Std. errors: heteroskedasticity-robust (HC3). Δ = Male − Female. CI = confidence interval; LL = lower limit; UL = upper limit.By default every table_* helper formats to APA
conventions. To match a journal instead, pass style = —
"jama", "nejm", "lancet",
"annals", "apa", "aer";
spicy_style_names() lists them all. Set it once for a whole
document with options(spicy.style = "nejm"), or compose a
custom variant with spicy_style(), for example
spicy_style("lancet", ci_sep = " to "). The same styles
apply to the regression tables below. French typographic conventions are
not a journal style: they come with the language, through
options(spicy.language = "fr").
For detailed guidance, see the dedicated articles on
table_categorical(), table_continuous(), and
table_continuous_lm(), and the final reporting overview
tying the summary tables together.
The items of a “check all that apply” question are answered by the
same respondents, and a plain categorical table misreads them. Stacked
into one variable, the choices outnumber the respondents, and the
chi-squared test of the items-by-groups table is not valid, since one
respondent contributes to several cells.
table_multiresponse() reads the items as one question: one
row per item, with the number of respondents who chose it and their
percentage, a row for those who chose none, then a row with the number
of respondents. In sochealth, the five logical
info_* items answer such a question: the sources of health
information used in the past 12 months. labels = TRUE shows
the item wording:
table_multiresponse(
sochealth,
select = starts_with("info_"),
labels = TRUE
)
#> Multiple responses
#>
#> Variable │ n %
#> ───────────────────────────────────────┼────────────────
#> Doctor or other health professional │ 721 61.9
#> Pharmacist │ 256 22.0
#> Internet, websites or apps │ 523 44.9
#> Family or friends │ 360 30.9
#> Television, radio or newspapers │ 141 12.1
#> (No item chosen) │ 67 5.8
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Respondents │ 1164
#>
#> Note. Multiple answers possible: percentages can sum to more than 100%. Base: respondents with data for the question (excluded: 36 with no data for the question).With by, the counts and percentages are given per group
and in a Total margin, each item gets a chi-squared test
(Fisher’s exact test when an expected count is below 5), and the
question a global test that allows for several answers per respondent,
both explained in the Multiple-response
questions article. p_value = FALSE leaves a descriptive
table:
table_multiresponse(
sochealth,
select = starts_with("info_"),
by = sex,
p_value = FALSE
)
#> Multiple responses by sex
#>
#> Variable │ Female n Female % Male n Male % Total n Total %
#> ──────────────────┼──────────────────────────────────────────────────────
#> info_doctor │ 377 62.9 344 60.9 721 61.9
#> info_pharmacist │ 137 22.9 119 21.1 256 22.0
#> info_internet │ 268 44.7 255 45.1 523 44.9
#> info_family │ 219 36.6 141 25.0 360 30.9
#> info_media │ 78 13.0 63 11.2 141 12.1
#> (No item chosen) │ 31 5.2 36 6.4 67 5.8
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> Respondents │ 599 565 1164
#>
#> Note. Multiple answers possible: percentages can sum to more than 100%. Base: respondents with data for the question (excluded: 36 with no data for the question).table_regression() reports the full coefficient table
for one or several fitted models — 39 classes, from lm() /
glm() to mixed-effects, ordinal, survival and Bayesian
engines (the full map is the Supported
models article) — with APA formatting by default and the same
journal styles on demand (style = "jama",
"nejm", "lancet" and more), factor grouping
with reference rows, robust variance, standardized coefficients, average
marginal effects, hierarchical comparisons, and side-by-side multi-model
layouts:
fit <- lm(wellbeing_score ~ age + sex + smoking, data = sochealth)
table_regression(fit)
#> Linear regression: wellbeing_score
#>
#> Variable │ B SE 95% CI p
#> ─────────────────┼──────────────────────────────────────
#> (Intercept) │ 65.20 1.66 [61.95, 68.45] <.001
#> age │ 0.05 0.03 [-0.01, 0.11] .130
#> sex: │
#> Female (ref.) │ – – – –
#> Male │ 3.86 0.91 [ 2.08, 5.63] <.001
#> smoking: │
#> No (ref.) │ – – – –
#> Yes │ -1.72 1.11 [-3.89, 0.45] .121
#> ╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌┼╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌
#> n │ 1175
#> R² │ 0.02
#> Adj. R² │ 0.02
#>
#> Note. Linear regression. Std. errors: classical (OLS). CI = confidence interval.For detailed guidance, start with the core
table_regression() article; the Learn more section below
lists the companion articles by model family.
mean_n(), sum_n(), and
count_n() compute row-wise statistics across selected
columns, with explicit control over missing values. By default, a row
mean or sum requires every selected value to be valid;
min_valid relaxes this to a minimum count (or proportion)
of valid values. The rows selected below include respondents with
missing satisfaction items so the difference is visible:
mean_sat (no min_valid) turns NA
as soon as one item is missing, sum_sat is computed from
the valid items as long as at least three are present (and turns
NA in the last row, which has only two), and
count_n(special = "NA") counts the missing items
themselves:
sochealth |>
dplyr::mutate(
mean_sat = mean_n(select = starts_with("life_sat")),
sum_sat = sum_n(select = starts_with("life_sat"), min_valid = 3),
n_missing = count_n(select = starts_with("life_sat"), special = "NA")
) |>
dplyr::select(starts_with("life_sat"), mean_sat, sum_sat, n_missing) |>
dplyr::slice(c(1, 2, 43, 82, 455)) |>
as.data.frame()
#> life_sat_health life_sat_work life_sat_relationships life_sat_standard
#> 1 5 3 5 5
#> 2 4 4 5 5
#> 3 2 NA 3 4
#> 4 1 NA 2 2
#> 5 5 NA 4 NA
#> mean_sat sum_sat n_missing
#> 1 4.5 18 0
#> 2 4.5 18 0
#> 3 NA 9 1
#> 4 NA 5 1
#> 5 NA NA 2This is the vignette the package ships. Every walk-through below is an article on the spicy website, kept alongside the reference pages and rebuilt with each release.
Explore and tabulate
varlist(), code_book(),
label_from_names()).freq() and
cross_tab() in depth (weights, simulation, labelled
data).Summary tables
style =).survey design: design-based standard errors, degrees of
freedom and tests (table_continuous_svy(),
table_categorical_svy()).Regression tables
table_regression_uv()).Putting it together
table_*() functions together.?mean_n, ?sum_n, ?count_n —
row-wise summaries with minimum-valid-values rules.