Using guess

guess: Adjust Estimates of Learning for Guessing

Over informative processes, naive estimator of learning—difference between post and pre process scores—underestimates actual learning. A heuristic account for why the naive estimator is negatively biased is as follows: people know as much or more after exposed to an informative process than before it. And the less people know, the larger the number of items they don’t know. And greater the opportunity to guess.

Guessing, even when random, only increases the proportion correct. Thus, bias due to guessing for naive measures of knowledge is always positive. On average, thus, there is more positive bias in the pre-process scores than post-process scores. And naturally, subtracting pre-process scores from post-process provides an attenuated estimate of actual learning. For a more complete treatment of the issue, read this paper by Ken Cor and Gaurav Sood.

We provide a few different ways to adjust estimates of learning for guessing. For now, we limit our attention to cases where the same battery of knowledge questions has been asked in both the pre- and the post-process wave. And to cases where closed-ended questions have been asked. (Guessing is not a serious issue on open-ended items. See more evidence for that in DK Means DK by Robert Luskin and John Bullock.) More generally, the package implements the methods to adjust learning for guessing discussed in this paper.

Measuring Learning:

Estimand

Proportion of people who learned a particular piece of information over the course of an informative process.

Other Issues

Measurement of knowledge is fundamentally reactive – we must probe to learn. But probing is not without its problems. For instance, people who don’t know the answer try to triangulate based on the cues in the question itself. For another, people are remarkably averse to confessing to their ignorance. So on a closed ended question, lots of people who don’t know the right answer, guess. Here are some pertinent issues that relate to how we analyze the data:

  1. Dealing with Missing Data
    If you assume that the data are missing completely at random, you can simply ignore them. Generally, however, respondents tend to skip items they don’t know. So missing responses on knowledge questions typically indicate ignorance. (Of course, it is important to investigate other potential reasons behind missing data. And we encourage researchers to take all precautions.) In our treatment, however, for simplicity sake, we treat missing as indicators of ignorance.

  2. Dealing with Don’t Knows
    We now know a little bit about Don’t Know. One generally strategy is to treat Don’t Know responses as ignorance. But research suggests that on average there is approximately 3% hidden knowledge behind Don’t Know responses. See DK Means DK by Robert Luskin and John Bullock. Thus one can also choose to replace Don’t Know responses with .03.

  3. Related Knowledge
    People either know a particular piece of information or they don’t. On an open-ended question, they may struggle to remember it but those kinds of concerns don’t apply to closed-ended questions where the task is simply to identify the correct answer. What does on occassion happen on closed-ended questions is that people have related cognitions. So for instance, asked to identify the prime minister of France, the respondents sometimes know that one of the options is the king of Sudan and may guess randomly between the remaining options. But that isn’t the same as knowing the prime minister of France.

Standard Correction for Guessing

The standard correction for guessing assumes that people guess randomly. And that people either know or don’t know. Using this assumption, it then uses total number of incorrect answers to estimate the total number of items that the person guessed on. For instance, let us assume there are 4 options on a multiple choice question. Say we have data from 100 respondents. Say there are 70 incorrect answers and 30 correct. Incorrect answers reflect attempts of guessing. (We also assume that people aren’t misinformed.) This means we can triangulat the total number of questions respondents guessed on – 70*(4/3). This means that the proportion of people who know the piece of information is roughly .067. Do it for the pre and the post wave and you have estimate of learning adjusted for guessing using the standard correction.

Latent Class Correction for Guessing

See the paper for details.

Installation

To get the current CRAN version:

install.packages("guess")

To get the current development version from GitHub:

# install.packages("devtools") # nolint
library(devtools)
# devtools::install_github("soodoku/guess") # nolint

Usage

Standard Correction for Guessing

To adjust estimates of learning for standard correction of guessing, use stnd_cor. The function requires takes pre test and post test data frames containing responses to the items on the pre- and the post-test, and a lucky vector that contains the probability of getting an item correct when guessing randomly. Under standard guessing correction, it is taken to be inverse of total number of options.

