Looking at relationship between vocabulary and grammar.
Seamus Donnelly
2022-01-04
1 Preliminaries
library(tidyverse)
library(brms)
library(patchwork)
library(rstatix)1.1 Open datasets and generate summary statistics.
These data sets contain the vocabulary and grammar factor scores from the IRT models.
vocab <- read_csv("files/vocab.csv")
grammar <- read_csv("files/grammar.csv")
complete <- left_join(grammar, vocab, by="data_id")
complete %>%
dplyr::select(-"data_id") %>%
get_summary_stats(type="common") # vocab only open classed items (nouns and predicates)One Grammar LV wasn’t matched to Vocabulary (they were missing a value for 1 word; see script on vocabulary) I doubt it makes a difference.
complete <- complete %>%
filter(!is.na(Total_Vocab))grammar_dens <- complete %>%
dplyr::select("M1_F", "M2_F") %>%
rename(M1 = M1_F, M2=M2_F) %>%
pivot_longer(everything(), names_to="Model", values_to="Theta") %>%
ggplot(aes(x=Theta, fill=Model)) + geom_density(alpha = .2) +
ggtitle("Distribution of Theta") + theme_minimal()
vocab_dens <- complete %>%
rename(Theta=Vocab_F) %>%
ggplot(aes(x=Theta)) + geom_density(alpha = .2) +
ggtitle("Distribution of Theta") + theme_minimal()grammar_se <- complete %>%
dplyr::select("M1_F", "M2_F", "M1_SE", "M2_SE") %>%
pivot_longer(everything(),
names_to = c("Model", ".value"),
names_pattern = "(.+)_(.+)"
) %>%
rename(Theta=`F`) %>%
ggplot(aes(x=Theta, y=SE, colour=Model)) + geom_point() + ggtitle("Standard Errors") +theme_minimal()
vocab_se <- complete %>%
dplyr::select(Theta= "Vocab_F", SE = "Vocab_SE") %>%
ggplot(aes(x=Theta, y=SE)) + geom_point() + ggtitle("Standard Errors") +theme_minimal()
p1 <- grammar_dens + grammar_se + plot_layout(guides = "collect")
ggsave("plots/thetaSE_g.png", last_plot(), width=6, height=4, units="in", dpi=400)
p2 <- vocab_dens + vocab_se
ggsave("plots/thetaSE_v.png", last_plot(), width=6, height=4, units="in", dpi=400)2 Fit Models to Raw data
Fit linear model to the raw data.
m1 <- brm(total ~ Total_Vocab, data=complete, file="bayes_models/m1.rds")
summary(m1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: total ~ Total_Vocab
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -5.15 0.32 -5.76 -4.54 1.00 3945 3195
## Total_Vocab 0.07 0.00 0.06 0.07 1.00 5112 2956
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 7.21 0.11 7.00 7.43 1.00 1961 1955
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).The posterior predictive distribution fits poorly.
Add smoothing.
m2 <- brm(total ~ s(Total_Vocab), data=complete, file="bayes_models/m2.rds")
summary(m2)## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: total ~ s(Total_Vocab)
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sTotal_Vocab_1) 10.51 3.48 5.57 19.19 1.00 1373 1799
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 11.31 0.14 11.04 11.58 1.00 3847 3182
## sTotal_Vocab_1 76.40 11.54 54.74 99.85 1.00 1720 1946
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 6.55 0.10 6.36 6.75 1.00 4374 2878
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).This still fits poorly, though it’s definitely a bit better.
Let’s compare the two models.
loo(m1, m2)## Output of model 'm1':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -7425.7 33.6
## p_loo 3.0 0.1
## looic 14851.4 67.2
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm2':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -7216.4 45.4
## p_loo 8.2 0.5
## looic 14432.8 90.9
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m2 0.0 0.0
## m1 -209.3 22.6The model with the spline fits better.
Change likelihood to poisson
m3 <- brm(total ~ s(Total_Vocab), family="poisson", data=complete, file="bayes_models/m3.rds")
summary(m3)## Warning: There were 12 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup## Family: poisson
## Links: mu = log
## Formula: total ~ s(Total_Vocab)
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sTotal_Vocab_1) 1.78 0.65 0.88 3.45 1.01 665 1256
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.87 0.01 1.84 1.89 1.00 2775 2346
## sTotal_Vocab_1 11.78 1.80 8.66 15.66 1.00 997 1522
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Plot predictions from all three models against the raw data
p3 <- plot(conditional_effects(m1), points=TRUE)
p4 <- plot(conditional_effects(m2), points=TRUE)
p5 <- plot(conditional_effects(m3), points=TRUE)
3 Fit models full dataset:
3.1 2pl model
Plot factors against one another
ggplot(complete, aes(x=Vocab_F, y=M1_F)) + geom_point() + stat_smooth(method="loess")## `geom_smooth()` using formula 'y ~ x'
Factor scores – no error (just for completeness)
m4 <- brm(M1_F ~ Vocab_F, data=complete, file="bayes_models/m4.rds")
summary(m4) ## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F ~ Vocab_F
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.51 0.01 -0.54 -0.48 1.00 3916 3042
## Vocab_F 0.78 0.01 0.75 0.80 1.00 4086 2906
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.56 0.01 0.54 0.58 1.00 3933 3057
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Add error for grammar factor score.
m5 <- brm(M1_F | mi(M1_SE) ~ Vocab_F, data=complete, save_mevars = TRUE, file="bayes_models/m5.rds")
summary(m5)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ Vocab_F
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.39 0.01 -0.41 -0.36 1.00 2687 2463
## Vocab_F 0.71 0.01 0.68 0.73 1.00 2906 2830
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.43 0.01 0.41 0.45 1.00 2311 2106
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Make shrinkage dataset for plots (see below)
d1 <- tibble(
ID = complete %>% pull(data_id),
G1_obs = complete %>% pull(M1_F),
V1_obs = complete %>% pull(Vocab_F),
G1_est = posterior_summary(m5) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Linear")I can, in principle, add measurement error to both latent variables, but I can’t get this to run yet. Maybe try higher adapt delta, which seems to clean things up a bit, but will take some time to run.
#m6 <- brm(M1_F | mi(M1_SE) ~ me(Vocab_F, Vocab_SE), data=complete, save_mevars = TRUE, control=list(adapt_delta=.999, max_treedepth=15), file="bayes_models/m6.rds")
#summary(m6)
#posterior_summary(m6) %>%
# round(2)Fit model with spline for comparison.
m7 <- brm(M1_F | mi(M1_SE) ~ s(Vocab_F), data=complete, save_mevars = TRUE, file="bayes_models/m7.rds")
summary(m7)## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ s(Vocab_F)
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 1.23 0.41 0.65 2.25 1.00 1260 2392
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.07 0.01 0.05 0.09 1.00 2658 2900
## sVocab_F_1 0.58 1.30 -2.13 3.06 1.00 1434 2275
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.01 0.40 0.44 1.00 1865 2988
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m7)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Create other dataset and make plot comparing shrinkage across the two models.
d2 <- tibble(
ID = complete %>% pull(data_id),
G1_obs = complete %>% pull(M1_F),
V1_obs = complete %>% pull(Vocab_F),
G1_est = posterior_summary(m7) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Spine")
d3 <- bind_rows(d1, d2)
p6 <- ggplot(d3, aes(x=Vocab, y = Grammar)) +
geom_line(aes(group = ID),
size = 1/4) +
geom_point(aes(color = Value)) +
facet_wrap(~Model) + theme_minimal()
ggsave("plots/measure1_gauss.png", last_plot(), height=4, width=6, unit="in", dpi=300)loo(m5, m7)## Output of model 'm5':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2086.2 63.3
## p_loo 12.7 0.8
## looic 4172.4 126.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm7':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2011.1 63.5
## p_loo 25.4 1.5
## looic 4022.1 127.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m7 0.0 0.0
## m5 -75.1 15.9Spline model fits substantially better.
3.1.1 Add Age to see if that adds more information around the ceiling and the floor.
m5age <- brm(M1_F | mi(M1_SE) ~ Vocab_F + age, data=complete, save_mevars = TRUE, file="bayes_models/m5age.rds")
summary(m5age)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ Vocab_F + age
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.51 0.07 -1.65 -1.36 1.00 3943 3440
## Vocab_F 0.58 0.01 0.56 0.61 1.00 3805 3576
## age 0.05 0.00 0.04 0.06 1.00 4215 3392
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.40 0.01 0.39 0.42 1.00 2678 2544
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).m7age <- brm(M1_F | mi(M1_SE) ~ s(Vocab_F) + age, data=complete, save_mevars = TRUE, file="bayes_models/m7age.rds")
summary(m7age)## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ s(Vocab_F) + age
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 1.07 0.38 0.53 2.05 1.00 1341 1973
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.05 0.08 -1.19 -0.89 1.00 4577 3088
## age 0.05 0.00 0.04 0.05 1.00 4693 3040
## sVocab_F_1 0.55 1.24 -2.03 2.82 1.00 1707 2172
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.39 0.01 0.38 0.41 1.00 2339 2385
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).loo(m5age, m7age)## Output of model 'm5age':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2009.8 67.4
## p_loo 17.6 1.0
## looic 4019.5 134.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm7age':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1948.2 67.1
## p_loo 30.0 1.8
## looic 3896.5 134.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 2186 100.0% 603
## (0.5, 0.7] (ok) 1 0.0% 1422
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m7age 0.0 0.0
## m5age -61.5 13.5Roughly the same results.
3.2 EH Model
ggplot(complete, aes(x=Vocab_F, y=M2_F)) + geom_point() + stat_smooth(method="loess")## `geom_smooth()` using formula 'y ~ x'
m8 <- brm(M2_F ~ Vocab_F, data=complete, file="bayes_models/m8.rds")
summary(m8) ## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F ~ Vocab_F
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.38 0.03 -1.44 -1.31 1.00 4763 3298
## Vocab_F 1.71 0.03 1.66 1.76 1.00 4387 3222
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.28 0.02 1.24 1.32 1.00 4371 3152
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).m9 <- brm(M2_F | mi(M2_SE) ~ Vocab_F, save_mevars = TRUE, data=complete, file="bayes_models/m9.rds")
summary(m9)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ Vocab_F
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.13 0.03 -1.20 -1.07 1.00 2444 2795
## Vocab_F 1.58 0.03 1.53 1.64 1.00 2969 3027
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.03 0.02 0.99 1.07 1.00 2164 2669
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).m10 <- brm(M2_F | mi(M2_SE) ~ s(Vocab_F), save_mevars = TRUE, data=complete, file="bayes_models/m10.rds")
summary(m10)## Warning: There were 6 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ s(Vocab_F)
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 2.53 0.81 1.38 4.45 1.00 1639 2291
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.11 0.03 -0.16 -0.06 1.00 5122 3118
## sVocab_F_1 1.46 2.75 -4.24 6.55 1.00 2325 2688
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.00 0.02 0.96 1.04 1.00 2932 3056
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).d1 <- tibble(
ID = complete %>% pull(data_id),
G1_obs = complete %>% pull(M1_F),
V1_obs = complete %>% pull(Vocab_F),
G1_est = posterior_summary(m9) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Linear")d2 <- tibble(
ID = complete %>% pull(data_id),
G1_obs = complete %>% pull(M1_F),
V1_obs = complete %>% pull(Vocab_F),
G1_est = posterior_summary(m10) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Spline")
d3 <- bind_rows(d1, d2)
p7 <- ggplot(d3, aes(x=Vocab, y = Grammar)) +
geom_line(aes(group = ID),
size = 1/4) +
geom_point(aes(color = Value)) +
facet_wrap(~Model) + theme_minimal()
ggsave("plots/measure1_EH.png", last_plot(), height=4, width=6, unit="in", dpi=300)loo(m9, m10)## Output of model 'm9':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -3813.1 59.9
## p_loo 11.0 0.7
## looic 7626.2 119.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm10':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -3751.9 61.7
## p_loo 22.1 1.4
## looic 7503.8 123.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m10 0.0 0.0
## m9 -61.2 14.63.2.1 Add age
m9age <- brm(M2_F | mi(M2_SE) ~ Vocab_F + age, save_mevars = TRUE, data=complete, file="bayes_models/m9age.rds")
summary(m9age)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ Vocab_F + age
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.71 0.17 -4.04 -3.38 1.00 3397 2977
## Vocab_F 1.29 0.03 1.24 1.36 1.00 2981 3033
## age 0.12 0.01 0.10 0.13 1.00 3395 3018
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.96 0.02 0.92 1.00 1.00 2216 2686
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).m10age <- brm(M2_F | mi(M2_SE) ~ s(Vocab_F) + age, save_mevars = TRUE, data=complete, file="bayes_models/m10age.rds")
summary(m10age)## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ s(Vocab_F) + age
## Data: complete (Number of observations: 2187)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 2.21 0.75 1.14 4.10 1.00 1465 1767
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -2.66 0.18 -3.00 -2.32 1.00 5185 3569
## age 0.11 0.01 0.09 0.12 1.00 5509 3634
## sVocab_F_1 1.36 2.58 -4.03 5.94 1.00 2052 2516
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.94 0.02 0.90 0.98 1.00 2558 2958
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m10age)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo(m9age, m10age)## Output of model 'm9age':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -3736.1 64.2
## p_loo 15.4 0.9
## looic 7472.2 128.4
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm10age':
##
## Computed from 4000 by 2187 log-likelihood matrix
##
## Estimate SE
## elpd_loo -3688.7 65.6
## p_loo 26.2 1.6
## looic 7377.5 131.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m10age 0.0 0.0
## m9age -47.4 12.14 Drop floor and ceiling cases.
It’s not clear if the curvilinear relationship above is due to actual curvilinearity or if the ceiling and floor are still pulling the trend around. I’ll try removing observations right at the ceiling and floor and seeing that the results change.
complete %>%
summarise(
min = min(M1_F),
max = max(M1_F)
)4.1 2PL model
m11 <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M1_F | mi(M1_SE) ~ Vocab_F, ., save_mevars = TRUE, file = "bayes_models/m11.rds")
summary(m11)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ Vocab_F
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.18 0.02 -0.21 -0.15 1.00 2796 2782
## Vocab_F 0.55 0.01 0.52 0.57 1.00 4185 3060
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.37 0.01 0.36 0.39 1.00 3078 3092
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).d1 <- tibble(
ID = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(data_id),
G1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(M1_F),
V1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(Vocab_F),
G1_est = posterior_summary(m11) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Linear")m12 <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M1_F | mi(M1_SE) ~ s(Vocab_F), ., save_mevars = TRUE, file = "bayes_models/m12.rds")
summary(m12)## Warning: There were 75 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ s(Vocab_F)
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 1.19 0.45 0.62 2.55 1.01 439 170
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.30 0.01 0.28 0.32 1.00 4092 3072
## sVocab_F_1 0.30 1.08 -1.80 2.48 1.00 2940 2794
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.36 0.01 0.34 0.37 1.00 912 175
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).d2 <- tibble(
ID = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(data_id),
G1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(M1_F),
V1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(Vocab_F),
G1_est = posterior_summary(m12) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "Observed", "Posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Spline")
d3 <- bind_rows(d1, d2)
p8 <- ggplot(d3, aes(x=Vocab, y = Grammar)) +
geom_line(aes(group = ID),
size = 1/4) +
geom_point(aes(color = Value)) +
facet_wrap(~Model) + theme_minimal()
ggsave("plots/measure2_gauss.png", last_plot(), height=4, width=6, unit="in", dpi=300)loo(m11, m12)## Output of model 'm11':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1023.0 38.6
## p_loo 6.4 0.4
## looic 2046.0 77.1
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm12':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -956.8 40.0
## p_loo 16.1 1.3
## looic 1913.6 79.9
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 1644 99.9% 478
## (0.5, 0.7] (ok) 1 0.1% 229
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m12 0.0 0.0
## m11 -66.2 14.94.1.1 Add age
m11age <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M1_F | mi(M1_SE) ~ Vocab_F + age, ., save_mevars = TRUE, file = "bayes_models/m11age.rds")
summary(m11age)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ Vocab_F + age
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.07 0.07 -1.21 -0.93 1.00 4541 3446
## Vocab_F 0.46 0.01 0.44 0.49 1.00 4013 3036
## age 0.04 0.00 0.03 0.05 1.00 4566 3480
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.35 0.01 0.34 0.37 1.00 2801 3137
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m11age)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
m12age <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M1_F | mi(M1_SE) ~ s(Vocab_F) + age, ., save_mevars = TRUE, file = "bayes_models/m12age.rds")
summary(m12age)## Warning: There were 18 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M1_F | mi(M1_SE) ~ s(Vocab_F) + age
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 1.08 0.39 0.56 2.12 1.00 1434 1212
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.57 0.08 -0.72 -0.42 1.00 3657 1613
## age 0.04 0.00 0.03 0.04 1.00 3738 1526
## sVocab_F_1 0.52 1.00 -1.47 2.52 1.00 2756 2941
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.34 0.01 0.33 0.35 1.00 3743 3513
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m12age)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo(m11age, m12age)## Output of model 'm11age':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -960.9 40.5
## p_loo 8.8 0.5
## looic 1921.8 81.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm12age':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -905.5 40.9
## p_loo 19.0 1.8
## looic 1811.1 81.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 1644 99.9% 1320
## (0.5, 0.7] (ok) 1 0.1% 196
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m12age 0.0 0.0
## m11age -55.4 13.54.2 EH Model
m13 <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M2_F | mi(M2_SE) ~ Vocab_F, ., save_mevars = TRUE, file = "bayes_models/m13.rds")
summary(m13)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ Vocab_F
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.35 0.03 -0.41 -0.29 1.00 3602 2852
## Vocab_F 1.01 0.02 0.96 1.05 1.00 3953 2747
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.69 0.01 0.66 0.72 1.00 2690 2892
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).m14 <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M2_F | mi(M2_SE) ~ s(Vocab_F), ., save_mevars = TRUE, file = "bayes_models/m14.rds")
summary(m14)## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ s(Vocab_F)
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 2.02 0.62 1.12 3.58 1.00 1480 2155
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.54 0.02 0.50 0.57 1.00 4749 3304
## sVocab_F_1 0.53 2.00 -3.34 4.60 1.00 2507 2868
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.66 0.01 0.63 0.69 1.00 3759 3210
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).d1 <- tibble(
ID = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(data_id),
G1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(M1_F),
V1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(Vocab_F),
G1_est = posterior_summary(m13) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "observed", "posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Linear")
d2 <- tibble(
ID = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(data_id),
G1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(M1_F),
V1_obs = filter(complete, M1_F > -1 & M1_F < 2) %>% pull(Vocab_F),
G1_est = posterior_summary(m14) %>%
data.frame() %>%
rownames_to_column("term") %>%
filter(str_detect(term, "Yl")) %>%
pull(Estimate)) %>%
pivot_longer(-c(ID, V1_obs), values_to="g") %>%
mutate(
Value = ifelse(name == "G1_obs", "observed", "posterior")
) %>%
rename(Grammar = g, Vocab = V1_obs) %>%
mutate(Model = "Spline")
d3 <- bind_rows(d1, d2)
p9 <- ggplot(d3, aes(x=Vocab, y = Grammar)) +
geom_line(aes(group = ID),
size = 1/4) +
geom_point(aes(color = Value)) +
facet_wrap(~Model) + theme_minimal()
ggsave("plots/measure2_EH.png", last_plot(), height=4, width=6, unit="in", dpi=300)loo(m13, m14)## Output of model 'm13':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2167.7 44.6
## p_loo 8.2 0.5
## looic 4335.3 89.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm14':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2108.7 46.0
## p_loo 19.7 1.6
## looic 4217.4 92.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m14 0.0 0.0
## m13 -59.0 15.64.2.1 Add age
m13age <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M2_F | mi(M2_SE) ~ Vocab_F + age, ., save_mevars = TRUE, file = "bayes_models/m13age.rds")
summary(m13age)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ Vocab_F + age
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -2.02 0.13 -2.27 -1.75 1.00 3794 3115
## Vocab_F 0.85 0.03 0.80 0.90 1.00 3001 3044
## age 0.07 0.01 0.06 0.08 1.00 3695 2982
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.01 0.62 0.68 1.00 2506 2974
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m13age)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
m14age <- complete %>%
filter(M1_F > -1 & M1_F < 2) %>%
brm(M2_F | mi(M2_SE) ~ s(Vocab_F) + age, ., save_mevars = TRUE, file = "bayes_models/m14age.rds")
summary(m14age)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: M2_F | mi(M2_SE) ~ s(Vocab_F) + age
## Data: . (Number of observations: 1645)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sVocab_F_1) 1.89 0.61 1.03 3.40 1.00 1850 2465
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.09 0.14 -1.36 -0.81 1.00 6038 3293
## age 0.07 0.01 0.05 0.08 1.00 6260 3290
## sVocab_F_1 0.95 1.90 -2.73 4.83 1.00 2493 2965
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.62 0.01 0.60 0.65 1.00 3252 3340
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m14age)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
loo(m13age, m14age)## Output of model 'm13age':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2111.2 47.1
## p_loo 10.9 0.6
## looic 4222.3 94.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm14age':
##
## Computed from 4000 by 1645 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2065.3 47.7
## p_loo 23.6 2.1
## looic 4130.5 95.4
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 1644 99.9% 1293
## (0.5, 0.7] (ok) 1 0.1% 121
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m14age 0.0 0.0
## m13age -45.9 14.3#m13 <- complete %>%
# filter(M1_F > -1 & M1_F < 2) %>%
# brm(M1_F | mi(M1_SE) ~ me(Vocab_F, Vocab_SE), ., file = "bayes_models/m13.rds",
# control=list(adapt_delta = .95, max_treedepth = 15))