Constant terms a and b were calculated by regressing the predicted number of fixations onto log2(N-1).
Here, we simulate single attribute decisions with threshold increase 𝜃Δ = .01, search sensitivity α = 5 and sampling noise 𝜎 = 1, although MASC consistently predicts this logarithmic relationship across a range of parameter sets.
# masc_hicks_law.R - Creates Figure 11 reproduction from
# Gluth, S., Deakin, J., & Rieskamp, J. (2026). A theory of multiattribute search and choice. Psychological Review.
library(ggplot2)
library(dplyr)
library(tibble)
library(purrr)
library(patchwork)
library(masc)
# Function to simulate relationship between number of options and fixations
simulate_hicks_law <- function(
n_trials = 200, # Number of trials per condition
n_subjects = 100, # Number of simulated participants
max_options = 20, # Maximum number of options to test
min_options = 2, # Minimum number of options to test
n_attributes = 1, # Use single attribute decisions
sampling_noise = 1, # Fixed sampling noise
alpha = 5, # Fixed search sensitivity
delta = 0.01, # Fixed threshold increase
theta = 0.01, # Fixed initial threshold
include_approx_rt = FALSE, # Whether to include approx RT calculation
plot_while_running = FALSE # Whether to update plot as simulation runs
) {
# Pre-allocate results
options_to_test <- min_options:max_options
mean_fixations <- numeric(length(options_to_test))
mean_approx_rt <- numeric(length(options_to_test))
# ExGaussian parameters (from paper)
mu <- 139.5530
sigma <- 63.1433
tau <- 91.0133
# Create color
start_color <- c(194, 218, 184) / 255
# Create and initialize plot if requested
if (plot_while_running) {
plot(1, 1, type = "n",
xlim = c(min_options-1, max_options+1),
ylim = c(0, 100), # Will adjust as data comes in
xlab = "Number of Options (N)",
ylab = "Number of Fixations",
main = "Hick's Law in Multi-Alternative Decisions")
}
# Generate weights
weights <- rep(1, n_attributes) # Equal weights for single attribute
# Loop through different numbers of options
for (i in seq_along(options_to_test)) {
n_opts <- options_to_test[i]
cat(sprintf("Processing %d options (%d/%d)\n",
n_opts, i, length(options_to_test)))
# Pre-allocate result arrays
all_fixations <- matrix(0, nrow = n_trials, ncol = n_subjects)
all_approx_rt <- matrix(0, nrow = n_trials, ncol = n_subjects)
# Process each subject
for (s in 1:n_subjects) {
# Generate attribute values for all trials
trial_data <- do.call(rbind, lapply(1:n_trials, function(t) {
# Generate valid attribute values
# For single attribute, we don't need to check for dominance
attr_vals <- matrix(rnorm(n_opts * n_attributes), nrow = n_opts)
as.data.frame(matrix(attr_vals, nrow = 1))
}))
# Rename columns appropriately for rMASC
colnames(trial_data) <- c(
paste0("opt", rep(1:n_opts, each = n_attributes),
"_att", rep(1:n_attributes, n_opts))
)
# Run rMASC for this subject with current parameters
result <- rMASC(
data = trial_data,
w = weights,
sigma = sampling_noise,
alpha = alpha,
delta = delta,
theta = theta,
n_options = n_opts,
n_attributes = n_attributes,
max_steps = 100
)
# Extract fixation counts
all_fixations[, s] <- result$results$rt
# Calculate approximate RT if requested
if (include_approx_rt) {
for (t in 1:n_trials) {
curr_rt <- result$results$rt[t]
gauss_comp <- rnorm(curr_rt, mu, sigma)
exp_comp <- rexp(curr_rt, 1/tau) # R uses rate parameter (1/scale)
all_approx_rt[t, s] <- sum(gauss_comp + exp_comp)
}
}
}
# Calculate means
mean_fixations[i] <- mean(all_fixations)
if (include_approx_rt) {
mean_approx_rt[i] <- mean(all_approx_rt)
}
# Update plot if requested
if (plot_while_running) {
# Fit Hick's Law model
x_vals <- log2(options_to_test[1:i] - 1)
y_vals <- mean_fixations[1:i]
model_fit <- lm(y_vals ~ x_vals)
a <- model_fit$coefficients[1]
b <- model_fit$coefficients[2]
# Generate prediction
x_pred <- 1:max_options
y_pred <- a + b * log2(pmax(x_pred - 1, 1)) # Avoid log(0)
# Clear plot
plot(options_to_test[1:i], mean_fixations[1:i],
pch = 21, bg = "gray50", col = "transparent",
xlim = c(min_options-1, max_options+1),
ylim = c(0, max(y_pred, mean_fixations[1:i], na.rm = TRUE) * 1.1),
xlab = "Number of Options (N)",
ylab = "Number of Fixations",
main = "Hick's Law in Multi-Alternative Decisions")
# Add prediction line
lines(x_pred, y_pred, col = rgb(start_color[1], start_color[2], start_color[3]), lwd = 2)
# Add legend
legend("topleft",
legend = c("Number of Fixations", "Hick's Law: a + b * log₂(N-1)"),
pch = c(21, NA),
lty = c(NA, 1),
pt.bg = c("gray50", NA),
col = c("transparent", rgb(start_color[1], start_color[2], start_color[3])),
lwd = c(NA, 2))
# Force update
Sys.sleep(0.1)
}
}
# Prepare results
results <- tibble(
n_options = options_to_test,
mean_fixations = mean_fixations
)
if (include_approx_rt) {
results$mean_approx_rt <- mean_approx_rt
}
return(results)
}
# Function to plot Hick's Law results
plot_hicks_law <- function(results, include_approx_rt = FALSE) {
# Fit Hick's Law model to fixation data
x_vals <- log2(results$n_options - 1)
y_vals <- results$mean_fixations
model_fit <- lm(y_vals ~ x_vals)
a <- model_fit$coefficients[1]
b <- model_fit$coefficients[2]
# Generate prediction
x_pred <- 1:max(results$n_options)
y_pred <- a + b * log2(pmax(x_pred - 1, 1)) # Avoid log(0)
# Start color from paper
start_color <- c(194, 218, 184) / 255
# Create plot data
plot_data <- tibble(
n_options = results$n_options,
mean_fixations = results$mean_fixations
)
pred_data <- tibble(
n_options = x_pred,
predicted = y_pred
)
# Create fixation plot
p1 <- ggplot() +
geom_point(data = plot_data,
aes(x = n_options, y = mean_fixations),
size = 3, shape = 21, fill = "gray50", color = "black") +
geom_line(data = pred_data,
aes(x = n_options, y = predicted),
color = rgb(start_color[1], start_color[2], start_color[3]),
size = 1.2) +
labs(
title = "Relationship between Number of Fixations and Number of Options",
x = "Number of Options (N)",
y = "Number of Fixations"
) +
theme_classic() +
theme(
plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
# Create RT plot if requested
if (include_approx_rt && "mean_approx_rt" %in% names(results)) {
# Fit Hick's Law model to RT data
rt_model <- lm(results$mean_approx_rt ~ x_vals)
rt_a <- rt_model$coefficients[1]
rt_b <- rt_model$coefficients[2]
# Generate RT prediction
rt_pred <- rt_a + rt_b * log2(pmax(x_pred - 1, 1))
# Create RT plot data
rt_plot_data <- tibble(
n_options = results$n_options,
mean_rt = results$mean_approx_rt
)
rt_pred_data <- tibble(
n_options = x_pred,
predicted = rt_pred
)
p2 <- ggplot() +
geom_point(data = rt_plot_data,
aes(x = n_options, y = mean_rt),
size = 3, shape = 21, fill = "gray50", color = "black") +
geom_line(data = rt_pred_data,
aes(x = n_options, y = predicted),
color = rgb(start_color[1], start_color[2], start_color[3]),
size = 1.2) +
labs(
title = "Relationship between Reaction Time and Number of Options",
x = "Number of Options (N)",
y = "Reaction Time (ms)"
) +
theme_classic() +
theme(
plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
# Combine plots
return(p1 + p2 + plot_layout(ncol = 2))
}
# Otherwise return just the fixation plot
return(p1)
}# Run simulation (note: this can take a while)
set.seed(2025)
results <- simulate_hicks_law(
n_trials = 15,
n_subjects = 6,
max_options = 8,
min_options = 2,
include_approx_rt = FALSE,
plot_while_running = FALSE
)
#> Processing 2 options (1/7)
#> Processing 3 options (2/7)
#> Processing 4 options (3/7)
#> Processing 5 options (4/7)
#> Processing 6 options (5/7)
#> Processing 7 options (6/7)
#> Processing 8 options (7/7)# Plot results
plot <- plot_hicks_law(results)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
print(plot)
# Calculate and report model fit
x_vals <- log2(results$n_options - 1)
model_fit <- lm(results$mean_fixations ~ x_vals)
summary(model_fit)
#>
#> Call:
#> lm(formula = results$mean_fixations ~ x_vals)
#>
#> Residuals:
#> 1 2 3 4 5 6 7
#> -0.5671 1.1500 -0.7250 0.4893 -0.2507 0.5698 -0.6662
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 11.5671 0.6590 17.554 1.1e-05 ***
#> x_vals 2.1607 0.3327 6.495 0.00129 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8051 on 5 degrees of freedom
#> Multiple R-squared: 0.894, Adjusted R-squared: 0.8728
#> F-statistic: 42.19 on 1 and 5 DF, p-value: 0.001291
cat(sprintf("\nHick's Law equation: Number of fixations = %.2f + %.2f * log₂(N-1)\n",
model_fit$coefficients[1], model_fit$coefficients[2]))
#>
#> Hick's Law equation: Number of fixations = 11.57 + 2.16 * log₂(N-1)
cat(sprintf("R² = %.4f\n", summary(model_fit)$r.squared))
#> R² = 0.8940