Structure of the Input Data:

  1. For current purposes, we assume missing responses to indicate ignorance. Thus functions internally code missing responses as 0.
  2. If the items offer an option to mark Don't know, code all Don't Know responses as ‘d’.
# Load library
library(guess)

# Generate some data without DK
pre_test <- data.frame(item1 = c(1, 0, 0, 1, 0), item2 = c(1, NA, 0, 1, 0))
pst_test <- pre_test + cbind(c(0, 1, 1, 0, 0), c(0, 1, 0, 0, 1))
lucky <- rep(.25, 2)

# Unadjusted Effect
# Treating missing responses as ignorance
pre_ignorance <- replace(pre_test, is.na(pre_test), 0)
pst_ignorance <- replace(pst_test, is.na(pst_test), 0)
colMeans(pst_ignorance - pre_ignorance)
## item1 item2 
##   0.4   0.2
# MCAR
colMeans(pst_test - pre_test, na.rm = TRUE)
## item1 item2 
##  0.40  0.25
# Adjusted Effect
stnd_cor(pre_test, pst_test, lucky)
## $pre
##     item1     item2 
## 0.2000000 0.2666667 
## 
## $pst
##     item1     item2 
## 0.7333333 0.5333333 
## 
## $learn
##     item1     item2 
## 0.5333333 0.2666667

Transition Matrix

# Without Don't Know
pre_test_var <- c(1, 0, 0, 1, 0, 1, 0)
pst_test_var <- c(1, 0, 1, 1, 0, 1, 1)
count_transitions(pre_test_var, pst_test_var)
## x00 x01 x10 x11 
##   2   2   0   3
# With Don't Know
pre_test_var <- c(1, 0, NA, 1, "d", "d", 0, 1, 0)
pst_test_var <- c(1, 0, 1, "d", 1, 0, 1, 1, "d")
count_transitions(pre_test_var, pst_test_var)
## x00 x01 x0d x10 x11 x1d xd0 xd1 xdd 
##   1   1   1   0   2   1   1   2   0

Adjusting Using the Latent Class Model

# Create example data for demonstration
set.seed(123)
nitems <- 10 # Number of knowledge questions
npeople <- 100 # Number of respondents

# Generate simulated pre-test and post-test responses (0=incorrect, 1=correct)
pre_test <- as.data.frame(replicate(nitems, rbinom(npeople, 1, 0.4)))
post_test <- as.data.frame(replicate(nitems, rbinom(npeople, 1, 0.6)))
names(pre_test) <- names(post_test) <- paste0("item", seq_len(nitems))

transmatrix <- count_item_transitions(pre_test, post_test)

res <- fit_item_lca_counts(transmatrix)

round(res$params[, 1:4], 3)
##       item1 item2 item3 item4
## gg     0.55 0.570 0.602 0.525
## gk     0.45 0.385 0.364 0.475
## kk     0.00 0.044 0.034 0.000
## gamma  0.40 0.351 0.368 0.390
round(res$learning[1:4], 3)
## item1 item2 item3 item4 
## 0.450 0.385 0.364 0.475
# LCA Correction with Don't Know responses
# Generate data with Don't Know responses (coded as "d")
pre_test_dk <- as.data.frame(replicate(
  nitems, sample(c("0", "1", "d"), npeople, replace = TRUE, prob = c(0.5, 0.35, 0.15))
))
post_test_dk <- as.data.frame(replicate(
  nitems, sample(c("0", "1", "d"), npeople, replace = TRUE, prob = c(0.4, 0.45, 0.15))
))
names(pre_test_dk) <- names(post_test_dk) <- paste0("item", seq_len(nitems))

transmatrix <- count_item_transitions(pre_test_dk, post_test_dk)
res_dk <- fit_item_lca_counts(transmatrix)

