Replication package for the paper Measuring affective states from technical debt
R. Torkar, J. Olsson, T. Besker, and A. Martini
First created on 2020-06-20. Updated on 2020-09-15.
1 Data preparation
d <- read.csv("../data/data.csv")
dat <- list(
V = d$V,
A = d$A,
D = d$D,
EDU = factor(d$EDU, levels = c("1","2","3","4","5","6","7"), ordered=TRUE), # ordinal 1-7 int (since it's monotonic)
EXAMPLE_idx = d$EXAMPLE, # nominal 1-10 factor
ID_idx = as.integer(d$ID), # subject ID
MAJOR_idx = d$MAJOR, # nominal 1-6 factor
ROLE_idx = d$ROLE, # nominal 1-5 factor
LANG_idx = d$LANG, # nominal 1-10 factor
EXP = scale(d$EXP), # 1-35 years of experience
ENTITIES = scale(d$ENTITIES) # 2-13 logical abstractions
)
str(dat)## List of 11
## $ V : int [1:200] 7 9 5 7 9 9 3 2 9 5 ...
## $ A : int [1:200] 6 4 5 7 5 7 8 3 1 5 ...
## $ D : int [1:200] 7 4 5 3 1 7 7 9 9 7 ...
## $ EDU : Ord.factor w/ 7 levels "1"<"2"<"3"<"4"<..: 5 5 5 5 5 5 5 5 5 5 ...
## $ EXAMPLE_idx: chr [1:200] "DL" "EH" "AH" "CL" ...
## $ ID_idx : int [1:200] 1 1 1 1 1 2 2 2 2 2 ...
## $ MAJOR_idx : chr [1:200] "CS" "CS" "CS" "CS" ...
## $ ROLE_idx : chr [1:200] "D" "D" "D" "D" ...
## $ LANG_idx : chr [1:200] "J" "J" "J" "J" ...
## $ EXP : num [1:200, 1] 0.895 0.895 0.895 0.895 0.895 ...
## ..- attr(*, "scaled:center")= num 10.1
## ..- attr(*, "scaled:scale")= num 7.74
## $ ENTITIES : num [1:200, 1] -0.891 -0.597 2.045 -1.185 1.164 ...
## ..- attr(*, "scaled:center")= num 6.04
## ..- attr(*, "scaled:scale")= num 3.41V, A, D are codes for emotions, feelings, and moods, respectively. For each of the outcomes, \(\{V,A,D\}\), we have recorded a number of predictors (covariates):
- EDU is ordinal (7 levels), and represents participants’ educational level (from none to PhD).
- MAJOR is nominal (6), and represents participants’ education background, e.g. Business Administration, Computer Science, Electrical Engineering, etc.
- ROLE nominal (5), represents participants’ work title: Architect, Developer, Manager, etc.
- EXAMPLE nominal (10), represents which scenario A-E and, additionally, level (i.e., High or Low technical debt).
- LANG nominal (10), programming language most experience in, e.g., Ada, C, C#, etc.
- ID nominal (40). ID of participant.
- EXP ratio (1,35) professional experience in programming (years, rounded down). As is evident above we’ve scaled the variable.
- ENTITIES ratio (2,13). Number of logical abstractions in the example. As is evident above we’ve scaled the variable.
2 Descriptive statistics
Let’s plot our outcomes,
simplehist(dat$V, xlim=c(1,9), ylim=c(0,50),
xlab = "response for V", bty="n")
simplehist(dat$A, xlim=c(1,9), ylim=c(0,50),
xlab = "response for A", bty="n")
simplehist(dat$D, xlim=c(1,9), ylim=c(0,50),
xlab = "response for D", bty="n")
Let’s also look at the cumulative proportions,
# discrete proportion of each response value in V, A, and D
pr_k_v <- table(dat$V) / nrow(d)
pr_k_a <- table(dat$A) / nrow(d)
pr_k_d <- table(dat$D) / nrow(d)
# cumsum converts to cumulative proportions
cum_pr_k_v <- cumsum(pr_k_v)
cum_pr_k_a <- cumsum(pr_k_a)
cum_pr_k_d <- cumsum(pr_k_d)
plot(1:9, cum_pr_k_v, type="b", xlab="response for V",
ylab="", ylim=c(0,1), bty="n", xaxt="n")
axis(1, seq(1,9,1))
plot(1:9, cum_pr_k_a, type="b", xlab="response for A",
ylab="", ylim=c(0,1), bty="n", xaxt="n")
axis(1, seq(1,9,1))
plot(1:9, cum_pr_k_d, type="b", xlab="response for D",
ylab="", ylim=c(0,1), bty="n", xaxt="n")
axis(1, seq(1,9,1))
and log cumulative odds,
# log-cumulative-odds
lco_v <- logit_scaled(cum_pr_k_v)
lco_a <- logit_scaled(cum_pr_k_a)
lco_d <- logit_scaled(cum_pr_k_d)
plot(1:9, lco_v, type="b", ylim = c(-4,4), xlim = c(1,9), xlab="response for V",
ylab="log-cumulative-odds for Emotions (V)", xaxt="n", bty = "n")
axis(1, seq(1,9,1))
plot(1:9, lco_a, type="b", ylim = c(-4,4), xlim = c(1,9), xlab="response for A",
ylab="log-cumulative-odds for Feelings (A)", xaxt="n", bty = "n")
axis(1, seq(1,9,1))
plot(1:9, lco_d, type="b", ylim = c(-4,4), xlim = c(1,9), xlab="response for D",
ylab="log-cumulative-odds for Moods (D)", xaxt="n", bty = "n")
axis(1, seq(1,9,1))
So now we have a feeling for the response variables.
3 Prior sensitivity analysis
Let’s look at the priors now. Later below we’ll do prior predictive checks.
So for each linear model in our multivariate model we set the following priors for now:
- mo(EDU) - Dirichlet prior (we set {2,2,2,2,2} which is very weak)
- EXAMPLE_idx - N(0,1) times 10 factors. Means (1*10)^2 in variance
- MAJOR_idx - N(0,1) times 6 factors.
- ROLE_idx - N(0,1) times 5 factors.
- LANG_idx - N(0,1) times 10 factors.
- ENTITIES - N(0,1)
- EXP - N(0,1)
- ID_idx - varying intercept
We have basically this in our prior sensitivity analysis,
N <- 100
SD <- 0.5 # Changed to 0.1 and 1 also and 0.5 is definitely enough.
intercept <- rnorm(N, 0, 1)
Edu <- rdirichlet(n = N, alpha = c(2,2,2,2,2,2))
Example1 <- rnorm(N, 0, SD)
Example2 <- rnorm(N, 0, SD)
Example3 <- rnorm(N, 0, SD)
Example4 <- rnorm(N, 0, SD)
Example5 <- rnorm(N, 0, SD)
Example6 <- rnorm(N, 0, SD)
Example7 <- rnorm(N, 0, SD)
Example8 <- rnorm(N, 0, SD)
Example9 <- rnorm(N, 0, SD)
Example10 <- rnorm(N, 0, SD)
Major1 <- rnorm(N, 0, SD)
Major2 <- rnorm(N, 0, SD)
Major3 <- rnorm(N, 0, SD)
Major4 <- rnorm(N, 0, SD)
Major5 <- rnorm(N, 0, SD)
Major6 <- rnorm(N, 0, SD)
Role1 <- rnorm(N, 0, SD)
Role2 <- rnorm(N, 0, SD)
Role3 <- rnorm(N, 0, SD)
Role4 <- rnorm(N, 0, SD)
Role5 <- rnorm(N, 0, SD)
Lang1 <- rnorm(N, 0, SD)
Lang2 <- rnorm(N, 0, SD)
Lang3 <- rnorm(N, 0, SD)
Lang4 <- rnorm(N, 0, SD)
Lang5 <- rnorm(N, 0, SD)
Lang6 <- rnorm(N, 0, SD)
Lang7 <- rnorm(N, 0, SD)
Lang8 <- rnorm(N, 0, SD)
Lang9 <- rnorm(N, 0, SD)
Lang10 <- rnorm(N, 0, SD)
Entities1 <- rnorm(N, 0, SD)
Entities2 <- rnorm(N, 0, SD)
Entities3 <- rnorm(N, 0, SD)
Entities4 <- rnorm(N, 0, SD)
Entities5 <- rnorm(N, 0, SD)
Entities6 <- rnorm(N, 0, SD)
Entities7 <- rnorm(N, 0, SD)
Exp <- rnorm(N, 0, SD)
plot(NULL, xlim=c(-2,2), ylim=c(0,1))
for(i in 1:N)
curve(inv_logit_scaled(intercept[i]*x + Edu[i]*x + Example1[i]*x + Example2[i]*x +
Example3[i]*x + Example4[i]*x + Example5[i]*x +
Example6[i]*x +Example7[i]*x + Example8[i]*x +
Example9[i]*x + Example10[i]*x +
Major1[i]*x + Major2[i]*x + Major3[i]*x +
Major4[i]*x + Major5[i]*x + Major6[i]*x +
Role1[i]*x + Role2[i]*x + Role3[i]*x + Role4[i]*x +
Role5[i]*x +
Lang1[i]*x + Lang2[i]*x + Lang3[i]*x + Lang4[i]*x +
Lang5[i]*x + Lang6[i]*x + Lang7[i]*x + Lang8[i]*x +
Lang9[i]*x + Lang10[i]*x +
Entities1[i]*x + Entities2[i]*x + Entities3[i]*x +
Entities4[i]*x + Entities5[i]*x + Entities6[i]*x +
Entities7[i]*x + Exp[i]*x),
add=TRUE)
In short, the above shows how our combined priors look like in the probability space (\(y\)-axis). Due to the gap in the middle one could say we are skeptical of any extreme curves around 0. In the end, prior predictive checks will allow us to investigate this on the outcome scale.
Above we used a Dirichlet prior on ordered categorical data. The Dirichlet is a generalization of the Beta distribution, modeling probability \((0,1)\), but instead of having one probability being estimated, we can have many (Blei, Ng, and Jordan 2003). How many depends on the number of categories \(K\), which means we need to set the prior with \(K-1\) parameters (the last category we can calculate from the previous category).
If we plot different Dirichlet priors we might better understand what it does,
x <- seq(0, 1, 0.001)
p6 <- data.frame(x, y = dbeta(x, 6, 30)) %>%
ggplot(aes(x, y)) + geom_smooth(stat = "identity", size=0.25, col = "black") +
xlab("Dirichlet(6,6,6,6,6)") + ylab("") + coord_cartesian(ylim=c(0, 7)) +
theme_tufte()
p4 <- data.frame(x, y = dbeta(x, 4, 20)) %>%
ggplot(aes(x, y)) + geom_smooth(stat = "identity", size=0.25, col = "black") +
xlab("Dirichlet(4,4,4,4,4)") + ylab("") + coord_cartesian(ylim=c(0, 7)) +
theme_tufte() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p2 <- data.frame(x, y = dbeta(x, 2, 10)) %>%
ggplot(aes(x, y)) + geom_smooth(stat = "identity", size=0.25, col = "black") +
xlab("Dirichlet(2,2,2,2,2)") + ylab("") + coord_cartesian(ylim=c(0, 7)) +
theme_tufte() +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p2 / p4 / p6
So all of the above are for estimating \(5+1\) categories using \(5\) parameters. It’s evident that the higher numbers we set for our parameters (\(2 \rightarrow 4 \rightarrow 6\)) the more mass is on lower probabilities. Having 2’s implies an almost flat prior, i.e., we really don’t have any good prior knowledge on how this simplex would look like so this is it.
4 Model comparisons and prior predictive checks
Since we use the cumulative likelihood we’ll have a hard time using any principled ways of conducting variable selection. Luckily, in our case, we are really not interested in arriving to as small a set of predictors as possible. Hence, we use all predictors and instead do model comparisons while adding them one by one.
First, we execute one model at a time (using default priors only), and then we use LOO to compare their relative out of sampled predictions capabilities (Vehtari, Gelman, and Gabry 2017).
m1 <- brm(mvbind(D,A,V) ~ mo(EDU) + (1|ID_idx),
data=dat, family=cumulative)
m2 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + (1|ID_idx),
data=dat, family=cumulative)
m3 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + (1|ID_idx),
data=dat, family=cumulative)
m4 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
(1|ID_idx),
data=dat, family=cumulative)
m5 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + (1|ID_idx),
data=dat, family=cumulative)
m6 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + ENTITIES + (1|ID_idx),
data=dat, family=cumulative)
m7 <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + ENTITIES + EXP + (1|ID_idx),
data=dat, family=cumulative)loo_compare(loo(m1),loo(m2),loo(m3),loo(m4),loo(m5),loo(m6),loo(m7))## elpd_diff se_diff
## m2 0.0 0.0
## m5 -6.6 8.7
## m3 -6.8 3.0
## m7 -8.0 8.8
## m6 -9.1 8.9
## m4 -9.7 4.8
## m1 -11.3 10.0So, in short, very small differences between the models and if we would consider out of sample prediction capabilities only, then one could opt for the smaller model, i.e., \(\mathcal{M}_2\). However, we are really interested in looking at all independent variables. Hence, \(\mathcal{M}_7\) is our preferred model, \(\mathcal{M}\), which we’ll now run and sample only from the priors (liek we did above in the prior sensitivity analysis).
First, set priors on our parameters we want to estimate and then sample only from the priors. We set \(\mathcal{N}(0, 0.5)\) on \(\beta\), and \(\mathcal{N}(0, 5)\) on \(\alpha\) (since we have a cumulative model we have several \(\alpha\), one for each border between responses, i.e., \(9\) responses means \(8\) \(\alpha\)). We have \(\mathrm{Weibull}(2,1)\) as prior for modeling \(\sigma\) (allowing only positive real numbers) and, finally, our Dirichlet prior on the predictor EDU, which is modeled as monotonic. We recommend you to change the priors and see how this affects the prior probability distribution.
p <- get_prior(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + ENTITIES + EXP + (1|ID_idx),
family = cumulative, data = dat)
p$prior[c(3,47,91)] <- "normal(0,0.5)" # 30 betas for each outcome
p$prior[c(34,78,122)] <- "normal(0,5)" # 8 alphas for each outcome
p$prior[c(43,87,131)] <- "cauchy(0,2)" # 1 sigma for each outcome
p$prior[c(46,90,134)] <- "dirichlet(2,2,2,2,2)" # 1 monotonic for each outcome
M_prio <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + ENTITIES + EXP + (1|ID_idx), prior = p, data=dat,
family=cumulative, sample_prior = "only", seed = SEED)Now plot our priors on the outcome scale,
pp_check(M_prio, type="bars", resp = "V", nsamples = 1000)
The \(y\)-axis are the counts while the \(x\)-axis is the Likert scale outcome for the response \(V\). As we see, the medians for our prior samples (\(y_{\mathrm{rep}}\)) are evenly distributed and the credible intervals are wide (the bars are our empirical data \(y\)). This looks to be a good setup for our model. Now, let’s combine our prior with our data to arrive to the posterior probability distribution,
M <- brm(mvbind(D,A,V) ~ mo(EDU) + EXAMPLE_idx + MAJOR_idx + ROLE_idx +
LANG_idx + ENTITIES + EXP + (1|ID_idx),
prior = p, data=dat, family=cumulative, seed = SEED)5 Diagnostics and posterior predictive checks
A minimum of three diagnostics should preferably always be checked before one puts any trust in a model (and you should never trust a model completely\(\ldots\)) First we can use plot(M) to visualize how well the chains have mixed (they should look like fat hairy caterpillar plots). We won’t plot them here since they will take up a lot of space.
Next we check \(\widehat{R}\) to check the between- and within-variance of the chains. If all is well they should converge towards \(1.0\) when \(N \rightarrow \infty\). As a rule of thumb, if \(\widehat{R} > 1.01\) it’s an indication that the chains have not converged to a stationary posterior probability distribution.
# the discrimination params are NAs, as they should, so we exclude them
max(rhat(M), na.rm = TRUE) ## [1] 1.00442The effective sample size should be larger than \(10\)% of the total sample size,
min(neff_ratio(M), na.rm = TRUE)## [1] 0.1855241So, diagnostics look good. Let’s also do some posterior predictive checks.
pp_check(M, type = "bars", resp = "V", nsamples = NULL)
pp_check(M, type = "bars", resp = "A", nsamples = NULL)
pp_check(M, type = "bars", resp = "D", nsamples = NULL)
According to the function description for ?ppc_bars:
Rootograms allow for diagnosing problems in count data models such as overdispersion or excess zeros. They consist of a histogram of \(y\) with the expected counts based on yrep overlaid as a line along with uncertainty intervals. The \(y\)-axis represents the square roots of the counts to approximately adjust for scale differences and thus ease comparison between observed and expected counts.
pp_check(M, resp="V", type="rootogram", nsamples=NULL, prob=0.89) + scale_x_discrete(limits=c(seq(1:9)))## Warning: Continuous limits supplied to discrete scale.
## Did you mean `limits = factor(...)` or `scale_*_continuous()`?
pp_check(M, resp="A", type="rootogram", nsamples=NULL, prob=0.89) + scale_x_discrete(limits=c(seq(1:9)))## Warning: Continuous limits supplied to discrete scale.
## Did you mean `limits = factor(...)` or `scale_*_continuous()`?
pp_check(M, resp="D", type="rootogram", nsamples=NULL, prob=0.89) + scale_x_discrete(limits=c(seq(1:9)))## Warning: Continuous limits supplied to discrete scale.
## Did you mean `limits = factor(...)` or `scale_*_continuous()`?
Finally, we can check \(\sigma\) in the posterior and compare to the prior,
post <- posterior_samples(M)
p1 <- tibble(x = seq(from = 0, to = 10, by = .5)) %>%
ggplot() +
# the prior
geom_ribbon(aes(x = x, ymin = 0, ymax = dhcauchy(x, sigma = 2)),
fill = "black", alpha = 1/10) +
# the posterior
geom_density(data = post,
aes(x = sd_ID_idx__V_Intercept),
fill = "black", alpha = 1/5, size = 0) +
xlab(expression(sigma[V])) +
theme_tufte()
# A
p2 <- tibble(x = seq(from = 0, to = 10, by = .5)) %>%
ggplot() +
# the prior
geom_ribbon(aes(x = x, ymin = 0, ymax = dhcauchy(x, sigma = 2)),
fill = "black", alpha = 1/10) +
# the posterior
geom_density(data = post,
aes(x = sd_ID_idx__A_Intercept),
fill = "black", alpha = 1/5, size = 0) +
xlab(expression(sigma[A])) +
theme_tufte()
# D
p3 <- tibble(x = seq(from = 0, to = 10, by = .5)) %>%
ggplot() +
# the prior
geom_ribbon(aes(x = x, ymin = 0, ymax = dhcauchy(x, sigma = 2)),
fill = "black", alpha = 1/10) +
# the posterior
geom_density(data = post,
aes(x = sd_ID_idx__D_Intercept),
fill = "black", alpha = 1/5, size = 0) +
xlab(expression(sigma[D])) +
theme_tufte()
p1 / p2 / p3
The dark densities are the estimated standard deviations. We can see that they differ from our priors and that they are quite precisely estimated.
6 Inference
In the manuscript we report \(95\)% credible intervals. The following parameters were significant on the \(95\)%-level:
- Outcome \(V\), covariate EXAMPLE DL
- Outcome \(V\), covariate EXAMPLE CL
- Outcome \(V\), covariate EXAMPLE BH
- Outcome \(D\), covariate EXAMPLE BL
Additionally, we wanted to examine Experience’s effect:
- Outcome \(V\), covariate EXPERIENCE
Let’s plot Experience as a conditional effect (for all three of the outcomes \(\{V,A,D\}\)), i.e., keep all covariates at their mean, or the reference level for categorical covariates, and vary only Experience (as we do in the manuscript).
conditional_effects(M, effects = "EXP", resp = "V")
conditional_effects(M, effects = "EXP", resp = "A")
conditional_effects(M, effects = "EXP", resp = "D")
We could also set some fixed values for our covariates and see what this implies. Below we define a function to do just that,
# Remember, our model had six fixed effects:
# EXAMPLE_idx, MAJOR_idx, ROLE_idx, LANG_idx, ENTITIES, and EXP
# Let's use a function to create fitted() response values for A, V,
# and D (where each consist of a Likert 1-9).
avg_resp <- function(resp, example){
if (resp != "A" & resp != "D" & resp != "V")
stop("Only outcomes A, V or D are supported.")
nd <-
tibble(EXAMPLE_idx = example,
MAJOR_idx = c("SE"), # software engineer (largest group, n=90)
EDU = 3, # Bachelor (largest group, n=85)
ROLE_idx = c("D"), # developer (largest group, n=125)
LANG_idx = c("CS"), # C# (largest group, n=100)
ENTITIES = 0, # we've scaled this already
EXP = 0 # same here
)
max_iter <- 100
resp <- resp
fitted(M,
newdata = nd,
subset = 1:max_iter,
re_formula = NA,
resp = resp,
summary = FALSE) %>%
as_tibble() %>%
gather() %>%
mutate(iter = rep(1:max_iter, times = 18)) %>%
select(iter, everything()) %>%
separate(key, into = c("intention", "rating")) %>%
mutate(intention = intention %>% as.double(),
rating = rating %>% as.double()) %>%
mutate(intention = intention -1) %>%
rename(pk = value) %>%
mutate(`pk:rating` = pk * rating) %>%
group_by(iter, intention) %>%
arrange(iter, intention, rating) %>%
mutate(probability = cumsum(pk)) %>%
filter(rating < 9) %>%
summarise(mean_rating = sum(`pk:rating`))
}Next, we use the above function. In our case, we have decided that we are interested in the largest category or median for each predictor, i.e.,
- major: SE \(n=90\)
- edu: bachelor \(n=85\)
- role: developer \(n=125\)
- lang: C# \(n=100\)
- entities: \(5\) (median)
- exp: \(8\) (median)
So, this is just another view we can take on this. Let’s analyze in particular outcome \(D\) and the artifact \(B\), Low and High.
example = c("BL", "BH")
avg_resp(resp = "D",
example = example) %>%
ggplot(aes(x = intention, y = mean_rating, group = iter)) +
geom_line(alpha = 5 / 10, color = "black") +
scale_x_continuous("outcome D", breaks = 0:1,
labels = c(example[1], example[2])) +
scale_y_continuous("", limits=c(2,6)) +
theme_hc() +
theme(plot.subtitle = element_text(size = 10),
axis.title = element_text(size = 9)
) 
What we see above are \(100\) draws from the posterior probability distribution. This is concretely what a significant effect looks like in this case.
6.1 Effect sizes
If we want to consider effects sizes then that’s very easy since we have a posterior probability distribution. In short, we simply count the number of times something is better than something else when we change the levels on factors. Let’s continue by using \(D\) and \(B\), Low/High.
# need a df and not a list
df <- as.data.frame(dat)
# set as factors
df$EXAMPLE_idx <- as.factor(df$EXAMPLE_idx)
# make sure to remove unused levels
df$EDU <- droplevels(df$EDU)
# For level BL vs BH in D
# create a representative sample
data_1 <- df[df$EXAMPLE_idx == levels(df$EXAMPLE_idx)[4], ] # BL
# do posterior predictions
PPD_1 <- posterior_predict(M, newdata = data_1)
# Now do the same for level BH
data_2 <- data_1
data_2$EXAMPLE_idx[1:nrow(data_2)] <- levels(df$EXAMPLE_idx)[3] # BH
PPD_2 <- posterior_predict(M, newdata = data_2)
# What is the difference?
PPD_diff <- PPD_2 - PPD_1
# remember, we had three outcomes D, A, and V, so store D separately
outcome_d <- PPD_diff[,,1] # D
# Simple check means of the rows
summary(rowMeans(outcome_d))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.5789 -0.1053 0.4211 0.4337 0.9474 3.0000Let’s do the same with the covariate years. First, create a function we can use, and then set two settings for years (\(5\), \(30\)), and a level \(B\) High.
comp_exp <- function(level, yrs){
if (level > 10 | level < 1 | yrs < 0)
stop("Set correct level and/or use only Z^+ for years.")
# pick out all rows where EXAMPLE_idx = level
# > levels(df$EXAMPLE_idx)
# [1] "AH" "AL" "BH" "BL" "CH" "CL" "DH" "DL" "EH" "EL"
dat1yrs <- df[df$EXAMPLE_idx == levels(df$EXAMPLE_idx)[level], ]
# Change EXP to yrs (check mean and sd on original dataset, since we
# standardized)
dat1yrs$EXP[1:nrow(dat1yrs)] <- (yrs-mean(d$EXP))/sd(d$EXP)
# predict responses based on the fitted model
posterior_predict(M, newdata = dat1yrs)
}Now we use the above function,
yrs1 <- 5
yrs2 <- 30
lvl <- 3 # BH
# BH
yrs_low <- comp_exp(lvl, yrs1)
yrs_high <- comp_exp(lvl,yrs2)
# what's the difference?
PPD_diff <- yrs_high - yrs_low
# D Below implies that yrs_low is larger by ~0.3
summary(rowMeans(PPD_diff[,,1]))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.6667 -0.9524 -0.3810 -0.3491 0.2857 3.4762As you can imagine, you can do as many of these checks as you want, and the good thing is that you usually don’t have to worry about multiple testing. Let’s focus on the Bayesian hypothesis testing approach next.
7 Hypothesis testing
Even though it might feel unnecessary to conduct hypothesis testing, many readers are comfortable with the concept and as such we show below how this can be done. Worth mentioning here is that if you conduct hypothesis testing in a Bayesian framework, then you need to have conducted prior predictive checks and sensitivity analysis beforehand. Setting good priors is key when in comes to using Bayes Factor, which is Bayesian hypothesis testing.
Below follows the code for the plots in the manuscript,
# First V
A <- hypothesis(M, "V_EXAMPLE_idxAL < 0", class = "b")
B <- hypothesis(M, "V_EXAMPLE_idxBL < V_EXAMPLE_idxBH", class = "b")
C <- hypothesis(M, "V_EXAMPLE_idxCL < V_EXAMPLE_idxCH", class = "b")
D <- hypothesis(M, "V_EXAMPLE_idxDL < V_EXAMPLE_idxDH", class = "b")
E <- hypothesis(M, "V_EXAMPLE_idxEL < V_EXAMPLE_idxEH", class = "b")
df <- data.frame(A=A$samples$H1, B=B$samples$H1, C=C$samples$H1,
D=D$samples$H1, E=E$samples$H1)
df <- gather(df, 'A', 'B', 'C', 'D', 'E', key = "Category", value = "Contrast")
ggplot(df, aes(x = Contrast,
y = fct_relevel(Category, levels = "C", "E", "B", "A", "D"))) +
geom_density_ridges(quantile_lines = TRUE, quantiles = c(0.025, 0.975),
scale = 2, rel_min_height = 0.01, show.legend = FALSE,
alpha = 0.7) +
ylab("Hypotheses") +
geom_vline(xintercept = 0, size = 0.2) +
annotate("text", x = 1.2, y = 'D',
label = deparse(bquote("Strong evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = 1.2, y = 'A',
label = deparse(bquote("Strong evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = 1.2, y = 'B',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = -2.6, y = 'E',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[0])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
annotate("text", x = -2.6, y = 'C',
label = deparse(bquote("Strong evi. for"~italic(H)[0])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
theme_tufte()## Picking joint bandwidth of 0.0751
# Next we got for A
A<- hypothesis(M, "A_EXAMPLE_idxAL < 0", class = "b")
B <- hypothesis(M, "A_EXAMPLE_idxBL < A_EXAMPLE_idxBH", class = "b")
C <- hypothesis(M, "A_EXAMPLE_idxCL < A_EXAMPLE_idxCH", class = "b")
D <- hypothesis(M, "A_EXAMPLE_idxDL < A_EXAMPLE_idxDH", class = "b")
E <- hypothesis(M, "A_EXAMPLE_idxEL < A_EXAMPLE_idxEH", class = "b")
df <- data.frame(A=A$samples$H1, B=B$samples$H1, C=C$samples$H1,
D=D$samples$H1, E=E$samples$H1)
df <- gather(df, 'A', 'B', 'C', 'D', 'E', key = "Category", value = "Contrast")
ggplot(df, aes(x = Contrast,
y = fct_relevel(Category, levels = "D", "B", "E", "A", "C"))) +
geom_density_ridges(quantile_lines = TRUE, quantiles = c(0.025, 0.975),
scale = 2, rel_min_height = 0.01, show.legend = FALSE,
alpha = 0.7) +
ylab("Hypotheses") +
geom_vline(xintercept = 0, size = 0.2) +
annotate("text", x = -2.6, y = 'D',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[0])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = -2.6, y = 'A',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = -2.6, y = 'B',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[0])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = -2.6, y = 'E',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
annotate("text", x = -2.6, y = 'C',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
theme_tufte()## Picking joint bandwidth of 0.0762
# Finally we got D
A<- hypothesis(M, "D_EXAMPLE_idxAL < 0", class = "b")
B <- hypothesis(M, "D_EXAMPLE_idxBL < D_EXAMPLE_idxBH", class = "b")
C <- hypothesis(M, "D_EXAMPLE_idxCL < D_EXAMPLE_idxCH", class = "b")
D <- hypothesis(M, "D_EXAMPLE_idxDL < D_EXAMPLE_idxDH", class = "b")
E <- hypothesis(M, "D_EXAMPLE_idxEL < D_EXAMPLE_idxEH", class = "b")
df <- data.frame(A=A$samples$H1, B=B$samples$H1, C=C$samples$H1,
D=D$samples$H1, E=E$samples$H1)
df <- gather(df, 'A', 'B', 'C', 'D', 'E', key = "Category", value = "Contrast")
ggplot(df, aes(x = Contrast,
y = fct_relevel(Category, levels = "D", "A", "E", "C", "B"))) +
geom_density_ridges(quantile_lines = TRUE, quantiles = c(0.025, 0.975),
scale = 2, rel_min_height = 0.01, show.legend = FALSE,
alpha = 0.7) +
ylab("Hypotheses") +
geom_vline(xintercept = 0, size = 0.2) +
annotate("text", x = -1.3, y = 'D',
label = deparse(bquote("Moderate evi. for"~italic(H)[0])),
parse = TRUE, hjust = "right", family = "serif", size=2) +
annotate("text", x = 1.2, y = 'A',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = 1.2, y = 'B',
label = deparse(bquote("Moderate evi. for"~italic(H)[1])),
parse = TRUE, hjust = 0, family = "serif", size=2) +
annotate("text", x = 1.2, y = 'E',
label = deparse(bquote("Anecdotal evi. for"~italic(H)[1])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
annotate("text", x = 1.2, y = 'C',
label = deparse(bquote("Moderate evi. for"~italic(H)[1])),
parse = TRUE, hjust = "left", family = "serif", size=2) +
theme_tufte()## Picking joint bandwidth of 0.0755
8 Computational environment
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))## CXX14FLAGS=-O3 -mtune=native
## CXX14FLAGS += -arch x86_64 -ftemplate-depth-256 -Wno-unused-local-typedefs -Wno-unused-function
## FC = /usr/local/bin/gfortran
## FLIBS = -L/usr/local/lib/gcc/10 -L/usr/local/lib -lgfortran -lquadmath -lmsessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] forcats_0.5.0 ggridges_0.5.2 extraDistr_1.9.1
## [4] tidyr_1.1.2 patchwork_1.0.1 ggthemes_4.2.0
## [7] dplyr_1.0.2 rethinking_2.13 rstan_2.21.2
## [10] ggplot2_3.3.2 StanHeaders_2.21.0-6 brms_2.13.9
## [13] Rcpp_1.0.5
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-149 matrixStats_0.56.0 xts_0.12.1
## [4] threejs_0.3.3 tools_4.0.2 backports_1.1.9
## [7] R6_2.4.1 DT_0.15 colorspace_1.4-1
## [10] withr_2.2.0 tidyselect_1.1.0 gridExtra_2.3
## [13] prettyunits_1.1.1 processx_3.4.4 Brobdingnag_1.2-6
## [16] curl_4.3 compiler_4.0.2 cli_2.0.2
## [19] shinyjs_2.0.0 labeling_0.3 colourpicker_1.1.0
## [22] bookdown_0.20 scales_1.1.1 dygraphs_1.1.1.6
## [25] mvtnorm_1.1-1 callr_3.4.4 stringr_1.4.0
## [28] digest_0.6.25 rmarkdown_2.3 base64enc_0.1-3
## [31] pkgconfig_2.0.3 htmltools_0.5.0 fastmap_1.0.1
## [34] htmlwidgets_1.5.1 rlang_0.4.7 shiny_1.5.0
## [37] farver_2.0.3 shape_1.4.5 generics_0.0.2
## [40] zoo_1.8-8 jsonlite_1.7.1 crosstalk_1.1.0.1
## [43] gtools_3.8.2 inline_0.3.16 magrittr_1.5
## [46] loo_2.3.1 bayesplot_1.7.2 Matrix_1.2-18
## [49] munsell_0.5.0 fansi_0.4.1 abind_1.4-5
## [52] lifecycle_0.2.0 stringi_1.5.3 yaml_2.2.1
## [55] MASS_7.3-53 pkgbuild_1.1.0 plyr_1.8.6
## [58] grid_4.0.2 promises_1.1.1 crayon_1.3.4
## [61] miniUI_0.1.1.1 lattice_0.20-41 knitr_1.29
## [64] ps_1.3.4 pillar_1.4.6 igraph_1.2.5
## [67] markdown_1.1 shinystan_2.5.0 reshape2_1.4.4
## [70] codetools_0.2-16 stats4_4.0.2 rstantools_2.1.1
## [73] glue_1.4.2 evaluate_0.14 V8_3.2.0
## [76] RcppParallel_5.0.2 vctrs_0.3.4 httpuv_1.5.4
## [79] gtable_0.3.0 purrr_0.3.4 assertthat_0.2.1
## [82] xfun_0.17 mime_0.9 xtable_1.8-4
## [85] coda_0.19-3 later_1.1.0.1 rsconnect_0.8.16
## [88] tibble_3.0.3 shinythemes_1.1.2 cmdstanr_0.1.2
## [91] ellipsis_0.3.1 bridgesampling_1.0-0References
Blei, D. M., A. Y. Ng, and M. I. Jordan. 2003. “Latent Dirichlet Allocation.” Journal of Machine Learning Research 3 (Jan): 993–1022.
Vehtari, A., A. Gelman, and J. Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.