We plot the data and can see that there is no obvious large difference between the versions with high and low debt.
d.both_completed %>%
ggplot(aes(x=time/60, fill=high_debt_version)) +
geom_boxplot() +
labs(
title = "Distribution of time measurements for the different debt levels",
subtitle = "Notice! Log10 x-scale",
x ="Time (min)"
) +
scale_y_continuous(breaks = NULL) +
scale_x_log10() +
scale_fill_manual(
name = "Debt level",
labels = c("High debt", "Low debt"),
values = c("#7070FF", "lightblue"),
guide = guide_legend(reverse = TRUE)
)
d.both_completed %>%
pull(time) %>%
summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 442.0 937.5 1358.0 1842.2 2423.8 7298.0
sprintf("Variance: %.0f", var(pull(d.both_completed, time)))
## [1] "Variance: 1880646"
As the variance is much greater than the mean we will use a negative binomial family that allows us to model the variance separately.
We include high_debt_verison
as a predictor in our model as this variable represent the very effect we want to measure. We also include a varying intercept for each individual to prevent the model from learning too much from single participants with extreme measurements.
We iterate over the model until we have sane priors. In this case an intercept that could reasonably fit the data with a decent amount of uncertainty to allow flexibility of the model.
time.with <- extendable_model(
base_name = "time",
base_formula = "time ~ 1 + high_debt_version + (1 | session)",
base_priors = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = d.both_completed,
base_control = list(adapt_delta = 0.95)
)
prior_summary(time.with(only_priors= TRUE))
prior_summary(time.with(sample_prior = "only"))
pp_check(time.with(sample_prior = "only"), nsamples = 200) + scale_x_log10()
We choose a beta prior that allows for large effects (+-25 minutes) but is skeptical to any effects larger than +-10 minutes.
sim.size <- 1000
sim.intercept <- rnorm(sim.size, 7.5, 1)
sim.beta <- rnorm(sim.size, 0, 0.5)
sim.beta.diff <- exp(sim.intercept + sim.beta) - exp(sim.intercept)
sim.beta.diff.min <- sim.beta.diff / 60
data.frame(x = sim.beta.diff.min) %>%
ggplot(aes(x)) +
geom_density() +
xlim(-50, 50) +
labs(
title = "Beta parameter prior influence",
x = "Time difference (min)",
y = "Density"
)
We check the posterior distribution and can see that the model seems to have been able to fit the data well. Sampling seems to also have worked well as Rhat values are close to 1 and the sampling plots look nice.
pp_check(time.with(), nsamples = 200) + scale_x_log10()
summary(time.with())
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(data) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.42 0.15 0.09 0.72 1.01 805 758
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.53 0.15 7.24 7.82 1.00 2118
## high_debt_versionfalse -0.17 0.16 -0.48 0.16 1.00 3907
## Tail_ESS
## Intercept 2363
## high_debt_versionfalse 2982
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.69 1.08 1.93 6.13 1.01 1101 2280
##
## Samples were drawn using sample(hmc). 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(time.with(), ask = FALSE)
# default prior for monotonic predictor
edlvl_prior <- prior(dirichlet(2), class = "simo", coef = "moeducation_level1")
We use loo
to check some possible extensions on the model.
loo_result <- loo(
# Benchmark model(s)
time.with(),
# New model(s)
time.with("work_domain"),
time.with("work_experience_programming.s"),
time.with("work_experience_java.s"),
time.with("education_field"),
time.with("mo(education_level)", edlvl_prior),
time.with("workplace_peer_review"),
time.with("workplace_td_tracking"),
time.with("workplace_pair_programming"),
time.with("workplace_coding_standards"),
time.with("scenario"),
time.with("group")
)
loo_result[2]
## $diffs
## elpd_diff se_diff
## time.with("scenario") 0.0 0.0
## time.with("education_field") -0.5 3.4
## time.with() -0.8 2.4
## time.with("mo(education_level)", edlvl_prior) -0.9 3.5
## time.with("work_experience_java.s") -1.3 2.3
## time.with("workplace_peer_review") -1.3 2.8
## time.with("workplace_pair_programming") -1.4 2.5
## time.with("workplace_coding_standards") -1.7 2.4
## time.with("work_experience_programming.s") -1.9 2.4
## time.with("group") -2.7 2.5
## time.with("workplace_td_tracking") -2.7 1.8
## time.with("work_domain") -3.3 3.1
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.0 7.1
## p_loo 14.2 2.8
## looic 734.0 14.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 714
## (0.5, 0.7] (ok) 10 22.7% 879
## (0.7, 1] (bad) 4 9.1% 39
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_domain")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -369.5 6.9
## p_loo 16.4 2.8
## looic 739.0 13.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 27 61.4% 632
## (0.5, 0.7] (ok) 11 25.0% 307
## (0.7, 1] (bad) 4 9.1% 43
## (1, Inf) (very bad) 2 4.5% 17
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_experience_programming.s")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.1 7.1
## p_loo 15.2 2.8
## looic 736.2 14.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 25 56.8% 595
## (0.5, 0.7] (ok) 16 36.4% 70
## (0.7, 1] (bad) 2 4.5% 132
## (1, Inf) (very bad) 1 2.3% 54
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_experience_java.s")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.5 7.0
## p_loo 14.4 2.7
## looic 735.0 13.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 32 72.7% 750
## (0.5, 0.7] (ok) 8 18.2% 79
## (0.7, 1] (bad) 4 9.1% 150
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("education_field")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.7 6.3
## p_loo 13.0 2.1
## looic 733.5 12.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 33 75.0% 786
## (0.5, 0.7] (ok) 7 15.9% 151
## (0.7, 1] (bad) 4 9.1% 130
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("mo(education_level)", edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.1 6.7
## p_loo 13.0 2.4
## looic 734.3 13.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 619
## (0.5, 0.7] (ok) 15 34.1% 157
## (0.7, 1] (bad) 2 4.5% 65
## (1, Inf) (very bad) 1 2.3% 85
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_peer_review")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.5 6.5
## p_loo 13.1 2.3
## looic 735.0 13.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 33 75.0% 432
## (0.5, 0.7] (ok) 9 20.5% 127
## (0.7, 1] (bad) 2 4.5% 134
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_td_tracking")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.9 8.2
## p_loo 16.2 4.0
## looic 737.8 16.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 772
## (0.5, 0.7] (ok) 11 25.0% 189
## (0.7, 1] (bad) 4 9.1% 68
## (1, Inf) (very bad) 1 2.3% 4
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_pair_programming")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.7 7.1
## p_loo 14.9 2.9
## looic 735.3 14.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 849
## (0.5, 0.7] (ok) 9 20.5% 678
## (0.7, 1] (bad) 5 11.4% 81
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_coding_standards")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.9 6.9
## p_loo 14.4 2.7
## looic 735.8 13.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 625
## (0.5, 0.7] (ok) 11 25.0% 462
## (0.7, 1] (bad) 3 6.8% 47
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.2 8.1
## p_loo 17.4 3.8
## looic 732.4 16.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 789
## (0.5, 0.7] (ok) 12 27.3% 211
## (0.7, 1] (bad) 1 2.3% 760
## (1, Inf) (very bad) 3 6.8% 9
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("group")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.9 6.8
## p_loo 15.5 2.7
## looic 737.7 13.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 22 50.0% 631
## (0.5, 0.7] (ok) 18 40.9% 153
## (0.7, 1] (bad) 4 9.1% 74
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
loo_result <- loo(
# Benchmark model(s)
time.with(),
time.with("scenario"),
time.with("education_field"),
time.with("mo(education_level)", edlvl_prior),
# New model(s)
time.with(c("scenario", "education_field")),
time.with(c("scenario", "mo(education_level)"), edlvl_prior),
time.with(c("education_field", "mo(education_level)"), edlvl_prior)
)
loo_result[2]
## $diffs
## elpd_diff
## time.with(c("scenario", "education_field")) 0.0
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) -0.1
## time.with("scenario") -0.5
## time.with("education_field") -1.0
## time.with() -1.3
## time.with("mo(education_level)", edlvl_prior) -1.4
## time.with(c("education_field", "mo(education_level)"), edlvl_prior) -2.0
## se_diff
## time.with(c("scenario", "education_field")) 0.0
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 1.8
## time.with("scenario") 1.6
## time.with("education_field") 2.6
## time.with() 2.3
## time.with("mo(education_level)", edlvl_prior) 2.9
## time.with(c("education_field", "mo(education_level)"), edlvl_prior) 2.8
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.0 7.1
## p_loo 14.2 2.8
## looic 734.0 14.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 714
## (0.5, 0.7] (ok) 10 22.7% 879
## (0.7, 1] (bad) 4 9.1% 39
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.2 8.1
## p_loo 17.4 3.8
## looic 732.4 16.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 789
## (0.5, 0.7] (ok) 12 27.3% 211
## (0.7, 1] (bad) 1 2.3% 760
## (1, Inf) (very bad) 3 6.8% 9
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("education_field")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.7 6.3
## p_loo 13.0 2.1
## looic 733.5 12.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 33 75.0% 786
## (0.5, 0.7] (ok) 7 15.9% 151
## (0.7, 1] (bad) 4 9.1% 130
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("mo(education_level)", edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.1 6.7
## p_loo 13.0 2.4
## looic 734.3 13.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 619
## (0.5, 0.7] (ok) 15 34.1% 157
## (0.7, 1] (bad) 2 4.5% 65
## (1, Inf) (very bad) 1 2.3% 85
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "education_field"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.7 7.2
## p_loo 15.8 3.0
## looic 731.4 14.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 31 70.5% 305
## (0.5, 0.7] (ok) 10 22.7% 349
## (0.7, 1] (bad) 2 4.5% 70
## (1, Inf) (very bad) 1 2.3% 11
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 7.4
## p_loo 15.3 3.1
## looic 731.6 14.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 65.9% 242
## (0.5, 0.7] (ok) 12 27.3% 149
## (0.7, 1] (bad) 2 4.5% 111
## (1, Inf) (very bad) 1 2.3% 15
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("education_field", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.7 6.6
## p_loo 13.8 2.5
## looic 735.4 13.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 575
## (0.5, 0.7] (ok) 11 25.0% 419
## (0.7, 1] (bad) 3 6.8% 38
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
loo_result <- loo(
# Benchmark model(s)
time.with(),
time.with("scenario"),
time.with(c("scenario", "education_field")),
time.with(c("scenario", "mo(education_level)"), edlvl_prior),
# New model(s)
time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior)
)
loo_result[2]
## $diffs
## elpd_diff
## time.with(c("scenario", "education_field")) 0.0
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) -0.1
## time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior) -0.2
## time.with("scenario") -0.5
## time.with() -1.3
## se_diff
## time.with(c("scenario", "education_field")) 0.0
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 1.8
## time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior) 1.6
## time.with("scenario") 1.6
## time.with() 2.3
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.0 7.1
## p_loo 14.2 2.8
## looic 734.0 14.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 714
## (0.5, 0.7] (ok) 10 22.7% 879
## (0.7, 1] (bad) 4 9.1% 39
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.2 8.1
## p_loo 17.4 3.8
## looic 732.4 16.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 789
## (0.5, 0.7] (ok) 12 27.3% 211
## (0.7, 1] (bad) 1 2.3% 760
## (1, Inf) (very bad) 3 6.8% 9
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "education_field"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.7 7.2
## p_loo 15.8 3.0
## looic 731.4 14.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 31 70.5% 305
## (0.5, 0.7] (ok) 10 22.7% 349
## (0.7, 1] (bad) 2 4.5% 70
## (1, Inf) (very bad) 1 2.3% 11
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 7.4
## p_loo 15.3 3.1
## looic 731.6 14.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 65.9% 242
## (0.5, 0.7] (ok) 12 27.3% 149
## (0.7, 1] (bad) 2 4.5% 111
## (1, Inf) (very bad) 1 2.3% 15
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.9 7.0
## p_loo 15.7 2.7
## looic 731.8 13.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 65.9% 363
## (0.5, 0.7] (ok) 12 27.3% 247
## (0.7, 1] (bad) 3 6.8% 25
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
We pick some of our top performing models as candidates and inspect them closer.
The candidate models are named and listed in order of complexity.
We select the simplest model as a baseline.
time0 <- brm(
"time ~ 1 + high_debt_version + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time0",
file_refit = "on_change",
seed = 20210421
)
summary(time0)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.42 0.15 0.09 0.72 1.01 805 758
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.53 0.15 7.24 7.82 1.00 2118
## high_debt_versionfalse -0.17 0.16 -0.48 0.16 1.00 3907
## Tail_ESS
## Intercept 2363
## high_debt_versionfalse 2982
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.69 1.08 1.93 6.13 1.01 1101 2280
##
## Samples were drawn using sample(hmc). 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).
ranef(time0)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.11669471 0.2912992 -0.691603450 0.4837491
## 6033d90a5af2c702367b3a96 0.30711870 0.2865791 -0.199448750 0.9029927
## 6034fc165af2c702367b3a98 0.31371085 0.2888299 -0.199371675 0.9129603
## 603500725af2c702367b3a99 -0.50344864 0.3793934 -1.283707500 0.1683801
## 603f97625af2c702367b3a9d -0.18807193 0.2909735 -0.767492400 0.3754829
## 603fd5d95af2c702367b3a9e 0.23546051 0.2878865 -0.276078525 0.8424215
## 60409b7b5af2c702367b3a9f 0.36291817 0.2989148 -0.156544550 0.9821946
## 604b82b5a7718fbed181b336 -0.32248465 0.3219711 -0.968571825 0.2755393
## 6050c1bf856f36729d2e5218 0.11419944 0.2803647 -0.409474000 0.7149726
## 6050e1e7856f36729d2e5219 0.07260779 0.2756596 -0.452757625 0.6871243
## 6055fdc6856f36729d2e521b -0.07518964 0.2821304 -0.623761075 0.4730638
## 60589862856f36729d2e521f 0.03947516 0.2755802 -0.500501450 0.5974026
## 605afa3a856f36729d2e5222 -0.07249313 0.2768240 -0.626371350 0.4960323
## 605c8bc6856f36729d2e5223 -0.34994406 0.3333765 -1.027536250 0.2724078
## 605f3f2d856f36729d2e5224 -0.11444664 0.2896482 -0.691027850 0.4698477
## 605f46c3856f36729d2e5225 -0.24739263 0.3076050 -0.869087925 0.3130423
## 60605337856f36729d2e5226 0.22880427 0.2834237 -0.289846375 0.8274207
## 60609ae6856f36729d2e5228 0.31922489 0.2973446 -0.208554525 0.9363410
## 6061ce91856f36729d2e522e -0.04307201 0.2846309 -0.605775850 0.5326793
## 6061f106856f36729d2e5231 -0.36589979 0.3323397 -1.027052500 0.2233296
## 6068ea9f856f36729d2e523e 0.65513541 0.3416769 0.008649201 1.3429870
## 6075ab05856f36729d2e5247 -0.29401594 0.3150361 -0.919818325 0.2887394
plot(time0, ask = FALSE)
pp_check(time0, nsamples = 200) + scale_x_log10()
We select the best performing model with one variable.
time1 <- brm(
"time ~ 1 + high_debt_version + scenario + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time1",
file_refit = "on_change",
seed = 20210421
)
summary(time1)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.46 0.14 0.14 0.74 1.00 963 849
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.63 0.16 7.32 7.96 1.00 2872
## high_debt_versionfalse -0.14 0.15 -0.44 0.16 1.00 4792
## scenariotickets -0.28 0.16 -0.58 0.03 1.00 4967
## Tail_ESS
## Intercept 2454
## high_debt_versionfalse 2950
## scenariotickets 2646
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.27 1.29 2.19 7.09 1.00 1187 1947
##
## Samples were drawn using sample(hmc). 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).
ranef(time1)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.115583610 0.2830331 -0.66478810 0.4511145
## 6033d90a5af2c702367b3a96 0.341514218 0.2832043 -0.17783672 0.9206645
## 6034fc165af2c702367b3a98 0.321268916 0.2871339 -0.19579905 0.9213725
## 603500725af2c702367b3a99 -0.587579744 0.3745017 -1.29527900 0.1316325
## 603f97625af2c702367b3a9d -0.181424732 0.2990057 -0.77470498 0.3992026
## 603fd5d95af2c702367b3a9e 0.229892791 0.2775679 -0.27245007 0.7920238
## 60409b7b5af2c702367b3a9f 0.444775895 0.3052817 -0.08956376 1.0837740
## 604b82b5a7718fbed181b336 -0.380108197 0.3214076 -1.00503625 0.2315264
## 6050c1bf856f36729d2e5218 0.182313790 0.2831587 -0.34961053 0.7680736
## 6050e1e7856f36729d2e5219 0.092090669 0.2812478 -0.42394600 0.6795051
## 6055fdc6856f36729d2e521b -0.065346367 0.2865710 -0.62094467 0.5018658
## 60589862856f36729d2e521f 0.007216917 0.2834994 -0.54521260 0.5822809
## 605afa3a856f36729d2e5222 -0.088390301 0.2857242 -0.66521785 0.4863975
## 605c8bc6856f36729d2e5223 -0.425262685 0.3363294 -1.07589375 0.2135765
## 605f3f2d856f36729d2e5224 -0.074144524 0.2956296 -0.66622725 0.5315819
## 605f46c3856f36729d2e5225 -0.280579700 0.3136488 -0.89609735 0.3224625
## 60605337856f36729d2e5226 0.266747824 0.2774766 -0.23797647 0.8315732
## 60609ae6856f36729d2e5228 0.369926537 0.2933434 -0.15128560 0.9724369
## 6061ce91856f36729d2e522e -0.048775492 0.2776012 -0.58712727 0.5234116
## 6061f106856f36729d2e5231 -0.397066789 0.3202059 -1.01565575 0.2025224
## 6068ea9f856f36729d2e523e 0.805526257 0.3434811 0.09021778 1.4881535
## 6075ab05856f36729d2e5247 -0.351681195 0.3161716 -0.97611828 0.2520101
plot(time1, ask = FALSE)
pp_check(time1, nsamples = 200) + scale_x_log10()
We select the best performing model with two variables.
time2 <- brm(
"time ~ 1 + high_debt_version + scenario + mo(education_level) + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape"),
prior(dirichlet(2), class = "simo", coef = "moeducation_level1")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time2",
file_refit = "on_change",
seed = 20210421
)
summary(time2)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + mo(education_level) + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.33 0.16 0.02 0.64 1.01 352 949
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 8.21 0.36 7.55 8.96 1.00 1655
## high_debt_versionfalse -0.14 0.16 -0.45 0.19 1.00 4879
## scenariotickets -0.30 0.16 -0.60 0.02 1.00 5242
## moeducation_level -0.23 0.11 -0.45 -0.00 1.00 1603
## Tail_ESS
## Intercept 2207
## high_debt_versionfalse 2613
## scenariotickets 2457
## moeducation_level 2143
##
## Simplex Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## moeducation_level1[1] 0.34 0.16 0.06 0.66 1.00 3197
## moeducation_level1[2] 0.16 0.11 0.02 0.43 1.00 3526
## moeducation_level1[3] 0.17 0.12 0.02 0.48 1.00 3267
## moeducation_level1[4] 0.32 0.15 0.07 0.62 1.00 3538
## Tail_ESS
## moeducation_level1[1] 2541
## moeducation_level1[2] 2341
## moeducation_level1[3] 2551
## moeducation_level1[4] 2814
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.02 1.22 2.14 6.91 1.01 673 2254
##
## Samples were drawn using sample(hmc). 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).
ranef(time2)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.07799140 0.2459130 -0.5870776 0.4102855
## 6033d90a5af2c702367b3a96 0.15129263 0.2568286 -0.2956168 0.7387291
## 6034fc165af2c702367b3a98 0.24592994 0.2712374 -0.1890150 0.8442728
## 603500725af2c702367b3a99 -0.47241977 0.4083291 -1.3551785 0.1242327
## 603f97625af2c702367b3a9d -0.10958302 0.2592202 -0.6531221 0.3989434
## 603fd5d95af2c702367b3a9e 0.05909038 0.2458195 -0.4085982 0.6090834
## 60409b7b5af2c702367b3a9f 0.34405974 0.2937080 -0.1035594 0.9744387
## 604b82b5a7718fbed181b336 -0.14955346 0.2768897 -0.7727246 0.3532587
## 6050c1bf856f36729d2e5218 0.15407439 0.2545243 -0.2858805 0.7288006
## 6050e1e7856f36729d2e5219 -0.03091432 0.2504365 -0.5365994 0.5047280
## 6055fdc6856f36729d2e521b -0.09465205 0.2544276 -0.6322817 0.3856729
## 60589862856f36729d2e521f 0.14675789 0.2624308 -0.3375379 0.7067854
## 605afa3a856f36729d2e5222 -0.04880437 0.2369269 -0.5407007 0.4302414
## 605c8bc6856f36729d2e5223 -0.17737440 0.2755683 -0.7863273 0.3013425
## 605f3f2d856f36729d2e5224 -0.04445211 0.2446502 -0.5406578 0.4555884
## 605f46c3856f36729d2e5225 -0.06821017 0.2633283 -0.6390652 0.4657307
## 60605337856f36729d2e5226 0.21621305 0.2658024 -0.2005363 0.8220544
## 60609ae6856f36729d2e5228 0.22241851 0.2666259 -0.2195490 0.8274027
## 6061ce91856f36729d2e522e -0.13184735 0.2626390 -0.6978558 0.3636380
## 6061f106856f36729d2e5231 -0.35758087 0.3401545 -1.0930725 0.1588716
## 6068ea9f856f36729d2e523e 0.33244322 0.3398458 -0.1933651 1.0666855
## 6075ab05856f36729d2e5247 -0.12851155 0.2746072 -0.7264079 0.3938674
plot(time2, ask = FALSE)
pp_check(time2, nsamples = 200) + scale_x_log10()
We select the second best performing model with two variables.
time3 <- brm(
"time ~ 1 + high_debt_version + scenario + education_field + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time3",
file_refit = "on_change",
seed = 20210421
)
summary(time3)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + education_field + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.37 0.16 0.04 0.69 1.00 654 1004
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 7.65 0.22 7.24 8.10 1.00
## high_debt_versionfalse -0.15 0.15 -0.44 0.16 1.00
## scenariotickets -0.30 0.16 -0.60 0.01 1.00
## education_fieldInteractionDesign -0.29 0.39 -1.05 0.47 1.00
## education_fieldNone 0.63 0.38 -0.13 1.35 1.00
## education_fieldSoftwareEngineering -0.02 0.23 -0.49 0.41 1.00
## Bulk_ESS Tail_ESS
## Intercept 2391 2360
## high_debt_versionfalse 4297 2550
## scenariotickets 3662 2730
## education_fieldInteractionDesign 3097 2537
## education_fieldNone 2331 2114
## education_fieldSoftwareEngineering 1851 1954
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.14 1.26 2.22 6.97 1.00 983 2607
##
## Samples were drawn using sample(hmc). 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).
ranef(time3)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.10769731 0.2802576 -0.6898158 0.4575978
## 6033d90a5af2c702367b3a96 0.29285080 0.2861934 -0.1731885 0.9262193
## 6034fc165af2c702367b3a98 0.26715952 0.2771937 -0.2065673 0.8464921
## 603500725af2c702367b3a99 -0.45806557 0.3653824 -1.2124855 0.1292843
## 603f97625af2c702367b3a9d -0.15447936 0.2813858 -0.7488767 0.3671758
## 603fd5d95af2c702367b3a9e 0.19417344 0.2674688 -0.2758848 0.7738807
## 60409b7b5af2c702367b3a9f 0.37195187 0.3039285 -0.1191115 1.0007145
## 604b82b5a7718fbed181b336 -0.29808745 0.3135078 -0.9463931 0.2457523
## 6050c1bf856f36729d2e5218 0.15781615 0.2635294 -0.3135425 0.7373020
## 6050e1e7856f36729d2e5219 0.08536319 0.2579874 -0.4051105 0.6485428
## 6055fdc6856f36729d2e521b -0.06246553 0.2699551 -0.6186508 0.4766975
## 60589862856f36729d2e521f 0.01726464 0.2663124 -0.4930719 0.5755771
## 605afa3a856f36729d2e5222 -0.08152158 0.2672625 -0.6297803 0.4293375
## 605c8bc6856f36729d2e5223 -0.32157773 0.3156398 -0.9759616 0.2095555
## 605f3f2d856f36729d2e5224 -0.06939428 0.2671083 -0.6307085 0.4555254
## 605f46c3856f36729d2e5225 -0.21271110 0.2875932 -0.8098952 0.3154569
## 60605337856f36729d2e5226 0.22496683 0.2638384 -0.2326179 0.7869825
## 60609ae6856f36729d2e5228 0.29591347 0.2900156 -0.1822160 0.9207212
## 6061ce91856f36729d2e522e -0.03812971 0.2522389 -0.5482349 0.4754588
## 6061f106856f36729d2e5231 -0.31152541 0.3094848 -0.9629391 0.2278749
## 6068ea9f856f36729d2e523e 0.33673487 0.3661209 -0.2388898 1.1750208
## 6075ab05856f36729d2e5247 -0.15558361 0.3486059 -0.9060679 0.5214046
plot(time3, ask = FALSE)
pp_check(time3, nsamples = 200) + scale_x_log10()
All candidate models look nice, none is significantly better than the others, we will proceed the simplest model: time0
We will try a few different variations of the selected candidate model.
Some participants did only complete one scenario. Those has been excluded from the initial dataset to improve sampling of the models. We do however want to use all data we can and will therefore try to fit the model with the complete dataset.
time0.all <- brm(
"time ~ 1 + high_debt_version + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.completed),
control = list(adapt_delta = 0.95),
file = "fits/time0.all",
file_refit = "on_change",
seed = 20210421
)
summary(time0.all)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(d.completed) (Number of observations: 51)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 29)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.14 0.14 0.70 1.01 743 953
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.53 0.14 7.26 7.81 1.00 3293
## high_debt_versionfalse -0.19 0.15 -0.48 0.12 1.01 5731
## Tail_ESS
## Intercept 2815
## high_debt_versionfalse 2987
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.74 1.07 2.05 6.18 1.01 932 2112
##
## Samples were drawn using sample(hmc). 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).
ranef(time0.all)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033c6fc5af2c702367b3a93 0.10410442 0.3151836 -0.48023627 0.7730021
## 6033d69a5af2c702367b3a95 -0.10868422 0.2866712 -0.66207277 0.4576609
## 6033d90a5af2c702367b3a96 0.32835057 0.3068929 -0.21397340 0.9810650
## 6034fc165af2c702367b3a98 0.32933741 0.2922653 -0.19479717 0.9538842
## 603500725af2c702367b3a99 -0.51782398 0.3760148 -1.26685100 0.1698214
## 603f84f15af2c702367b3a9b -0.37150250 0.4122591 -1.23816400 0.3610192
## 603f97625af2c702367b3a9d -0.18481010 0.2922115 -0.74686020 0.4005106
## 603fd5d95af2c702367b3a9e 0.24770264 0.2845793 -0.28087505 0.8394540
## 60409b7b5af2c702367b3a9f 0.38807890 0.3019859 -0.14703702 1.0172445
## 604b82b5a7718fbed181b336 -0.33546966 0.3268299 -0.97231285 0.2820819
## 604f1239a7718fbed181b33f -0.02264884 0.3374812 -0.69265603 0.6616463
## 6050c1bf856f36729d2e5218 0.12196576 0.2715580 -0.36594708 0.6914848
## 6050e1e7856f36729d2e5219 0.07211696 0.2728541 -0.45671998 0.6409706
## 6055fdc6856f36729d2e521b -0.07050498 0.2822568 -0.61090290 0.4958888
## 60579f2a856f36729d2e521e -0.05580746 0.3367583 -0.71705032 0.6252733
## 60589862856f36729d2e521f 0.04952168 0.2816452 -0.47759367 0.6338071
## 605a30a7856f36729d2e5221 -0.20353520 0.3658427 -0.94769972 0.4936266
## 605afa3a856f36729d2e5222 -0.06548215 0.2860643 -0.62135292 0.5254119
## 605c8bc6856f36729d2e5223 -0.35909178 0.3369657 -1.02624150 0.2814192
## 605f3f2d856f36729d2e5224 -0.10918111 0.2856132 -0.66585697 0.4825179
## 605f46c3856f36729d2e5225 -0.25262460 0.2960580 -0.83618012 0.3307429
## 60605337856f36729d2e5226 0.24323437 0.2845441 -0.26955008 0.8498517
## 60609ae6856f36729d2e5228 0.33646945 0.2956213 -0.17501068 0.9464874
## 6061ce91856f36729d2e522e -0.02992564 0.2887877 -0.58879340 0.5731751
## 6061f106856f36729d2e5231 -0.37765650 0.3213749 -1.01890375 0.2388714
## 60672faa856f36729d2e523c -0.11571495 0.3437219 -0.81626365 0.5768576
## 6068ea9f856f36729d2e523e 0.68947070 0.3350446 0.05149223 1.3585685
## 606db69d856f36729d2e5243 0.56246999 0.3622623 -0.05314622 1.3209190
## 6075ab05856f36729d2e5247 -0.30504090 0.3133298 -0.93250217 0.2894201
plot(time0.all, ask = FALSE)
pp_check(time0.all, nsamples = 200) + scale_x_log10()
As including all data points didn’t harm the model we will create this variant with all data points as well.
This variation includes work_experience_programming.s
predictors as it can give further insight into how experience play a factor in the effect we try to measure. This is especially important as our sampling shewed towards containing less experienced developer than the population at large.
time0.all.exp <- brm(
"time ~ 1 + high_debt_version + work_experience_programming.s + (1 | session)",
prior = c(
prior(normal(0, 0.5), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.completed),
control = list(adapt_delta = 0.95),
file = "fits/time0.all.exp",
file_refit = "on_change",
seed = 20210421
)
summary(time0.all.exp)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + work_experience_programming.s + (1 | session)
## Data: as.data.frame(d.completed) (Number of observations: 51)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 29)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.45 0.13 0.18 0.71 1.00 954 1083
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 7.52 0.14 7.25 7.81 1.00
## high_debt_versionfalse -0.19 0.16 -0.50 0.11 1.00
## work_experience_programming.s 0.04 0.12 -0.19 0.28 1.00
## Bulk_ESS Tail_ESS
## Intercept 2786 2903
## high_debt_versionfalse 3965 3087
## work_experience_programming.s 2179 2083
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.77 1.04 2.08 6.09 1.00 1167 1828
##
## Samples were drawn using sample(hmc). 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).
ranef(time0.all.exp)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033c6fc5af2c702367b3a93 0.124888272 0.3421160 -0.50859793 0.8153015
## 6033d69a5af2c702367b3a95 -0.103583827 0.2947649 -0.67491387 0.4948020
## 6033d90a5af2c702367b3a96 0.353909619 0.2968758 -0.19523363 0.9688395
## 6034fc165af2c702367b3a98 0.359055718 0.2996807 -0.18926650 0.9761628
## 603500725af2c702367b3a99 -0.533716370 0.3649502 -1.24922025 0.1524968
## 603f84f15af2c702367b3a9b -0.386483543 0.4149624 -1.20548575 0.3733449
## 603f97625af2c702367b3a9d -0.180252412 0.3004508 -0.76661150 0.4211721
## 603fd5d95af2c702367b3a9e 0.272772291 0.2911299 -0.24574110 0.8999739
## 60409b7b5af2c702367b3a9f 0.420570695 0.3029043 -0.12626617 1.0546725
## 604b82b5a7718fbed181b336 -0.343519413 0.3262546 -0.98508110 0.2932244
## 604f1239a7718fbed181b33f -0.018207584 0.3413515 -0.68189342 0.6554163
## 6050c1bf856f36729d2e5218 0.121943599 0.2814906 -0.40612295 0.7027546
## 6050e1e7856f36729d2e5219 0.090362926 0.2894689 -0.44775650 0.6866958
## 6055fdc6856f36729d2e521b -0.064477546 0.2944380 -0.63331800 0.5327611
## 60579f2a856f36729d2e521e -0.054081352 0.3583681 -0.76609610 0.6857650
## 60589862856f36729d2e521f 0.001986506 0.3392792 -0.65973057 0.6738390
## 605a30a7856f36729d2e5221 -0.218435513 0.3631870 -0.95567822 0.4729102
## 605afa3a856f36729d2e5222 -0.111097583 0.3238058 -0.73092248 0.5279057
## 605c8bc6856f36729d2e5223 -0.389694078 0.3373437 -1.04876575 0.2457953
## 605f3f2d856f36729d2e5224 -0.180838763 0.3934267 -1.01742100 0.5682156
## 605f46c3856f36729d2e5225 -0.253352831 0.3111771 -0.86839082 0.3532523
## 60605337856f36729d2e5226 0.271966547 0.2946398 -0.26571985 0.8729451
## 60609ae6856f36729d2e5228 0.362958794 0.3002797 -0.17412855 0.9797948
## 6061ce91856f36729d2e522e -0.027025459 0.2930384 -0.59551745 0.5924932
## 6061f106856f36729d2e5231 -0.389104894 0.3229276 -1.02496150 0.2317015
## 60672faa856f36729d2e523c -0.113602645 0.3509402 -0.81445963 0.5772171
## 6068ea9f856f36729d2e523e 0.726095083 0.3271202 0.07378635 1.3815093
## 606db69d856f36729d2e5243 0.563458647 0.3679828 -0.10679083 1.3267103
## 6075ab05856f36729d2e5247 -0.315713976 0.3179565 -0.94020947 0.2775328
loo(
time0.all,
time0.all.exp
)
## Output of model 'time0.all':
##
## Computed from 4000 by 51 log-likelihood matrix
##
## Estimate SE
## elpd_loo -425.4 7.8
## p_loo 17.7 3.1
## looic 850.7 15.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 56.9% 542
## (0.5, 0.7] (ok) 15 29.4% 133
## (0.7, 1] (bad) 7 13.7% 43
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'time0.all.exp':
##
## Computed from 4000 by 51 log-likelihood matrix
##
## Estimate SE
## elpd_loo -425.5 7.6
## p_loo 18.2 3.0
## looic 850.9 15.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 33 64.7% 820
## (0.5, 0.7] (ok) 12 23.5% 739
## (0.7, 1] (bad) 6 11.8% 52
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## time0.all 0.0 0.0
## time0.all.exp -0.1 0.5
plot(time0.all.exp, ask = FALSE)
pp_check(time0.all.exp, nsamples = 200) + scale_x_log10()
This means that our final model, with all data points and experience predictors, is time0.all.exp
To begin interpreting the model we look at how it’s parameters were estimated. As our research is focused on how the outcome of the model is effected we will mainly analyze the \(\beta\) parameters.
mcmc_areas(time0.all.exp, pars = c("b_high_debt_versionfalse", "b_work_experience_programming.s"), prob = 0.95) + scale_y_discrete() +
scale_y_discrete(labels=c("High debt version: false", "Professional programming experience")) +
ggtitle("Beta parameters densities in time model", subtitle = "Shaded region marks 95% of the density. Line marks the median")
We start by extracting posterior samples
scale_programming_experience <- function(x) {
(x - mean(d.completed$work_experience_programming))/ sd(d.completed$work_experience_programming)
}
unscale_programming_experience <- function(x) {
x * sd(d.completed$work_experience_programming) + mean(d.completed$work_experience_programming)
}
post_settings <- expand.grid(
high_debt_version = c("false", "true"),
session = NA,
work_experience_programming.s = sapply(c(0, 3, 10, 25, 40), scale_programming_experience)
)
post <- posterior_predict(time0.all.exp, newdata = post_settings) %>%
melt(value.name = "estimate", varnames = c("sample_number", "settings_id")) %>%
left_join(
rowid_to_column(post_settings, var= "settings_id"),
by = "settings_id"
) %>%
mutate(work_experience_programming = unscale_programming_experience(work_experience_programming.s)) %>%
select(
estimate,
high_debt_version,
work_experience_programming
) %>%
mutate(estimate = estimate/60)
ggplot(post, aes(x=estimate, fill = high_debt_version)) +
geom_density(alpha = 0.5) +
scale_x_log10() +
scale_fill_manual(
name = "Debt version",
labels = c("Low debt", "High debt"),
values = c("lightblue", "darkblue")
) +
facet_grid(rows = vars(work_experience_programming)) +
labs(
title = "Time to complete task / years of programming experience",
subtitle = "Notice! x-axis is log10 scaled.",
x = "Time (min)",
y = "Density"
)
post.diff <- post %>% filter(high_debt_version == "true")
post.diff$estimate = post.diff$estimate - filter(post, high_debt_version == "false")$estimate
ggplot(post.diff, aes(x=estimate, y = 0, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = hdi, vline_linetype = 2) +
scale_fill_manual(name = "HDI", labels = c("100%", "95%", "100%"), values = c("transparent", "lightblue", "transparent"),) +
xlim(-100, 100) +
facet_grid(rows = vars(work_experience_programming)) +
labs(
title = "Time to complete task / years of programming experience",
subtitle = "Difference as: high debt time - low debt time",
x = "Time (min)",
y = "Density"
)
As the effect is neglectable we will not compute any specific probabilities.