round(res_dk$params[, 1:4], 3)
##       item1 item2 item3 item4
## gg    0.672 0.514 0.410 0.692
## gk    0.038 0.174 0.236 0.048
## gd    0.130 0.090 0.130 0.100
## kk    0.000 0.081 0.014 0.000
## dg    0.050 0.110 0.117 0.099
## dk    0.090 0.000 0.063 0.031
## dd    0.020 0.030 0.030 0.030
## gamma 0.405 0.358 0.317 0.393
round(res_dk$learning[1:4], 3)
## item1 item2 item3 item4 
## 0.128 0.174 0.299 0.079

Adjust by Groups

Account for propensity to guess and item level attribute that a guess would be lucky.

pre_test_var <- data.frame(
  pre = c(1, 0, 1, 1, 0, "d", "d", 0, 1, NA, 0, 1, 1, 1, 1, 0, 0, "d", 0, 0)
)
pst_test_var <- data.frame(pst = c(1, 0, NA, 1, "d", 1, 0, 1, 1, "d", 1, 1, 1, 0, 1, 1, 0, 1, 1, 0))
grp <- c(rep(1, 10), rep(0, 10))

group_adj(
  pre_test_var, pst_test_var,
  guessing_probability = 0, knowledge_given_dont_know = 0
)$mean_learning
## pst 
## 0.2
group_adj(
  pre_test_var, pst_test_var,
  guessing_probability = .25, knowledge_given_dont_know = 0
)$mean_learning
##  pst 
## 0.25
stnd_cor(pre_test_var, pst_test_var, guessing_probability = .25)$learn
##  pre 
## 0.25
grp0_raw <- group_adj(
  subset(pre_test_var, grp == 0), subset(pst_test_var, grp == 0),
  guessing_probability = 0, knowledge_given_dont_know = 0
)$mean_learning
grp1_raw <- group_adj(
  subset(pre_test_var, grp == 1), subset(pst_test_var, grp == 1),
  guessing_probability = 0, knowledge_given_dont_know = 0
)$mean_learning

grp0_adj <- group_adj(
  subset(pre_test_var, grp == 0), subset(pst_test_var, grp == 0),
  guessing_probability = .25, knowledge_given_dont_know = 0
)$mean_learning
grp1_adj <- group_adj(
  subset(pre_test_var, grp == 1), subset(pst_test_var, grp == 1),
  guessing_probability = .25, knowledge_given_dont_know = 0
)$mean_learning

grp0_raw - grp1_raw
## pst 
## 0.2
grp0_adj - grp1_adj
##       pst 
## 0.2333333

Standard Errors

# Raw: standard error of the mean pre/post difference
change <- post_test - pre_test
sapply(change, function(x) sd(x) / sqrt(length(x)))
##      item1      item2      item3      item4      item5      item6      item7 
## 0.06942040 0.06723245 0.06794977 0.07005770 0.06614187 0.07086051 0.06256310 
##      item8      item9     item10 
## 0.07109533 0.06624107 0.06337319
# LCA model, bootstrapped over respondents
lca_stnd_err <- lca_se(pre_test, post_test, n_resamples = 10, seed = 123)
round(lca_stnd_err$learning_standard_error, 3)
##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
## lca 0.118 0.103 0.074 0.076 0.099 0.116 0.056 0.109 0.098 0.082 0.041
lca_dk_stnd_err <- lca_se(pre_test_dk, post_test_dk, n_resamples = 10, seed = 123)
round(lca_dk_stnd_err$learning_standard_error, 3)
##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8]  [,9] [,10] [,11]
## lca 0.061 0.101 0.057 0.058 0.095 0.048 0.08 0.059 0.081 0.088 0.033

The last column of learning_standard_error is the standard error of the aggregate estimate, which pools transitions across items.

Fit

fit <- assess_item_lca_fit(res, pre_test, post_test)

print(fit$statistics)
##        statistic df p_value n_observations
## item1         NA  0      NA            100
## item2         NA  0      NA            100
## item3         NA  0      NA            100
## item4         NA  0      NA            100
## item5         NA  0      NA            100
## item6         NA  0      NA            100
## item7         NA  0      NA            100
## item8         NA  0      NA            100
## item9         NA  0      NA            100
## item10        NA  0      NA            100
fit <- assess_item_lca_fit(res_dk, pre_test_dk, post_test_dk)

