The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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.8940These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.