Prediction by trial number
Conditions appear largely consistent in predictions across trial
number/across novel Gorps (novel Gorps were presented in fixed order in
prediction trials).
Even in the first trial, there does not appear to be a difference
between conditions.
glmer_infer_sport_trial_1 <-
glmer(infer_sport ~ condition + (1 | participant),
data = data_tidy %>% filter(infer_trial_num == 1),
family = binomial)
# condition difference?
glmer_infer_sport_trial_1 %>%
summary()
On trial 1, there was no evidence children made different predictions
in the skewed vs not skewed conditions about the sport preferences of
novel Gorps (b = 0.03, z = 0.1, p =
0.922).
Predictions: bimodality in responses
To test for bimodality, we treat each participant’s prediction trials
as draws from a Binomial(4, p) process, and ask whether
participants’ responses in the skewed condition are better described by
a single shared p (a unimodal population, with ordinary
binomial sampling noise across the 4 trials) or by a finite mixture of
subpopulations with different p’s (a truly bimodal
population).
Unimodality vs multimodality
library(flexmix) # bimodality analysis, loads modeltools, load after brms models so modeltools doesn't override brms::prior
set.seed(42)
d_mix <- data_infer_summary %>%
ungroup() %>%
mutate(fail = infer_non_na_count - infer_sport_count)
d_mix_skewed <- d_mix %>%
filter(condition == "skewed")
d_mix_not_skewed <- d_mix %>%
filter(condition == "not_skewed")
# fit binomial mixtures with 1, 2, and 3 latent classes
fm_not_skewed <- map(1:3, ~ flexmix(
cbind(infer_sport_count, fail) ~ 1,
data = d_mix_not_skewed,
k = .x,
model = FLXMRglm(family = "binomial")
))
bic_not_skewed <- tibble(
k = 1:3,
AIC = map_dbl(fm_not_skewed, AIC),
BIC = map_dbl(fm_not_skewed, BIC)
)
# class-specific probabilities & sizes for the 2-class model, ordered low to high
not_skewed_comp_p_raw <- plogis(parameters(fm_not_skewed[[2]])) %>% as.numeric()
not_skewed_comp_n_raw <- table(factor(clusters(fm_not_skewed[[2]]), levels = 1:2))
not_skewed_comp_stats <- tibble(p = not_skewed_comp_p_raw, n = as.numeric(not_skewed_comp_n_raw)) %>% arrange(p)
| Binomial mixture model comparison (not skewed condition) |
| k |
AIC |
BIC |
| 1 |
303.9 |
306.5 |
| 2 |
298.6 |
306.5 |
| 3 |
302.6 |
315.8 |
The data in the not skewed condition is best fit by a mixture of two
binomial classes (BIC: 306.5 for a single shared p vs. 306.5
for 2 classes, vs. 315.8 for 3 classes), with estimated classes at
p = 0.54 (n = 82) and p = 1 (n =
23).
# fit binomial mixtures with 1, 2, and 3 latent classes
fm_skewed <- map(1:3, ~ flexmix(
cbind(infer_sport_count, fail) ~ 1,
data = d_mix_skewed,
k = .x,
model = FLXMRglm(family = "binomial")
))
bic_skewed <- tibble(
k = 1:3,
AIC = map_dbl(fm_skewed, AIC),
BIC = map_dbl(fm_skewed, BIC)
)
# class-specific probabilities & sizes for the 2-class model, ordered low to high
skewed_comp_p_raw <- plogis(parameters(fm_skewed[[2]])) %>% as.numeric()
skewed_comp_n_raw <- table(factor(clusters(fm_skewed[[2]]), levels = 1:2))
skewed_comp_stats <- tibble(p = skewed_comp_p_raw, n = as.numeric(skewed_comp_n_raw)) %>% arrange(p)
| Binomial mixture model comparison (skewed condition) |
| k |
AIC |
BIC |
| 1 |
301.4 |
304.0 |
| 2 |
295.8 |
303.7 |
| 3 |
299.8 |
313.0 |
The data in the skewed condition is best fit by a mixture of two
binomial classes (BIC: 304 for a single shared p vs. 303.7 for
2 classes, vs. 313 for 3 classes), with estimated classes at p
= 0.58 (n = 76) and p = 1 (n = 28).
Best unimodal model vs best bimodal mixture of .5 and 1
To test the specific bimodality of exactly .5 (guessing/discounting)
and exactly 1 (generalizing) — we compare a single free-p
binomial model against a 2-class mixture with class probabilities
fixed at .5 and 1 where only the mixing weight is
estimated.
Because the regularity conditions for a standard chi-square
likelihood-ratio test do not hold when comparing mixture models with
different numbers of components, we obtain a null distribution via
parametric bootstrap (simulating from the fitted single-p null
model).
loglik_binom <- function(y, n, p) sum(dbinom(y, n, p, log = TRUE))
# M0: single free p shared by everyone
p_hat <- sum(d_mix_not_skewed$infer_sport_count) / sum(d_mix_not_skewed$infer_non_na_count)
ll_M0 <- loglik_binom(d_mix_not_skewed$infer_sport_count, 4, p_hat)
# M1: 2-class mixture, p fixed at .5 and 1, only mixing weight (pi) estimated
negloglik_mix_fixed <- function(logit_pi, y, n) {
pi_ <- plogis(logit_pi)
lik <- pi_ * dbinom(y, n, 0.5) + (1 - pi_) * dbinom(y, n, 1)
-sum(log(pmax(lik, 1e-300)))
}
opt_M1 <- optimize(negloglik_mix_fixed, interval = c(-10, 10),
y = d_mix_not_skewed$infer_sport_count, n = 4)
pi_hat <- plogis(opt_M1$minimum)
ll_M1 <- -opt_M1$objective
lrt_obs <- 2 * (ll_M1 - ll_M0)
# parametric bootstrap null distribution: simulate under M0, refit both models each time
set.seed(42)
n_boot <- 1000
n_obs <- nrow(d_mix_not_skewed)
lrt_boot <- map_dbl(1:n_boot, function(i) {
y_sim <- rbinom(n_obs, 4, p_hat)
ll0 <- loglik_binom(y_sim, 4, sum(y_sim) / (4 * n_obs))
opt1 <- optimize(negloglik_mix_fixed, interval = c(-10, 10), y = y_sim, n = 4)
2 * (-opt1$objective - ll0)
})
p_boot <- (sum(lrt_boot >= lrt_obs) + 1) / (n_boot + 1)
The fixed .5/1 mixture (mixing weight 0.83 at .5, 0.17 at 1) fits
substantially better than a single shared probability (p =
0.61, LRT = 7.4, parametric-bootstrap p = 0.003). This supports
treating the not skewed condition as a mixture of guessing/discounting
(0.5) and generalizing (1) participants, rather than a single
population.
loglik_binom <- function(y, n, p) sum(dbinom(y, n, p, log = TRUE))
# M0: single free p shared by everyone
p_hat <- sum(d_mix_skewed$infer_sport_count) / sum(d_mix_skewed$infer_non_na_count)
ll_M0 <- loglik_binom(d_mix_skewed$infer_sport_count, 4, p_hat)
# M1: 2-class mixture, p fixed at .5 and 1, only mixing weight (pi) estimated
negloglik_mix_fixed <- function(logit_pi, y, n) {
pi_ <- plogis(logit_pi)
lik <- pi_ * dbinom(y, n, 0.5) + (1 - pi_) * dbinom(y, n, 1)
-sum(log(pmax(lik, 1e-300)))
}
opt_M1 <- optimize(negloglik_mix_fixed, interval = c(-10, 10),
y = d_mix_skewed$infer_sport_count, n = 4)
pi_hat <- plogis(opt_M1$minimum)
ll_M1 <- -opt_M1$objective
lrt_obs <- 2 * (ll_M1 - ll_M0)
# parametric bootstrap null distribution: simulate under M0, refit both models each time
set.seed(42)
n_boot <- 1000
n_obs <- nrow(d_mix_skewed)
lrt_boot <- map_dbl(1:n_boot, function(i) {
y_sim <- rbinom(n_obs, 4, p_hat)
ll0 <- loglik_binom(y_sim, 4, sum(y_sim) / (4 * n_obs))
opt1 <- optimize(negloglik_mix_fixed, interval = c(-10, 10), y = y_sim, n = 4)
2 * (-opt1$objective - ll0)
})
p_boot <- (sum(lrt_boot >= lrt_obs) + 1) / (n_boot + 1)
The fixed .5/1 mixture (mixing weight 0.78 at .5, 0.22 at 1) fits
substantially better than a single shared probability (p =
0.65, LRT = 3.9, parametric-bootstrap p = 0.002). This supports
treating the skewed condition as a mixture of guessing/discounting (0.5)
and generalizing (1) participants, rather than a single population.
Condition differences in mixture
Finally, we test whether class membership (0.5 vs 1) depends on
condition, by fitting a single 2-class mixture across both conditions
with condition entered as a predictor of class membership (a
“concomitant variable” mixture model), and comparing it via
likelihood-ratio test to a version where condition does not predict
class membership. Because both models fix the number of classes at 2 and
differ only in the (ordinary, regular) concomitant logistic model, a
standard chi-square LRT is valid here.
set.seed(42)
fm_null_mix <- flexmix(cbind(infer_sport_count, fail) ~ 1,
data = d_mix, k = 2,
model = FLXMRglm(family = "binomial"),
concomitant = FLXPmultinom(~ 1))
set.seed(42)
fm_cond_mix <- flexmix(cbind(infer_sport_count, fail) ~ 1,
data = d_mix, k = 2,
model = FLXMRglm(family = "binomial"),
concomitant = FLXPmultinom(~ condition))
ll0 <- logLik(fm_null_mix)
ll1 <- logLik(fm_cond_mix)
lrt_cond <- as.numeric(2 * (ll1 - ll0))
df_cond <- attr(ll1, "df") - attr(ll0, "df")
p_cond <- pchisq(lrt_cond, df = df_cond, lower.tail = FALSE)
# identify which cluster index is the "chance-like" (lower p) one
cond_p <- plogis(parameters(fm_cond_mix)) %>% as.numeric()
chance_cluster <- which.min(cond_p)
d_mix$cluster <- clusters(fm_cond_mix)
cluster_by_condition <- d_mix %>%
count(condition, cluster) %>%
group_by(condition) %>%
mutate(prop = n / sum(n))
skewed_chance_prop <- cluster_by_condition %>%
filter(condition == "skewed", cluster == chance_cluster) %>% pull(prop)
notskewed_chance_prop <- cluster_by_condition %>%
filter(condition == "not_skewed", cluster == chance_cluster) %>% pull(prop)
| Mixture class membership by condition |
| cluster |
n |
prop |
| not_skewed |
| 1 |
82 |
78.1% |
| 2 |
23 |
21.9% |
| skewed |
| 1 |
76 |
73.1% |
| 2 |
28 |
26.9% |
Condition does not significantly predict response pattern (LRT(1) =
0.7, p = 0.398): 73.1% of skewed condition participants fall
into the 0.5 class, compared to 78.1% of not skewed condition
participants.
Prediction vs comparison forced-choice
Participants’ responses to the comparison forced-choice question
appear to be unrelated to their earlier responses to the prediction
trials.
The comparison forced-choice asked who likes the sample sport more:
Alex’s Gorp friends, Gorps on Gorp Planet, or if they like it the
same.
Since 8 out of 8 of Alex’s Gorp friends liked the sample sport,
matching the sample proportion would correspond with predictions of 1
(100%). As a result:
predictions of 0%, 25%, 50%, 75% should correspond with
responding “Alex’s Gorp friends” like the sample sport more
predictions of 100% should correspond with responding that Alex’s
Gorp friends and Gorps on Gorp Planet like the sample sport “the
same”
However, that’s not what we see: if anything, responses of “Alex’s
Gorp friends” liking the sample sport more appear to be more
common, not less common, after predictions of 100% than in other
predictions.
Likewise, if participants were responding consistently across the two
measures, we would expect that:
participants who responded that Alex’s Gorp friends like [sample
sport] “the same” as Gorps on Gorp Planet should show higher
levels of generalization on the prediction trials, since Alex’s Gorp
friends all liked soccer
participants who responded “Alex’s Gorp friends” like soccer more
than Gorps on Gorp Planet should show lower levels of
generalization on the prediction trials.
However, that’s not what we see; generalization appears to be largely
similar whether participants say “the same” or “Alex’s Gorp
friends”.