print(fit$statistics)
##         statistic df    p_value n_observations
## item1  0.56262030  1 0.45320641            100
## item2  1.69887946  1 0.19243458            100
## item3  0.65478212  1 0.41840812            100
## item4  1.57219278  1 0.20988885            100
## item5  0.13750940  1 0.71076967            100
## item6  0.54400369  1 0.46077819            100
## item7  0.02855475  1 0.86581119            100
## item8  4.24012099  1 0.03947935            100
## item9  3.85346374  1 0.04964337            100
## item10 2.77290922  1 0.09587177            100

Model Criticism: Perplexity and Cross-Validation

The package provides tools to evaluate model fit through perplexity and cross-validation. These help assess whether the LCA model adequately captures the data-generating process.

Perplexity measures how well the model predicts the data. Lower values indicate better fit. It can be interpreted as the effective number of equiprobable outcomes the model assigns to each observation.

Cross-validation provides held-out performance estimates, helping detect overfitting.

Two sets of functions are available: - *_items: Work with aggregated transition matrices, split by items - *_individuals: Work with individual-level data, split by individuals

# Generate data
set.seed(42)
nitems <- 10
npeople <- 200

pre_test <- as.data.frame(replicate(nitems, rbinom(npeople, 1, 0.4)))
post_test <- as.data.frame(
  pmin(as.matrix(pre_test) + replicate(nitems, rbinom(npeople, 1, 0.25)), 1)
)
names(pre_test) <- names(post_test) <- paste0("item", 1:nitems)

# Fit model using convenience wrapper
fit <- fit_item_lca(pre_test, post_test)

# Or equivalently:
# transmatrix <- count_item_transitions(pre_test, post_test) # nolint
# fit <- fit_item_lca_counts(transmatrix) # nolint
Scores from aggregated item data
transmatrix <- count_item_transitions(pre_test, post_test)
score <- score_item_lca(fit, transmatrix)

score$perplexity
## [1] 2.771211
score$item_scores
##        log_likelihood n_observations mean_log_likelihood perplexity
## item1       -195.1970            200          -0.9759851   2.653780
## item2       -209.8447            200          -1.0492234   2.855433
## item3       -204.1951            200          -1.0209756   2.775902
## item4       -204.9668            200          -1.0248340   2.786633
## item5       -203.1894            200          -1.0159470   2.761978
## item6       -209.0896            200          -1.0454479   2.844672
## item7       -208.6109            200          -1.0430543   2.837872
## item8       -206.6019            200          -1.0330096   2.809509
## item9       -209.0896            200          -1.0454479   2.844672
## item10      -187.7836            200          -0.9389178   2.557212
Scoring individual data
individual_score <- score_individual_lca(fit, pre_test, post_test)

individual_score$perplexity
## [1] 2.771211
hist(
  individual_score$individual_scores$perplexity,
  main = "Distribution of Individual Perplexity"
)

Cross-validation over individuals

Useful when you want to assess generalization to new respondents.

cv_result <- cv_individual_lca(pre_test, post_test, k = 5, seed = 123)

cv_result$perplexity # Held-out perplexity
## [1] 2.815828
cv_result$fold_results # Per-fold details
##   fold train_n test_n  train_ll   test_ll test_obs test_ll_per_obs
## 1    1     160     40 -1645.315 -396.2255      400      -0.9905637
## 2    2     160     40 -1624.361 -417.8041      400      -1.0445103
## 3    3     160     40 -1624.870 -417.2485      400      -1.0431211
## 4    4     160     40 -1615.723 -427.8618      400      -1.0696546
## 5    5     160     40 -1629.882 -411.3726      400      -1.0284314

Simulation and Parameter Recovery Validation

The package provides functions to simulate data from known parameters, useful for validating the estimation procedure and understanding model behavior.

Basic simulation (No-DK model)
# Simulate data with known learning rate of 30%
sim <- simulate_lca(
  n = 500, # 500 individuals
  n_items = 2, # 2 test items
  gg = 0.35, # 35% in guess->guess (stable ignorance)
  gk = 0.30, # 30% in guess->know (LEARNED)
  kk = 0.35, # 35% in know->know (stable knowledge)
  gamma = 0.25, # 25% chance of guessing correctly
  seed = 123
)

# Fit model and check recovery
fit <- fit_item_lca(sim$pre, sim$post)
fit$params["gk", ] # Should be close to 0.30 (true learning rate)
##     item1     item2 
## 0.2746930 0.2901333
Simulation with Don’t Know responses
# Simulate DK data with known parameters
sim_dk <- simulate_lca_dk(
  n = 500,
  n_items = 1,
  gg = 0.25, gk = 0.15, gd = 0.10,
  kk = 0.20, dg = 0.10, dk = 0.10,
  dd = 0.10, gamma = 0.25,
  seed = 456
)

# Fit DK model
dk_fit <- fit_item_lca(sim_dk$pre, sim_dk$post)
dk_fit$params["gk", ] # Should be close to 0.15
## [1] 0.1316666
Comprehensive parameter recovery validation

Use validate_recovery() for Monte Carlo studies of parameter recovery:

# Run 100 simulations to assess parameter recovery
results <- validate_recovery(
  true_params = c(gg = 0.35, gk = 0.30, kk = 0.35, gamma = 0.25),
  n = 500, # Sample size per simulation
  n_items = 2, # Number of items
  n_sims = 100, # Number of simulations
  seed = 789
)

# View results
# Bias, RMSE, and Monte Carlo SD for each parameter, averaged across items
print(results)
##   parameter true_value mean_estimate          bias       rmse estimate_sd
## 1        gg       0.35     0.3507703  0.0007703033 0.02348697  0.02359259
## 2        gk       0.30     0.2960969 -0.0039031214 0.02904593  0.02892749
## 3        kk       0.35     0.3531328  0.0031328181 0.02905283  0.02902894
## 4     gamma       0.25     0.2531744  0.0031744373 0.02595291  0.02588780
##          mcse
## 1 0.002359259
## 2 0.002892749
## 3 0.002902894
## 4 0.002588780

Statistical Efficiency: Sample Size and Parameter Recovery

As sample size increases, the item-level LCA produces more precise estimates. This section shows how precision changes with sample size, number of items, and the true parameters, using simulation.

We assess estimation quality using three metrics:

  • Bias: Mean estimate minus true value (should approach 0)
  • RMSE: Root mean squared error (reflects both bias and variance)
  • SD: Standard deviation of the estimates across simulations
Sample Size Effect

The most fundamental efficiency result is that RMSE decreases approximately with 1/√n. Here we compare parameter recovery across different sample sizes:

# Compare parameter recovery across sample sizes
sample_sizes <- c(100, 250, 500, 1000)
true_params <- c(gg = 0.35, gk = 0.30, kk = 0.35, gamma = 0.25)

results_by_n <- lapply(sample_sizes, function(n) {
  validate_recovery(true_params, n = n, n_items = 2, n_sims = 100, seed = 123)
})

# Combine into summary table
efficiency_table <- do.call(rbind, lapply(seq_along(sample_sizes), function(i) {
  r <- results_by_n[[i]]
  data.frame(
    n = sample_sizes[i],
    parameter = r$parameter,
    bias = round(r$bias, 4),
    rmse = round(r$rmse, 4),
    sd = round(r$estimate_sd, 4)
  )
}))

# Display results for the learning parameter (gk)
gk_results <- efficiency_table[efficiency_table$parameter == "gk", ]
print(gk_results)
##       n parameter    bias   rmse     sd
## 2   100        gk  0.0017 0.0580 0.0583
## 6   250        gk -0.0038 0.0399 0.0399
## 10  500        gk -0.0009 0.0261 0.0262
## 14 1000        gk  0.0021 0.0184 0.0184

If RMSE falls with 1/√n, doubling the sample size divides it by about √2:

round(gk_results$rmse[gk_results$n == 500] / gk_results$rmse[gk_results$n == 1000], 2)
## [1] 1.42
Number of Items Effect

The item-level model fits each item separately, so adding items does not make any one item’s estimates more precise. validate_recovery() averages each parameter across items, however, and the simulated items share the same true values. The RMSE of that average falls as items are added:

# Compare recovery with varying number of items
n_items_values <- c(1, 2, 5, 10)
true_params <- c(gg = 0.35, gk = 0.30, kk = 0.35, gamma = 0.25)

results_by_items <- lapply(n_items_values, function(ni) {
  validate_recovery(true_params, n = 500, n_items = ni, n_sims = 50, seed = 456)
})

# Combine results
items_table <- do.call(rbind, lapply(seq_along(n_items_values), function(i) {
  r <- results_by_items[[i]]
  data.frame(
    n_items = n_items_values[i],
    parameter = r$parameter,
    rmse = round(r$rmse, 4),
    sd = round(r$estimate_sd, 4)
  )
}))

# Show gamma results (most affected by number of items)
gamma_results <- items_table[items_table$parameter == "gamma", ]
print(gamma_results)
##    n_items parameter   rmse     sd
## 4        1     gamma 0.0393 0.0393
## 8        2     gamma 0.0302 0.0304
## 12       5     gamma 0.0171 0.0173
## 16      10     gamma 0.0108 0.0109

This gain comes from averaging across items. It says how well the model recovers an across-item average, not how well it estimates any single item.

Parameter Scenarios

Recovery performance varies with the true parameter configuration. Some scenarios are inherently harder to estimate:

# Test different true parameter configurations
scenarios <- list(
  balanced = c(gg = 0.35, gk = 0.30, kk = 0.35, gamma = 0.25),
  high_learning = c(gg = 0.20, gk = 0.50, kk = 0.30, gamma = 0.25),
  low_learning = c(gg = 0.50, gk = 0.10, kk = 0.40, gamma = 0.25),
  high_guessing = c(gg = 0.35, gk = 0.30, kk = 0.35, gamma = 0.40)
)

scenario_results <- lapply(scenarios, function(params) {
  validate_recovery(params, n = 500, n_items = 2, n_sims = 50, seed = 789)
})

# Compare learning (gk) recovery across scenarios
gk_comparison <- data.frame(
  scenario = names(scenarios),
  true_gk = sapply(scenarios, function(p) p["gk"]),
  bias = sapply(scenario_results, function(r) round(r$bias[r$parameter == "gk"], 4)),
  rmse = sapply(scenario_results, function(r) round(r$rmse[r$parameter == "gk"], 4)),
  sd = sapply(scenario_results, function(r) round(r$estimate_sd[r$parameter == "gk"], 4))
)
print(gk_comparison)
##                       scenario true_gk    bias   rmse     sd
## balanced.gk           balanced     0.3 -0.0019 0.0293 0.0296
## high_learning.gk high_learning     0.5 -0.0028 0.0317 0.0319
## low_learning.gk   low_learning     0.1 -0.0003 0.0221 0.0223
## high_guessing.gk high_guessing     0.3  0.0005 0.0369 0.0373
When Individual-Level Data Adds Information

The item-level model above uses only each item’s transition counts, so it ignores which responses came from the same person. fit_person_lca() uses that information. It assumes each person follows one gg, gk, or kk trajectory across all items, which makes a person’s responses correlated across items:

  • People in the know-know (kk) class answer every item correctly
  • People in the guess-guess (gg) class answer correctly only as often as their guesses succeed
  • A person who is correct on items 1, 2, and 3 is therefore more likely to be in kk than a lucky guesser

Aggregated transition counts cannot see this pattern. Individual responses can, provided the one-trajectory-per-person assumption holds.