Learning to write syntheses - Walkthrough

In this document we document the different steps in the analyses for the study “Learning to write syntheses”.

The following steps are documented:

• Setting the stage for the analyses;
• Methodology;
• Analyses and results

The datasets, Stan code and results (draws from the posterior) and the R scripts are available on OSF (Open Science Framework). To access these use the following link: https://osf.io/3sb7m/?view_only=b89adb32abaa46c59676cc7c6a47a62f.

In this walkthrough some R-code is shared as well. To be able to reproduce this file and the code, it is important to place this file in an RStudio project with the same folder structure as used in our project. The folder structure we used is the follwing:

- Data
- Figs
- Output
    - Baseline_models
    - Draws_Correlations
    - Fitted_Models
    - Loo
-  R_Scripts
-  Stan

Setting the stage for the analyses

Before we can perform the actual analyses, we have to set some preparatory steps. In what follows we share the code that is needed to load the necessary R packages, read the data and prepare the data for the statistical analyses.

Necessary packages

For this study we made use of the following packages:

• tidyverse (Wickham et al. 2019) which contains a set of packages that facilitate the tidy way of working in R
• foreign (2022) is a package that enables us to import SPSS datafiles
• readxl (Wickham and Bryan 2022) is a package that we use to import MS Excel worksheets
• here (Müller 2020) to create reproducible code that is not dependent on local paths of computers
• brms (Bürkner 2017) is a package to estimate the Bayesion models in R making use of Stanat the background
• ggdist (Kay 2022) is a package that has some great implementations of visualizations that support the interpretation of the results
• patchwork (Pedersen 2022) to combine plots to a single visualization
• kableExtra (Zhu 2021) to print the tables in a formatted way in this document
• bayestestR (Makowski, Ben-Shachar, and Lüdecke 2019) to access a set of functions to summarize the posterior draws

Code

library(tidyverse)
library(foreign)
library(readxl)
library(here)
library(brms)
library(bayestestR)
library(ggdist)
library(patchwork)
library(kableExtra)

Set a plotting theme

To make all our visualizations uniform we will set a theme for all plots generated with ggplot2. For consistency in the colours used in visualizations we create two named vectors that can be used to define colours for the three (or two) groups in the plot: the two conditions and the baseline. These vectors are called Conditions_colours (3 colour) and Conditions_colours2 (2 colours: baseline vs intervention).

Code

theme_set(theme_linedraw() +
            theme(panel.grid = element_blank()))

Conditions_colours <- c(
  "Feedback only" =  "#f1b53d",
  "Feedback + Observation" = "#82A1AD",
  "Baseline" = "lightgrey"
)

Conditions_colours2 <- c(
  "Intervention" =  "grey40",
  "Baseline" = "lightgrey"
)

Load the different datasets

The first dataset is the data with the process data of the writing processes of students in the experimental groups. The pseudonomized version of this dataset is stored in a file called ProcesData_Experiment.RDS in the Data folder.

Next we can also load the data file containing the information on text quality (TextQuality_Experiment.RDS) stored in the Data folder.

Code

ProcesData <- readRDS(
  file = here(
    "Data",
    "ProcesData_Experiment.RDS"
  )
)

TextQuality <- readRDS(
  file = here(
    "Data",
    "TextQuality_Experiment.RDS"
  )
)

Methodology

Bayesian Analyses

To estimate the impact of the intervention and test whether there is an effect of the intervention we will make use of Bayesian modelling.

For the analyses we will rely on three alternative models, that are an implementation of a multivariate model, assuming that the scores for each session are three separate dependent variables that will be modeled simultaneously.

A first model (Model1) is an unconditional model not taking into account the possible effects of the intervention. Formally this model can be written as,

\[\begin{equation} \begin{aligned} \begin{bmatrix} \text{Session1}_{i}\\ \text{Session2}_{i}\\ \text{Session3}_{i} \end{bmatrix} & \sim \textbf{ MVNormal} \left( \begin{bmatrix} \mu_{1} \\ \mu_{2} \\ \mu_{3} \end{bmatrix} , \Sigma_1 \right) \\ \mu_{1} & = \beta_{01} \\ \mu_{2} & = \beta_{02} \\ \mu_{3} & = \beta_{03} \\ \Sigma_1 & = \textbf{S}_1\textbf{R}_1\textbf{S}_1 \\ \textbf{S}_1 & = \begin{bmatrix} \sigma_{e_{1i}} & 0 & 0 \\ 0 & \sigma_{e_{2i}} & 0 \\ 0 & 0 & \sigma_{e_{3i}} \end{bmatrix}\\ \textbf{R}_1 & = \begin{bmatrix} 1 & \rho_{e_{1i}e_{2i}} & \rho_{e_{1i}e_{3i}} \\ \rho_{e_{2i}e_{1i}} & 1 & \rho_{e_{2i}e_{3i}} \\ \rho_{e_{3i}e_{1i}} & \rho_{e_{3i}e_{2i}} & 1 \end{bmatrix}\\ \end{aligned} \end{equation}\]

with \(\beta_{01} \rightarrow \beta_{03}\) being the expected mean values for each session and allowing the residuals \(e_{1i} \rightarrow e_{3i}\) to be correlated. In other words, this model assumes that there are 3 potentially different mean values for the dependent variable, one mean score per session. Further, no differences between both conditions are expected in this first model.

Now, of course our hypothesis is that students in the Feedback and Instruction condition are influenced by the intervention and will score different than the students in the Feedback only condition. This hypothesis is reflected in an expected difference between both conditions when we look at the scores for variables both at Session2 and Session3. Therefore, we introduce the effect of the dummy variable Condition which has value 1 for students in the Feedback and Instruction condition and value 0 for students in the Feedback only condition. The dummy variable is added to the model as a predictor for scores on Session2 and Session3 in Model2. Specifically, in Model 2 we assume that the effect of Condition is identical at both sessions. In other words, there is a single parameter expressing the effect of Condition (\(\beta_{1}\) in the following model), which is a parsimonious model.

\[\begin{equation} \begin{aligned} \begin{bmatrix} \text{Session1}_{i}\\ \text{Session2}_{i}\\ \text{Session3}_{i} \end{bmatrix} & \sim \textbf{ MVNormal} \left( \begin{bmatrix} \mu_{1} \\ \mu_{2} \\ \mu_{3} \end{bmatrix} , \Sigma_1 \right) \\ \mu_{1} & = \beta_{01} + \\ \mu_{2} & = \beta_{02} + \beta_{1} * \text{Condition} \\ \mu_{3} & = \beta_{03} + \beta_{1} * \text{Condition}\\ \Sigma_1 & = \textbf{S}_1\textbf{R}_1\textbf{S}_1 \\ \textbf{S}_1 & = \begin{bmatrix} \sigma_{e_{1i}} & 0 & 0 \\ 0 & \sigma_{e_{2i}} & 0 \\ 0 & 0 & \sigma_{e_{3i}} \end{bmatrix}\\ \textbf{R}_1 & = \begin{bmatrix} 1 & \rho_{e_{1i}e_{2i}} & \rho_{e_{1i}e_{3i}} \\ \rho_{e_{2i}e_{1i}} & 1 & \rho_{e_{2i}e_{3i}} \\ \rho_{e_{3i}e_{1i}} & \rho_{e_{3i}e_{2i}} & 1 \end{bmatrix}\\ \end{aligned} \end{equation}\]

A third model releases this assumption, estimating an effect of Condition that is not similar at Session 2 and Session 3, expressed in the following model by two separate parameters (\(\beta_{12}\) and \(\beta_{13}\)).

\[\begin{equation} \begin{aligned} \begin{bmatrix} \text{Session1}_{i}\\ \text{Session2}_{i}\\ \text{Session3}_{i} \end{bmatrix} & \sim \textbf{ MVNormal} \left( \begin{bmatrix} \mu_{1} \\ \mu_{2} \\ \mu_{3} \end{bmatrix} , \Sigma_1 \right) \\ \mu_{1} & = \beta_{01} + \\ \mu_{2} & = \beta_{02} + \beta_{12} * \text{Condition} \\ \mu_{3} & = \beta_{03} + \beta_{13} * \text{Condition}\\ \Sigma_1 & = \textbf{S}_1\textbf{R}_1\textbf{S}_1 \\ \textbf{S}_1 & = \begin{bmatrix} \sigma_{e_{1i}} & 0 & 0 \\ 0 & \sigma_{e_{2i}} & 0 \\ 0 & 0 & \sigma_{e_{3i}} \end{bmatrix}\\ \textbf{R}_1 & = \begin{bmatrix} 1 & \rho_{e_{1i}e_{2i}} & \rho_{e_{1i}e_{3i}} \\ \rho_{e_{2i}e_{1i}} & 1 & \rho_{e_{2i}e_{3i}} \\ \rho_{e_{3i}e_{1i}} & \rho_{e_{3i}e_{2i}} & 1 \end{bmatrix}\\ \end{aligned} \end{equation}\]

Results

Correlations

A first exploration we can make is inspecting the correlations between the different behavioural measures and text quality. We do this for each session separately. The correlations are modeled in a bayesian framework (code available on OSF and the rationale is explained in the paper). As a result we can plot draws from the posterior to envisage the uncertainties for the correlation parameters. The following code creates a figure that plots the 89% credible interval for each of the correlation estimates. The dotted lines indicate the rule of thumb of Cohen to decide on the strength of correlations.

Code

# Load the draws from the posterior probability distributions for the correlations

Estimates_Session1 <- readRDS(
  here("Output",
       "Draws_Correlations",
       "Correlations_Session1.RDS")
)


Estimates_Session2 <- readRDS(
  here("Output",
       "Draws_Correlations",
       "Correlations_Session2.RDS")
)

Estimates_Session3 <- readRDS(
  here("Output",
       "Draws_Correlations",
       "Correlations_Session3.RDS")
)


# Create the plots for each session

S1 <- Estimates_Session1 %>%
  select(starts_with("rescor_"))  %>%
  pivot_longer(everything()) %>%
  mutate(
    name = case_when(
      name == "rescor__Graderescaled__AWTint1" ~ "TQ_AWT",
      name == "rescor__Graderescaled__Speed" ~ "TQ_KPM",
      name == "rescor__Graderescaled__TimeinSources" ~ "TQ_TIS",
      name == "rescor__Graderescaled__SourcesSwitch" ~ "TQ_SBS",
      name == "rescor__Graderescaled__SourceTextSwitch" ~ "TQ_SST",
      name == "rescor__AWTint1__Speed" ~ "AWT_KPM",
      name == "rescor__AWTint1__TimeinSources" ~ "AWT_TIS",
      name == "rescor__Speed__TimeinSources" ~ "KPM_TIS",
      name == "rescor__AWTint1__SourcesSwitch" ~ "AWT_SBS",
      name == "rescor__Speed__SourcesSwitch" ~ "KPM_SBS",
      name == "rescor__TimeinSources__SourcesSwitch" ~ "TIS_SBS",
      name == "rescor__AWTint1__SourceTextSwitch" ~ "AWT_SST",
      name == "rescor__Speed__SourceTextSwitch" ~ "KPM_SST",
      name == "rescor__TimeinSources__SourceTextSwitch" ~ "TIS_SST",
      name == "rescor__SourcesSwitch__SourceTextSwitch" ~ "SBS_SST"
      )
    ) %>%
  mutate(
    Order = case_when(
      name == "KPM_TIS" ~ 1,
      name == "AWT_TIS" ~ 2,
      name == "KPM_SBS" ~ 3,
      name == "TIS_SST" ~ 4,
      name == "AWT_SBS" ~ 5,
      name == "TQ_TIS" ~  6,
      name == "SBS_SST" ~ 7,
      name == "TQ_SBS" ~  8,
      name == "TQ_AWT" ~  9,
      name == "TQ_KPM" ~ 10,
      name == "TQ_SST" ~ 11,
      name == "TIS_SBS" ~ 12,
      name == "KPM_SST" ~ 13,
      name == "AWT_SST" ~ 14,
      name == "AWT_KPM" ~15 
    )
  ) %>%
  ggplot(
    aes(x = value, y = reorder(name, Order))
  ) + 
  stat_pointinterval(point_interval = mean_qi, .width = .89, 
                     size = 1/6) +
  geom_vline(
    xintercept = 0,  
    size = 0.4) +
  geom_vline(
    xintercept = c(-0.8, -0.5, -0.2, 0.2, 0.5, 0.8),  
    size = 0.2,
    lty = 2) + 
  scale_x_continuous(
    limits = c(-1,1),
    breaks = seq(-1,1,0.2)
  ) +
  labs(
    title = "Session 1",
    y = "",
    x = "correlation"
  )

S2 <- Estimates_Session2 %>%
  select(starts_with("rescor_"))  %>%
  pivot_longer(everything()) %>%
  mutate(
    name = case_when(
      name == "rescor__Graderescaled__AWTint1" ~ "TQ_AWT",
      name == "rescor__Graderescaled__Speed" ~ "TQ_KPM",
      name == "rescor__Graderescaled__TimeinSources" ~ "TQ_TIS",
      name == "rescor__Graderescaled__SourcesSwitch" ~ "TQ_SBS",
      name == "rescor__Graderescaled__SourceTextSwitch" ~ "TQ_SST",
      name == "rescor__AWTint1__Speed" ~ "AWT_KPM",
      name == "rescor__AWTint1__TimeinSources" ~ "AWT_TIS",
      name == "rescor__Speed__TimeinSources" ~ "KPM_TIS",
      name == "rescor__AWTint1__SourcesSwitch" ~ "AWT_SBS",
      name == "rescor__Speed__SourcesSwitch" ~ "KPM_SBS",
      name == "rescor__TimeinSources__SourcesSwitch" ~ "TIS_SBS",
      name == "rescor__AWTint1__SourceTextSwitch" ~ "AWT_SST",
      name == "rescor__Speed__SourceTextSwitch" ~ "KPM_SST",
      name == "rescor__TimeinSources__SourceTextSwitch" ~ "TIS_SST",
      name == "rescor__SourcesSwitch__SourceTextSwitch" ~ "SBS_SST"
      )
    ) %>%
  mutate(
    Order = case_when(
      name == "KPM_TIS" ~ 1,
      name == "AWT_TIS" ~ 2,
      name == "KPM_SBS" ~ 3,
      name == "TIS_SST" ~ 4,
      name == "AWT_SBS" ~ 5,
      name == "TQ_TIS" ~  6,
      name == "SBS_SST" ~ 7,
      name == "TQ_SBS" ~  8,
      name == "TQ_AWT" ~  9,
      name == "TQ_KPM" ~ 10,
      name == "TQ_SST" ~ 11,
      name == "TIS_SBS" ~ 12,
      name == "KPM_SST" ~ 13,
      name == "AWT_SST" ~ 14,
      name == "AWT_KPM" ~15 
    )
  ) %>%
  ggplot(
    aes(x = value, y = reorder(name, Order))
  ) + 
  stat_pointinterval(point_interval = mean_qi, .width = .89, 
                     size = 1/6) +
  geom_vline(
    xintercept = 0,  
    size = 0.4) +
  geom_vline(
    xintercept = c(-0.8, -0.5, -0.2, 0.2, 0.5, 0.8),  
    size = 0.2,
    lty = 2) + 
  scale_x_continuous(
    limits = c(-1,1),
    breaks = seq(-1,1,0.2)
  ) +
  labs(
    title = "Session 2",
    y = "",
    x = "correlation"
  )

S3 <- Estimates_Session3 %>%
  select(starts_with("rescor_"))  %>%
  pivot_longer(everything()) %>%
  mutate(
    name = case_when(
      name == "rescor__Graderescaled__AWTint1" ~ "TQ_AWT",
      name == "rescor__Graderescaled__Speed" ~ "TQ_KPM",
      name == "rescor__Graderescaled__TimeinSources" ~ "TQ_TIS",
      name == "rescor__Graderescaled__SourcesSwitch" ~ "TQ_SBS",
      name == "rescor__Graderescaled__SourceTextSwitch" ~ "TQ_SST",
      name == "rescor__AWTint1__Speed" ~ "AWT_KPM",
      name == "rescor__AWTint1__TimeinSources" ~ "AWT_TIS",
      name == "rescor__Speed__TimeinSources" ~ "KPM_TIS",
      name == "rescor__AWTint1__SourcesSwitch" ~ "AWT_SBS",
      name == "rescor__Speed__SourcesSwitch" ~ "KPM_SBS",
      name == "rescor__TimeinSources__SourcesSwitch" ~ "TIS_SBS",
      name == "rescor__AWTint1__SourceTextSwitch" ~ "AWT_SST",
      name == "rescor__Speed__SourceTextSwitch" ~ "KPM_SST",
      name == "rescor__TimeinSources__SourceTextSwitch" ~ "TIS_SST",
      name == "rescor__SourcesSwitch__SourceTextSwitch" ~ "SBS_SST"
      )
    ) %>%
  mutate(
    Order = case_when(
      name == "KPM_TIS" ~ 1,
      name == "AWT_TIS" ~ 2,
      name == "KPM_SBS" ~ 3,
      name == "TIS_SST" ~ 4,
      name == "AWT_SBS" ~ 5,
      name == "TQ_TIS" ~  6,
      name == "SBS_SST" ~ 7,
      name == "TQ_SBS" ~  8,
      name == "TQ_AWT" ~  9,
      name == "TQ_KPM" ~ 10,
      name == "TQ_SST" ~ 11,
      name == "TIS_SBS" ~ 12,
      name == "KPM_SST" ~ 13,
      name == "AWT_SST" ~ 14,
      name == "AWT_KPM" ~15 
    )
  ) %>%
  ggplot(
    aes(x = value, y = reorder(name, Order))
  ) + 
  stat_pointinterval(point_interval = mean_qi, .width = .89, 
                     size = 1/6) +
  geom_vline(
    xintercept = 0,  
    size = 0.4) +
  geom_vline(
    xintercept = c(-0.8, -0.5, -0.2, 0.2, 0.5, 0.8),  
    size = 0.2,
    lty = 2) + 
  scale_x_continuous(
    limits = c(-1,1),
    breaks = seq(-1,1,0.2)
  ) +
  labs(
    title = "Session 3",
    y = "",
    x = "correlation"
  )

# Create an object to add to the figures with a legend on the variable names

Labels_var <-
  tibble(
    Abb = c("TQ", "AWT", "KPM", "SST", "TIS", "SBS"),
    Variable = c("Text quality", 
                 "Active writing time", 
                 "Speed",
                 "Switches sources-synthesis ",
                 "Source time",
                 "Switches sources")
  ) %>%
  column_to_rownames(var = "Abb") 

## Combine the plots in one plot with patchwork code

tt <- gridExtra::ttheme_minimal(
  base_size = 10,
  core = list(
    fg_params=list(hjust=0, x = 0.01)
    )
)

S1 + S2 + S3 + gridExtra::tableGrob(Labels_var, theme = tt, cols = NULL) + 
  plot_layout(ncol = 2)

## Store the plot in the Figs folder for publication

ggsave(
  device = "jpeg",
  file = here("Figs", "Fig1_Correlations.jpg"), 
  dpi = 300
  )

Figure 1: Posterior distribution for the correlations between the process variables and text quality at each session.

Effects of intervention on Text Quality

Before we perform the analyses we plot the data to inspect the data visually.

Code

TextQuality %>%
  filter(!is.na(Condition)) %>%
  ggplot(aes(x = Session, y = Grade_rescaled, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Text Quality Scores") 

Figure 2: Text Quality in the first interval of the writing process for each session and both conditions

First we start by comparing the model fit of the three alternative models, comparing the expected log of the predicted density (\(\widehat{elpd}\)) for each model (which is a measure that can be used as a relative measure of out-of-sample predictive performance of a model).

As the results of the leave-one-out cross-validation are stored in files (see OSF), we can load these files and perform a model comparison.

The following table shows the difference in \(\widehat{elpd}\) (elpd_diff) between the best fitting model (model printed above in the list) and the other two models, the standard error for this difference (se_diff), the \(\widehat{elpd}\) value of the models (elpd_loo), the standard error for the \(\widehat{elpd}\) (se_elpd_loo), the “effective number of parameters” of the model (p_loo), the standard error for the effective number of parameters (se_p_loo), the LOO information criteria (looic) and the standard error for the LOO information criteria (se_looic). A model with the highest \(\widehat{elpd}\) actually is best at generalizing to new data. The LOO information criteria are nothing more than -2 times the \(\widehat{elpd}\), so that a lower looic refers to a better fit. Important to note is that these fit indices do not have any absolute value (the numbers themselves have no straightforward interpretation) but rather a relative value (you can use them to compare models). The standard error of the difference in \(\widehat{elpd}\) between two models is only reliably estimated when the difference in \(\widehat{elpd}\) is larger than 4 and the number of observations in the data is larger than 100 (Gelman, Hill, and Vehtari 2021, 178).

Code

loo_M1_TQ <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_TQ.RDS"
             )
     )

loo_M2_TQ <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_TQ.RDS"
             )
     )

loo_M3_TQ <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_TQ.RDS"
             )
     )

loo_comparison_TQ <- loo_compare(loo_M1_TQ, loo_M2_TQ, loo_M3_TQ)

kable(
  loo_comparison_TQ, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on Text Quality"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Text Quality

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model1	0.00	0.00	-2213.34	15.49	9.03	1.08	4426.67	30.98
model2	-1.08	0.26	-2214.42	15.49	10.13	1.17	4428.83	30.98
model3	-1.90	0.58	-2215.23	15.43	11.01	1.22	4430.47	30.85

The model comparison shows that the simplest model (Model 1 - with no effect of condition) fits best to the data (see lowest looic). This leads to the conclusion that we do not have convincing evidence to state that giving only feedback or accompanying feedback with instructions would result in different text quality scores. That doesn’t mean that both interventions might have had an effect on text quality. This is further explored in what follows.

The following table summarizes the information on the parameter estimates of model 1. From this table we can see that the text quality of texts written at the third session, is higher than the text quality of texts written at session one and two. The 89% most probable values for the difference between session three and one are between 13 and 28 points, while the 89% most probable estimates of the difference between session three and two are situated between 5 and 19 points. The posterior probability difference in text quality between session two and one has a median of 8 with a 89% credible interval ranging between 1 and 15, leading to the conclusion that there is also an improvement in text quality between session one and two (also shown in the fact that 97% of the posterior distribution is situated above 0, see column pd).

Code

Model1_TQ <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model1_TQ.RDS"
                       )
                     )


Draws_TQ_M1 <-  as_draws_df(Model1_TQ$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1 = b_TQ2_Intercept - b_TQ1_Intercept,
    diff_M3_M1 = b_TQ3_Intercept - b_TQ1_Intercept,
    diff_M3_M2 = b_TQ3_Intercept - b_TQ2_Intercept
  )


P <- describe_posterior(
  Draws_TQ_M1,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  )  %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_TQ1_Intercept" ~ "Intercept Session1",
      Parameter == "b_TQ2_Intercept" ~ "Intercept Session2",
      Parameter == "b_TQ3_Intercept" ~ "Intercept Session3",
      Parameter == "diff_M2_M1" ~ "Contrast: Session 2 - Session 1",
      Parameter == "diff_M3_M1" ~ "Contrast: Session 3 - Session 1",
      Parameter == "diff_M3_M2" ~ "Contrast: Session 3 - Session 2"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 1 (Text Quality)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 1 (Text Quality)

Parameter	Median	CI_low	CI_high	pd
Intercept Session1	102.654	97.059	108.319	1.000
Intercept Session2	110.688	105.412	116.038	1.000
Intercept Session3	122.947	117.651	128.222	1.000
Contrast: Session 2 - Session 1	8.045	1.292	15.082	0.969
Contrast: Session 3 - Session 1	20.337	13.061	27.686	1.000
Contrast: Session 3 - Session 2	12.273	5.358	19.178	0.998
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

To illustrate that there are no differences between both experimental conditions, we can plot the results of the posterior distribution based on model 2 (see next figure). This visualization shows a big overlap in the posteriors for both conditions at session two and session three.

Code

Model2_TQ <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model2_TQ.RDS"
                       )
                     )

as_draws_df(Model2_TQ$draws()) %>%
  select(b_TQ1_Intercept, 
         b_TQ2_Intercept, 
         b_TQ3_Intercept,
         `b_TQ[1]`) %>%
  mutate(Occ1_C0 = b_TQ1_Intercept,
         Occ2_C0 = b_TQ2_Intercept,
         Occ3_C0 = b_TQ3_Intercept,
         Occ1_C1 = b_TQ1_Intercept,
         Occ2_C1 = b_TQ2_Intercept + `b_TQ[1]`,
         Occ3_C1 = b_TQ3_Intercept + `b_TQ[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only",
                            "Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
  scale_color_manual(
    values = Conditions_colours[c(1,2)]) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 3: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Text Quality) model 2

Next, we want to test whether the effects of the interventions (as shown by an increase in text quality at session two and three) can be explained by the fact that different tasks are used at each session. Therefore, we make a comparison with information from the larger baseline study. In the following visualization we plot the posterior probability distribution of the intervention effect against the posterior probability distribution for the average text quality of similar tasks at the baseline study.

Code

Draws_baseline_TQ <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_TQ.RDS"
  )
)

Draws_baseline_TQ <- Draws_baseline_TQ %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_TQ <- as_draws_df(
  Model1_TQ$draws()
)

Draw_Exp_TQ %>%
  cbind(Draws_baseline_TQ$b_Sess1) %>%
  cbind(Draws_baseline_TQ$b_Sess2) %>%
  cbind(Draws_baseline_TQ$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_TQ$b_Sess1`,
    Baseline2 = `Draws_baseline_TQ$b_Sess2`,
    Baseline3 = `Draws_baseline_TQ$b_Sess3`,
  ) %>%
  select(
    b_TQ1_Intercept, 
    b_TQ2_Intercept, 
    b_TQ3_Intercept,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_TQ1_Intercept,
    Occ2_C0 = b_TQ2_Intercept,
    Occ3_C0 = b_TQ3_Intercept,
    Occ1_C1 = b_TQ1_Intercept,
    Occ2_C1 = b_TQ2_Intercept,
    Occ3_C1 = b_TQ3_Intercept,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Intervention",
      grepl("C1", name) ~ "Intervention",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  ylab("") + 
  xlab("Text Quality") +
  scale_fill_manual(
    values = Conditions_colours2
  ) +
  labs(fill = "Condition: ") 

ggsave(
  device = "jpeg",
  file = here("Figs", "Fig2_Posterior_TQ.jpg"),
  dpi = 300
  )

Figure 4: Posterior probability distributions for the average Text Quality at each session, for the two conditions and the baseline students.

We can also calculate whether the differences between the students in the intervention and the national baseline increases in each session. This is done in the following table.

Code

Differences_Baseline_Intervention <- 
  Draw_Exp_TQ %>%
  cbind(Draws_baseline_TQ$b_Sess1) %>%
  cbind(Draws_baseline_TQ$b_Sess2) %>%
  cbind(Draws_baseline_TQ$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_TQ$b_Sess1`,
    Baseline2 = `Draws_baseline_TQ$b_Sess2`,
    Baseline3 = `Draws_baseline_TQ$b_Sess3`,
  ) %>%
  select(
    b_TQ1_Intercept, 
    b_TQ2_Intercept, 
    b_TQ3_Intercept,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    diff_S1 = b_TQ1_Intercept - Baseline1,
    diff_S2 = b_TQ2_Intercept - Baseline2,
    diff_S3 = b_TQ3_Intercept - Baseline3
  ) %>%
  mutate(
    diff_S2_S1 = diff_S2 - diff_S1,
    diff_S3_S1 = diff_S3 - diff_S1,
    diff_S3_S2 = diff_S3 - diff_S2
  ) %>%
  select(
    starts_with("dif")
  )

P <- describe_posterior(
  Differences_Baseline_Intervention,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  )  %>%
  mutate(
    Parameter = case_when(
      Parameter == "diff_S1" ~ "diff. Int. - Base. S1",
      Parameter == "diff_S2" ~ "diff. Int. - Base. S2",
      Parameter == "diff_S3" ~ "diff. Int. - Base. S3",
      Parameter == "diff_S2_S1" ~ "Increase in diff. S2 vs S1",
      Parameter == "diff_S3_S1" ~ "Increase in diff. S3 vs S1",
      Parameter == "diff_S3_S2" ~ "Increase in diff. S3 vs S2"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Posterior probability for the differences between baseline and intervention students",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Posterior probability for the differences between baseline and intervention students

Parameter	Median	CI_low	CI_high	pd
diff. Int. - Base. S1	18.141	12.123	24.510	1.000
diff. Int. - Base. S2	23.942	18.049	29.698	1.000
diff. Int. - Base. S3	36.943	31.108	42.657	1.000
Increase in diff. S2 vs S1	5.800	-1.588	13.243	0.893
Increase in diff. S3 vs S1	18.769	11.234	26.672	1.000
Increase in diff. S3 vs S2	13.001	5.635	20.075	0.998
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

Effects of intervention on process data

We will document the analyses for the variable Active Writing Time more detailed as an example. The other variables are analysed in a similar way. At the end of this section we give an overview of the results for all process variables in an integrated figure (also used in the manuscript itself).

Active Writing Time

Before we perform the analyses we plot the data to inspect the data visually.

Code

ProcesData %>%
    mutate(
    To_remove = case_when(
      Product_ID == "B25" & Session == "1" & AWT_int1 == 0 ~ 1
    )
  ) %>%
  filter(!is.na(Condition) & is.na(To_remove)) %>%
  ggplot(aes(x = Session, y = AWT_int1, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Active Writing Time in Interval 1") 

Figure 5: Active Writing Time in the first interval of the writing process for each session and both conditions

The next step is estimating both models for this variable. This is done by the following code where we, at the one hand build a proper dataset and at the other hand use the Stan code file for the models. These Stan models are build, making use of the package brms to get a basic model expression in Stan code. This code is edited at two places. First, in all models we added code in the generated quantities block calculating the log-likelihood information needed to apply a model comparison based on a leave-one-out predictive performance approach (aka loo fit of the model). Second, we tweaked the Stan code of the second model to generate a model that has only one parameter for the effect of the variable Condition.

Note

To estimate the model in brms we first have to rearrange the dataset to a wide format with 3 columns for the scores of our variable, one column per session. This is the first piece of code in the following code chunk.

Code

ActiveWritingTime <- ProcesData %>%
  mutate(
    To_remove = case_when(
      Product_ID == "B25" & Session == "1" & AWT_int1 == 0 ~ 1,
    )
  ) %>%
  filter(!is.na(Conditie) & is.na(To_remove)) %>%
  select(Product_ID, Session, Condition, AWT_int1) %>%
  rename(AWT = AWT_int1) %>%
  mutate(ID = row_number()) %>%
  pivot_wider(id_cols = c(Product_ID, Condition), names_from = Session, values_from = AWT, names_prefix = "AWT.") %>% unchop(everything())

To estimate the models we can use the following code (also shared in the OSF repository).

To avoid long running times to create this document due to re-running the models each time we render this document, we stored the model results in separate files in the Output folder. These files will be loaded in the code.

Code

## MODEL 1: NO EFFECT OF INTERENTION

# Create the stan code via brms
# generates Stan code that cen be saved in a separate file
# this part of the code is de-activated as the model specification
# for stan is already saved

# make_stancode(
#  bf(AWT.1 ~ 1) +
#    bf(AWT.2 ~ 1) +
#    bf(AWT.3 ~ 1) + 
#    set_rescor(rescor = TRUE),
#  data = ActiveWritingTime
#  )

# Create a proper data list to match the model

M1_AW_data <- 
  make_standata(
    bf(AWT.1 ~ 1) +
      bf(AWT.2 ~ 1) +
      bf(AWT.3 ~ 1) + 
      set_rescor(rescor = TRUE),
    data = ActiveWritingTime,
  )


# Activate the model based on the stan model file

Mod1_AW <- cmdstan_model(here("Stan", "M1_NoEffect.stan"))

# fit the model with the data

fit_Mod1_AW <- Mod1_AW$sample(
  data = M1_AW_data,
  seed = 2021,
  chains = 6,
  iter_sampling = 5000,
  iter_warmup = 1500,
  parallel_chains = 6,
  refresh = 500
)

# Save the model results

Model1_AW <- fit_Mod1_AW$save_object(file = here("Output", "Fitted_Models","Model1_AW.RDS"))

# Perform loo validation

loo_M1_AW <- fit_Mod1_AW$loo()

# Save the loo validation result

saveRDS(loo_M1_AW, file = here("Output" , "Loo", "loo_M1_AW.RDS"))


## MODEL 2: ONE EFFECT OF INTERVENTION

# Create the stan code via brms
# generates Stan code that can be saved in a separate file
# in this file some adaptations have to be made to 
# model 1 effect of INTERVENTION on both sessions (2 and 3)
# this part of the code is de-activated as the model specification
# for stan is already saved

#make_stancode(
#  bf(AWT.1 ~ 1) +
#    bf(AWT.2 ~ 1 + Condition) +
#    bf(AWT.3 ~ 1 + Condition) + 
#    set_rescor(rescor = TRUE),
#  data = ActiveWritingTime,
# )

# Create a proper data list to match the model

M2_AW_data <- 
  make_standata(
    bf(AWT.1 ~ 1) +
      bf(AWT.2 ~ 1 + Condition) +
      bf(AWT.3 ~ 1 + Condition) + 
      set_rescor(rescor = TRUE),
    data = ActiveWritingTime,
  )

# Activate the model based on the stan model file

Mod2_AW <- cmdstan_model(here("Stan", "M2_OneEffect.stan"))

# fit the model with the data

fit_Mod2_AW <- Mod2_AW$sample(
  data = M2_AW_data,
  seed = 2021,
  chains = 6,
  iter_sampling = 5000,
  iter_warmup = 1500,
  parallel_chains = 6,
  refresh = 500
)


# Save the model results

Model2_AW <- fit_Mod2_AW$save_object(
  file = here("Output", 
              "Fitted_Models",
              "Model2_AW.RDS")
  )

# Perform loo validation

loo_M2_AW <- fit_Mod2_AW$loo()

# Save the loo validation result

saveRDS(loo_M2_AW, 
     file = here("Output", 
                 "Loo", 
                 "loo_M2_AW.RDS")
     )


## MODEL 3: TWO EFFECTS OF INTERVENTION

# Create the stan code via brms
# generates Stan code that can be saved in a separate file
# this part of the code is de-activated as the model specification
# for stan is already saved

#make_stancode(
#  bf(AWT.1 ~ 1) +
#    bf(AWT.2 ~ 1 + Condition) +
#    bf(AWT.3 ~ 1 + Condition) + 
#    set_rescor(rescor = TRUE),
#  data = ActiveWritingTime,
# )

# Create a proper data list to match the model

M3_AW_data <- 
  make_standata(
    bf(AWT.1 ~ 1) +
      bf(AWT.2 ~ 1 + Condition) +
      bf(AWT.3 ~ 1 + Condition) + 
      set_rescor(rescor = TRUE),
    data = ActiveWritingTime,
  )

# Activate the model based on the stan model file

Mod3_AW <- cmdstan_model(here("Stan", "M3_TwoEffects.stan"))

# fit the model with the data

fit_Mod3_AW <- Mod3_AW$sample(
  data = M3_AW_data,
  seed = 2021,
  chains = 6,
  iter_sampling = 5000,
  iter_warmup = 1500,
  parallel_chains = 6,
  refresh = 500
)

# Save the model results

Model3_AW <- fit_Mod3_AW$save_object(
  file = here("Output", 
              "Fitted_Models",
              "Model3_AW.RDS"
              )
  )

# Perform loo validation

loo_M3_AW <- fit_Mod3_AW$loo()

# Save the loo validation result

saveRDS(loo_M3_AW, 
     file = here("Output", 
                 "Loo", 
                 "loo_M3_AW.RDS"
                 )
     )

As the results of the leave-one-out cross-validation are stored in files, we can load these files and perform a model comparison.

Code

loo_M1_AW <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_AW.RDS"
             )
     )

loo_M2_AW <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_AW.RDS"
             )
     )

loo_M3_AW <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_AW.RDS"
             )
     )

loo_comparison_AW <- loo_compare(loo_M1_AW, loo_M2_AW, loo_M3_AW)

kable(
  loo_comparison_AW, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on  Active Writing Time"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Active Writing Time

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model2	0.00	0.00	186.20	15.83	11.04	2.04	-372.39	31.65
model3	-0.02	1.24	186.18	15.57	11.73	2.03	-372.36	31.14
model1	-1.14	2.08	185.06	16.04	10.13	2.02	-370.12	32.08

From the model comparison we learn that model 2, with a single effect of condition, outperforms the other models.

The following table summarizes the parameter estimates of the fixed part in the second model.

Code

Model2_AW <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model2_AW.RDS"
                       )
                     )


Draws_AW_M2 <-  as_draws_df(Model2_AW$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1_C0 = b_AWT2_Intercept - b_AWT1_Intercept,
    diff_M3_M1_C0 = b_AWT3_Intercept - b_AWT1_Intercept,
    diff_M3_M2_C0 = b_AWT3_Intercept - b_AWT2_Intercept,
    diff_M2_M1_C1 = b_AWT2_Intercept + `b_AWT[1]` - b_AWT1_Intercept,
    diff_M3_M1_C1 = b_AWT3_Intercept + `b_AWT[1]` - b_AWT1_Intercept,
    diff_M3_M2_C1 = (b_AWT3_Intercept + `b_AWT[1]`) - (b_AWT2_Intercept + `b_AWT[1]`)
  )

P <- describe_posterior(
  Draws_AW_M2,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_AWT1_Intercept" ~ "Intercept Session1",
      Parameter == "b_AWT2_Intercept" ~ "Intercept Session2",
      Parameter == "b_AWT3_Intercept" ~ "Intercept Session3",
      Parameter == "b_AWT[1]" ~ "Effect Condtion (1 = Instr.)",
      Parameter == "diff_M2_M1_C0" ~ "Contrast: Session 2 - Session 1 (F)",
      Parameter == "diff_M2_M1_C1" ~ "Contrast: Session 2 - Session 1 (FO)",
      Parameter == "diff_M3_M1_C0" ~ "Contrast: Session 3 - Session 1 (F)",
      Parameter == "diff_M3_M1_C1" ~ "Contrast: Session 3 - Session 1 (FO)",
      Parameter == "diff_M3_M2_C0" ~ "Contrast: Session 3 - Session 2 (F)",
      Parameter == "diff_M3_M2_C1" ~ "Contrast: Session 3 - Session 2 (FO)"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 2 (Active writing time) and contrasts",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 2 (Active writing time) and contrasts

Parameter	Median	CI_low	CI_high	pd
Effect Condtion (1 = Instr.)	0.043	0.009	0.075	0.982
Intercept Session1	0.639	0.611	0.666	1.000
Intercept Session2	0.707	0.680	0.733	1.000
Intercept Session3	0.725	0.697	0.752	1.000
Contrast: Session 2 - Session 1 (F)	0.069	0.036	0.098	1.000
Contrast: Session 3 - Session 1 (F)	0.086	0.053	0.117	1.000
Contrast: Session 3 - Session 2 (F)	0.017	-0.001	0.038	0.928
Contrast: Session 2 - Session 1 (FO)	0.112	0.080	0.142	1.000
Contrast: Session 3 - Session 1 (FO)	0.129	0.098	0.162	1.000
Contrast: Session 3 - Session 2 (FO)	0.017	-0.001	0.038	0.928
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

From this table we learn that the credible interval (expressed as an 89% highest density interval) of the effect of Instruction goes from 0.01 to 0.085, meaning that the 89% most probable values of the difference between both conditions in the population is situated between 1% and 8.5% more active writing time at session 2 and at session 3. The probability that the parameter estimate of the effect of condition is positive is high (pd = 0.976). In other words, given the data we can derive a high probability that students in the condition “Instruction” will score higher for active writing time in both the second and third session compared too students in the condition “Feedback only”.

Moreover, the ROPE_percentage shows that the effect of Condition has a rather high probability to be of practical significance¹ (see next table) given that the 90.0% most probable values for the effect of Condition in the population are at least equivalent to a small effect size or larger.

Code

Test <- Model2_AW$summary(
  "b_AWT[1]", 
   pr_lt_half = ~ mean(. <= 0.017))%>%
  rename(
    Parameter = variable
  ) %>%
  mutate(
    ROPE_Percentage = pr_lt_half *100
  ) %>%
  select(
    Parameter,
    ROPE_Percentage
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_AWT[1]" ~ "Effect Condtion (1 = Instruction)"
    )
  )

kable(Test, 
      digits = 3, 
      caption = "Test for practical equivalence for the effect of condition (Model 2 for Active Writing Time)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "ROPE = Region Of Practical Equivalence")

Test for practical equivalence for the effect of condition (Model 2 for Active Writing Time)

Parameter	ROPE_Percentage
Effect Condtion (1 = Instruction)	10.263
Note:
ROPE = Region Of Practical Equivalence

The following visualization shows the results of the second model. In this figure we plot the 50%, 89% and 95% credible intervals for the average scores for each session in both conditions, showing only small overlap of the credible intervals for the averages of both conditions at session 2 and 3.

Code

as_draws_df(Model2_AW$draws()) %>%
  select(b_AWT1_Intercept, 
         b_AWT2_Intercept, 
         b_AWT3_Intercept,
         `b_AWT[1]`) %>%
  mutate(Occ1_C0 = b_AWT1_Intercept,
         Occ2_C0 = b_AWT2_Intercept,
         Occ3_C0 = b_AWT3_Intercept,
         Occ1_C1 = b_AWT1_Intercept,
         Occ2_C1 = b_AWT2_Intercept + `b_AWT[1]`,
         Occ3_C1 = b_AWT3_Intercept + `b_AWT[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only",
                            "Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
  scale_color_manual(
    values = Conditions_colours[c(1,2)]) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 6: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Active Writing Time)

We can also plot the posterior probability distribution of the effect of the condition (see Figure below). This visualization supports the conclusion that we might infer that there is a rather high probability that students from both conditions differ in the amount of active writing time in session 2 and session 3.

Code

as_draws_df(Model2_AW$draws()) %>%
  select(`b_AWT[1]`) %>%
  pivot_longer(everything()) %>%
  mutate(
    Parameter = case_when(
      name == "b_AWT[1]" ~ "Effect of Instruction"
    )
  ) %>%
  ggplot(aes(x = value, 
             y = 0)) +
  stat_dotsinterval(.width = c(.50, .89, .95), 
    quantiles = 100) +
  xlab("Marginal posterior for difference between conditions in Active Writing Time") +
  ylab("") +
  facet_wrap(~Parameter)

Figure 7: Posterior probability distribution for the effect of the intervention at session 2 and 3 (Active Writing Time)

Finally, we can also compare the results with the results of a large baseline study. In the next visualization we plot the posterior probability distribution for the averages of each condition and the baseline study for each session. At session 1 students in both conditions score similar (as we modelled no difference between conditions at session 1). At each session we see that the students in the experiment demonstrated way more active writing in the first interval than the students in the baseline study.

Code

Draws_baseline_AWT <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_AWT.RDS"
  )
)

Draws_baseline_AWT <- Draws_baseline_AWT %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_AW <- as_draws_df(
  Model2_AW$draws()
)

Fig3_Posterior_AWT <- Draw_Exp_AW %>%
  cbind(Draws_baseline_AWT$b_Sess1) %>%
  cbind(Draws_baseline_AWT$b_Sess2) %>%
  cbind(Draws_baseline_AWT$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_AWT$b_Sess1`,
    Baseline2 = `Draws_baseline_AWT$b_Sess2`,
    Baseline3 = `Draws_baseline_AWT$b_Sess3`,
  ) %>%
  select(
    b_AWT1_Intercept, 
    b_AWT2_Intercept, 
    b_AWT3_Intercept,
    `b_AWT[1]`,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_AWT1_Intercept,
    Occ2_C0 = b_AWT2_Intercept,
    Occ3_C0 = b_AWT3_Intercept,
    Occ1_C1 = b_AWT1_Intercept,
    Occ2_C1 = b_AWT2_Intercept + `b_AWT[1]`,
    Occ3_C1 = b_AWT3_Intercept + `b_AWT[1]`,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Feedback only",
      grepl("C1", name) ~ "Feedback + Observation",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  ylab("") + 
  xlab("Active writing time") +
  scale_fill_manual(
    values = Conditions_colours
  ) +
  labs(fill = "Condition: ") 

Fig3_Posterior_AWT

ggsave(
  Fig3_Posterior_AWT,
  device = "jpeg",
  file = here("Figs", "Fig3_Posterior_AWT.jpg"), 
  dpi = 300
)

Figure 8: Posterior probability distributions for the average Active Writing Time at each session, for the two conditions and the baseline students.

Source time

Code

ProcesData %>%
  mutate(
    To_remove = case_when(
      is.na(Product_ID) ~ 1,
    )
  ) %>%
  filter(!is.na(Conditie) & is.na(To_remove)) %>%
  mutate(
    Condition = case_when(
      Conditie == "Feedback" ~ "Feedback_only",
      Conditie == "Instructie" ~ "Instruction"
    )
  ) %>%
  select(Product_ID, Session, Condition, Time_in_Sources) %>%
  rename(TIS = Time_in_Sources) %>%
  ggplot(aes(x = Session, y = TIS, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Time in Sources in Interval 1") 

Figure 9: Time in Sources in the first interval of the writing process for each session and both conditions

Following table summarizes the results of the model comparison.

Code

loo_M1_TIS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_TIS.RDS"
             )
     )

loo_M2_TIS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_TIS.RDS"
             )
     )

loo_M3_TIS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_TIS.RDS"
             )
     )

loo_comparison_TIS <- loo_compare(loo_M1_TIS, loo_M2_TIS, loo_M3_TIS)

kable(
  loo_comparison_TIS, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on Time In Sources"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Time In Sources

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model2	0.00	0.00	141.73	15.99	11.70	2.30	-283.47	31.98
model3	-0.64	0.73	141.09	15.78	12.43	2.29	-282.18	31.56
model1	-2.27	2.55	139.46	15.72	10.64	2.22	-278.92	31.44

From this model comparison we learn that the best fitting model is Model 2. Next we create a table summarizing the information on the parameter estimates of this model for Source time.

Code

Model2_TIS <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model2_TIS.RDS"
                       )
                     )


Draws_TIS_M2 <-  as_draws_df(Model2_TIS$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1_C0 = b_TIS2_Intercept - b_TIS1_Intercept,
    diff_M3_M1_C0 = b_TIS3_Intercept - b_TIS1_Intercept,
    diff_M3_M2_C0 = b_TIS3_Intercept - b_TIS2_Intercept,
    diff_M2_M1_C1 = b_TIS2_Intercept + `b_TIS[1]` - b_TIS1_Intercept,
    diff_M3_M1_C1 = b_TIS3_Intercept + `b_TIS[1]` - b_TIS1_Intercept,
    diff_M3_M2_C1 = (b_TIS3_Intercept + `b_TIS[1]`) - (b_TIS2_Intercept + `b_TIS[1]`)
  )

P <- describe_posterior(
  Draws_TIS_M2,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_TIS1_Intercept" ~ "Intercept Session1",
      Parameter == "b_TIS2_Intercept" ~ "Intercept Session2",
      Parameter == "b_TIS3_Intercept" ~ "Intercept Session3",
      Parameter == "b_TIS[1]" ~ "Effect Condtion (1 = Instr.)",
      Parameter == "diff_M2_M1_C0" ~ "Contrast: Session 2 - Session 1 (F)",
      Parameter == "diff_M2_M1_C1" ~ "Contrast: Session 2 - Session 1 (FO)",
      Parameter == "diff_M3_M1_C0" ~ "Contrast: Session 3 - Session 1 (F)",
      Parameter == "diff_M3_M1_C1" ~ "Contrast: Session 3 - Session 1 (FO)",
      Parameter == "diff_M3_M2_C0" ~ "Contrast: Session 3 - Session 2 (F)",
      Parameter == "diff_M3_M2_C1" ~ "Contrast: Session 3 - Session 2 (FO)"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 2 (Source time) and contrasts",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 2 (Source time) and contrasts

Parameter	Median	CI_low	CI_high	pd
Effect Condtion (1 = Instr.)	0.062	0.024	0.100	0.995
Intercept Session1	0.497	0.467	0.528	1.000
Intercept Session2	0.472	0.439	0.504	1.000
Intercept Session3	0.401	0.368	0.433	1.000
Contrast: Session 2 - Session 1 (F)	-0.025	-0.063	0.011	0.860
Contrast: Session 3 - Session 1 (F)	-0.096	-0.137	-0.056	1.000
Contrast: Session 3 - Session 2 (F)	-0.071	-0.094	-0.047	1.000
Contrast: Session 2 - Session 1 (FO)	0.037	0.004	0.071	0.962
Contrast: Session 3 - Session 1 (FO)	-0.033	-0.071	0.004	0.924
Contrast: Session 3 - Session 2 (FO)	-0.071	-0.094	-0.047	1.000
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

We can create a table with ROPE information.

Code

Test <- Model2_TIS$summary(
  c("b_TIS[1]"), 
   pr_lt_half = ~ mean(. <= 0.018))%>%
  rename(
    Parameter = variable
  ) %>%
  mutate(
    ROPE_Percentage = pr_lt_half *100
  ) %>%
  select(
    Parameter,
    ROPE_Percentage
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_TIS[1]" ~ "Effect Condtion (1 = Instruction)"
    )
  )

kable(Test, 
      digits = 3, 
      caption = "Test for practical equivalence for the effect of condition (Model 2 for Source time)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "ROPE = Region Of Practical Equivalence")

Test for practical equivalence for the effect of condition (Model 2 for Source time)

Parameter	ROPE_Percentage
Effect Condtion (1 = Instruction)	3.183
Note:
ROPE = Region Of Practical Equivalence

A visualization of the parameter estimates in the best fitting model.

Code

as_draws_df(Model2_TIS$draws()) %>%
  select(b_TIS1_Intercept, 
         b_TIS2_Intercept, 
         b_TIS3_Intercept,
         `b_TIS[1]`
         ) %>%
  mutate(Occ1_C0 = b_TIS1_Intercept,
         Occ2_C0 = b_TIS2_Intercept,
         Occ3_C0 = b_TIS3_Intercept,
         Occ1_C1 = b_TIS1_Intercept,
         Occ2_C1 = b_TIS2_Intercept + `b_TIS[1]`,
         Occ3_C1 = b_TIS3_Intercept + `b_TIS[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only","Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
  scale_color_manual(
    values = Conditions_colours[c(1,2)]
  ) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 10: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Time in Sources)

Plotting the posterior probability distribution of the effect of condition.

Code

as_draws_df(Model2_TIS$draws()) %>%
  select(`b_TIS[1]`) %>%
  pivot_longer(everything()) %>%
  mutate(
    Parameter = case_when(
      name == "b_TIS[1]" ~ "Effect of Instruction"
    )
  ) %>%
  ggplot(aes(x = value, 
             y = 0)) +
  stat_dotsinterval(.width = c(.50, .89, .95), 
    quantiles = 100) +
  xlab("Marginal posterior for difference between conditions for Time in Sources") +
  ylab("") +
  facet_wrap(~Parameter)

Figure 11: Posterior probability distribution for the effects of the intervention at session 2 and 3 (Time in Sources)

Creating the plot for the comparison with the baseline study.

Code

Draws_baseline_TIS <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_TIS.RDS"
  )
)

Draws_baseline_TIS <- Draws_baseline_TIS %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_TIS <- as_draws_df(
  Model2_TIS$draws()
)

Fig4_Posterior_TIS <- Draw_Exp_TIS %>%
  cbind(Draws_baseline_TIS$b_Sess1) %>%
  cbind(Draws_baseline_TIS$b_Sess2) %>%
  cbind(Draws_baseline_TIS$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_TIS$b_Sess1`,
    Baseline2 = `Draws_baseline_TIS$b_Sess2`,
    Baseline3 = `Draws_baseline_TIS$b_Sess3`,
  ) %>%
  select(
    b_TIS1_Intercept, 
    b_TIS2_Intercept, 
    b_TIS3_Intercept,
    `b_TIS[1]`,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_TIS1_Intercept,
    Occ2_C0 = b_TIS2_Intercept,
    Occ3_C0 = b_TIS3_Intercept,
    Occ1_C1 = b_TIS1_Intercept,
    Occ2_C1 = b_TIS2_Intercept + `b_TIS[1]`,
    Occ3_C1 = b_TIS3_Intercept + `b_TIS[1]`,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Feedback only",
      grepl("C1", name) ~ "Feedback + Observation",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  scale_fill_manual(
    values = Conditions_colours
  ) +
  ylab("") + 
  xlab("Source time") +
  labs(fill = "Condition: ") 

Fig4_Posterior_TIS

ggsave(
  Fig4_Posterior_TIS,
  device = "jpeg",
  file = here("Figs", "Fig4_Posterior_TIS.jpg"), 
  dpi = 300
  )

Figure 12: Posterior probability distributions for the average Time in Sources at each session, for the two conditions and the baseline students.

Keystrokes per minute

Code

ProcesData %>%
  mutate(
    To_remove = case_when(
      is.na(Product_ID) ~ 1,
    )
  ) %>%
  filter(!is.na(Conditie) & is.na(To_remove)) %>%
  mutate(
    Condition = case_when(
      Conditie == "Feedback" ~ "Feedback_only",
      Conditie == "Instructie" ~ "Instruction"
    )
  ) %>%
  select(Product_ID, Session, Condition, Speed) %>%
  rename(KPM = Speed) %>%
  ggplot(aes(x = Session, y = KPM, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Speed in Interval 1") 

Figure 13: Keystrokes per minute in the first interval of the writing process for each session and both conditions

Following table summarizes the results of the model comparison.

Code

loo_M1_KPM <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_KPM.RDS"
             )
     )

loo_M2_KPM <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_KPM.RDS"
             )
     )

loo_M3_KPM <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_KPM.RDS"
             )
     )

loo_comparison_KPM <- loo_compare(loo_M1_KPM, loo_M2_KPM, loo_M3_KPM)

kable(
  loo_comparison_KPM, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on Keystrokes per minute"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Keystrokes per minute

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model3	0.00	0.00	-1473.38	15.72	11.74	1.69	2946.75	31.45
model2	-0.34	1.54	-1473.71	15.92	10.95	1.69	2947.43	31.84
model1	-2.03	2.78	-1475.41	15.50	9.82	1.45	2950.82	30.99

From this model comparison we learn that the best fitting model is Model 3. Next we create a table with a summary of the parameter estimates information.

Code

Model3_KPM <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model3_KPM.RDS"
                       )
                     )


Draws_KPM_M3 <-  as_draws_df(Model3_KPM$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1_C0 = b_KPM2_Intercept - b_KPM1_Intercept,
    diff_M3_M1_C0 = b_KPM3_Intercept - b_KPM1_Intercept,
    diff_M3_M2_C0 = b_KPM3_Intercept - b_KPM2_Intercept,
    diff_M2_M1_C1 = b_KPM2_Intercept + `b_KPM2[1]` - b_KPM1_Intercept,
    diff_M3_M1_C1 = b_KPM3_Intercept + `b_KPM3[1]` - b_KPM1_Intercept,
    diff_M3_M2_C1 = (b_KPM3_Intercept + `b_KPM3[1]`) - (b_KPM2_Intercept + `b_KPM2[1]`)
  )

P <- describe_posterior(
  Draws_KPM_M3,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_KPM1_Intercept" ~ "Intercept Session1",
      Parameter == "b_KPM2_Intercept" ~ "Intercept Session2",
      Parameter == "b_KPM3_Intercept" ~ "Intercept Session3",
      Parameter == "b_KPM2[1]" ~ "Effect Condtion (1 = Instr.) Session 2",
      Parameter == "b_KPM3[1]" ~ "Effect Condtion (1 = Instr.) Session 3",
      Parameter == "diff_M2_M1_C0" ~ "Contrast: Session 2 - Session 1 (F)",
      Parameter == "diff_M2_M1_C1" ~ "Contrast: Session 2 - Session 1 (FO)",
      Parameter == "diff_M3_M1_C0" ~ "Contrast: Session 3 - Session 1 (F)",
      Parameter == "diff_M3_M1_C1" ~ "Contrast: Session 3 - Session 1 (FO)",
      Parameter == "diff_M3_M2_C0" ~ "Contrast: Session 3 - Session 2 (F)",
      Parameter == "diff_M3_M2_C1" ~ "Contrast: Session 3 - Session 2 (FO)"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 2 (Speed) and contrasts",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 2 (Speed) and contrasts

Parameter	Median	CI_low	CI_high	pd
Effect Condtion (1 = Instr.) Session 2	-7.468	-13.895	-1.541	0.973
Effect Condtion (1 = Instr.) Session 3	-13.112	-20.253	-5.773	0.998
Intercept Session1	45.105	40.831	49.213	1.000
Intercept Session2	53.390	48.826	58.279	1.000
Intercept Session3	70.444	65.073	75.998	1.000
Contrast: Session 2 - Session 1 (F)	8.289	3.465	13.321	0.996
Contrast: Session 3 - Session 1 (F)	25.340	19.873	30.984	1.000
Contrast: Session 3 - Session 2 (F)	17.072	12.963	21.226	1.000
Contrast: Session 2 - Session 1 (FO)	0.803	-3.944	5.743	0.608
Contrast: Session 3 - Session 1 (FO)	12.247	6.449	17.475	1.000
Contrast: Session 3 - Session 2 (FO)	11.404	7.383	15.585	1.000
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

Creating a table with the ROPE information.

Code

Test <- Model3_KPM$summary(
  c("b_KPM2[1]","b_KPM3[1]") , 
   pr_lt_half = ~ mean(. >= -2.74))%>%
  rename(
    Parameter = variable
  ) %>%
  mutate(
    ROPE_Percentage = pr_lt_half *100
  ) %>%
  select(
    Parameter,
    ROPE_Percentage
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_KPM2[1]" ~ "Effect Condtion Session 2 (1 = Instruction)",
      Parameter == "b_KPM3[1]" ~ "Effect Condtion Session 3 (1 = Instruction)",
    )
  )

kable(Test, 
      digits = 3, 
      caption = "Test for practical equivalence for the effect of condition (Model 3 for Speed)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "ROPE = Region Of Practical Equivalence")

Test for practical equivalence for the effect of condition (Model 3 for Speed)

Parameter	ROPE_Percentage
Effect Condtion Session 2 (1 = Instruction)	10.97
Effect Condtion Session 3 (1 = Instruction)	1.17
Note:
ROPE = Region Of Practical Equivalence

Plotting the parameter estimates of the best fitting model.

Code

as_draws_df(Model3_KPM$draws()) %>%
  select(b_KPM1_Intercept, 
         b_KPM2_Intercept, 
         b_KPM3_Intercept,
         `b_KPM2[1]`,
         `b_KPM3[1]`
         ) %>%
  mutate(Occ1_C0 = b_KPM1_Intercept,
         Occ2_C0 = b_KPM2_Intercept,
         Occ3_C0 = b_KPM3_Intercept,
         Occ1_C1 = b_KPM1_Intercept,
         Occ2_C1 = b_KPM2_Intercept + `b_KPM2[1]`,
         Occ3_C1 = b_KPM3_Intercept + `b_KPM3[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only","Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
  scale_colour_manual(
    values = Conditions_colours[c(1,2)]
  ) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 14: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Keystrokes per minute)

Plotting the effects of the condition.

Code

as_draws_df(Model3_KPM$draws()) %>%
  select(c(`b_KPM2[1]`,`b_KPM3[1]`)) %>%
  pivot_longer(everything()) %>%
  mutate(
    Parameter = case_when(
      name == "b_KPM2[1]" ~ "Effect of Instruction Session 2",
      name == "b_KPM3[1]" ~ "Effect of Instruction Session 3"
    )
  ) %>%
  ggplot(aes(x = value, 
             y = 0)) +
  stat_dotsinterval(.width = c(.50, .89, .95), 
    quantiles = 100) +
  xlab("Marginal posterior for difference between conditions in Speed") +
  ylab("") +
  facet_wrap(~Parameter)

Figure 15: Posterior probability distribution for the effects of the intervention at both sessions (Keystrokes per minute)

Making a visual comparison with the baseline study.

Code

Draws_baseline_KPM <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_SPM.RDS"
  )
)

Draws_baseline_KPM <- Draws_baseline_KPM %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_KPM <- as_draws_df(
  Model3_KPM$draws()
)

Fig5_Posterior_KPM <- Draw_Exp_KPM %>%
  cbind(Draws_baseline_KPM$b_Sess1) %>%
  cbind(Draws_baseline_KPM$b_Sess2) %>%
  cbind(Draws_baseline_KPM$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_KPM$b_Sess1`,
    Baseline2 = `Draws_baseline_KPM$b_Sess2`,
    Baseline3 = `Draws_baseline_KPM$b_Sess3`,
  ) %>%
  select(
    b_KPM1_Intercept, 
    b_KPM2_Intercept, 
    b_KPM3_Intercept,
    `b_KPM2[1]`,
    `b_KPM3[1]`,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_KPM1_Intercept,
    Occ2_C0 = b_KPM2_Intercept,
    Occ3_C0 = b_KPM3_Intercept,
    Occ1_C1 = b_KPM1_Intercept,
    Occ2_C1 = b_KPM2_Intercept + `b_KPM2[1]`,
    Occ3_C1 = b_KPM3_Intercept + `b_KPM3[1]`,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Feedback only",
      grepl("C1", name) ~ "Feedback + Observation",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  scale_fill_manual(
    values = Conditions_colours
  ) +
  ylab("") + 
  xlab("Speed") +
  labs(fill = "Condition: ") 

Fig5_Posterior_KPM

ggsave(
  Fig5_Posterior_KPM,
  device = "jpeg",
  file = here("Figs", "Fig5_Posterior_KPM.jpg"), 
  dpi = 300
  )

Figure 16: Posterior probability distributions for the average Keystrokes per minute at each session, for the two conditions and the baseline students.

Switches between sources

Code

ProcesData %>%
  mutate(
    To_remove = case_when(
      is.na(Product_ID) ~ 1,
    )
  ) %>%
  filter(!is.na(Conditie) & is.na(To_remove)) %>%
  mutate(
    Condition = case_when(
      Conditie == "Feedback" ~ "Feedback_only",
      Conditie == "Instructie" ~ "Instruction"
    )
  ) %>%
  select(Product_ID, Session, Condition, Sources_Switch) %>%
  rename(SBS = Sources_Switch) %>%
  ggplot(aes(x = Session, y = SBS, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Sources sources in Interval 1") 

Figure 17: Time in Sources in the first interval of the writing process for each session and both conditions

Following table summarizes the results of the model comparison.

Code

loo_M1_SBS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_SBS.RDS"
             )
     )

loo_M2_SBS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_SBS.RDS"
             )
     )

loo_M3_SBS <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_SBS.RDS"
             )
     )

loo_comparison_SBS <- loo_compare(loo_M1_SBS, loo_M2_SBS, loo_M3_SBS)

kable(
  loo_comparison_SBS, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on Switches sources"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Switches sources

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model3	0.00	0.00	-187.62	25.61	17.10	6.28	375.24	51.21
model2	-0.82	1.68	-188.44	25.72	16.43	6.29	376.88	51.43
model1	-1.14	2.36	-188.76	25.74	15.44	6.11	377.52	51.48

From this model comparison we learn that the best fitting model is Model 3. Next we create a table that summarizes the parameter estimates for this model.

Code

Model3_SBS <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model3_SBS.RDS"
                       )
                     )


Draws_SBS_M3 <-  as_draws_df(Model3_SBS$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1_C0 = b_SBS2_Intercept - b_SBS1_Intercept,
    diff_M3_M1_C0 = b_SBS3_Intercept - b_SBS1_Intercept,
    diff_M3_M2_C0 = b_SBS3_Intercept - b_SBS2_Intercept,
    diff_M2_M1_C1 = b_SBS2_Intercept + `b_SBS2[1]` - b_SBS1_Intercept,
    diff_M3_M1_C1 = b_SBS3_Intercept + `b_SBS3[1]` - b_SBS1_Intercept,
    diff_M3_M2_C1 = (b_SBS3_Intercept + `b_SBS3[1]`) - (b_SBS2_Intercept + `b_SBS2[1]`)
  )

P <- describe_posterior(
  Draws_SBS_M3,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_SBS1_Intercept" ~ "Intercept Session1",
      Parameter == "b_SBS2_Intercept" ~ "Intercept Session2",
      Parameter == "b_SBS3_Intercept" ~ "Intercept Session3",
      Parameter == "b_SBS2[1]" ~ "Effect Condtion (1 = Instr.) Session 2",
      Parameter == "b_SBS3[1]" ~ "Effect Condtion (1 = Instr.) Session 3",
      Parameter == "diff_M2_M1_C0" ~ "Contrast: Session 2 - Session 1 (F)",
      Parameter == "diff_M2_M1_C1" ~ "Contrast: Session 2 - Session 1 (FO)",
      Parameter == "diff_M3_M1_C0" ~ "Contrast: Session 3 - Session 1 (F)",
      Parameter == "diff_M3_M1_C1" ~ "Contrast: Session 3 - Session 1 (FO)",
      Parameter == "diff_M3_M2_C0" ~ "Contrast: Session 3 - Session 2 (F)",
      Parameter == "diff_M3_M2_C1" ~ "Contrast: Session 3 - Session 2 (FO)"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 2 (Switches sources) and contrasts",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 2 (Switches sources) and contrasts

Parameter	Median	CI_low	CI_high	pd
Effect Condtion (1 = Instr.) Session 2	0.050	-0.048	0.153	0.788
Effect Condtion (1 = Instr.) Session 3	0.204	0.074	0.342	0.993
Intercept Session1	0.572	0.483	0.661	1.000
Intercept Session2	0.430	0.355	0.505	1.000
Intercept Session3	0.493	0.393	0.591	1.000
Contrast: Session 2 - Session 1 (F)	-0.142	-0.253	-0.033	0.980
Contrast: Session 3 - Session 1 (F)	-0.079	-0.202	0.052	0.840
Contrast: Session 3 - Session 2 (F)	0.063	-0.042	0.167	0.836
Contrast: Session 2 - Session 1 (FO)	-0.092	-0.196	0.014	0.920
Contrast: Session 3 - Session 1 (FO)	0.125	0.008	0.248	0.951
Contrast: Session 3 - Session 2 (FO)	0.217	0.121	0.310	1.000
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

We can create a table with information on the ROPE.

Code

Test <- Model3_SBS$summary(
  c("b_SBS2[1]", "b_SBS3[1]"), 
   pr_lt_half = ~ mean(. <= 0.045))%>%
  rename(
    Parameter = variable
  ) %>%
  mutate(
    ROPE_Percentage = pr_lt_half *100
  ) %>%
  select(
    Parameter,
    ROPE_Percentage
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_SBS2[1]" ~ "Effect Condtion Session 2(1 = Instruction)",
      Parameter == "b_SBS3[1]" ~ "Effect Condtion Session 3(1 = Instruction)"
    )
  )

kable(Test, 
      digits = 3, 
      caption = "Test for practical equivalence for the effect of condition (Model 3 for Switches sources)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "ROPE = Region Of Practical Equivalence")

Test for practical equivalence for the effect of condition (Model 3 for Switches sources)

Parameter	ROPE_Percentage
Effect Condtion Session 2(1 = Instruction)	46.78
Effect Condtion Session 3(1 = Instruction)	2.96
Note:
ROPE = Region Of Practical Equivalence

Creating a plot that visualizes the parameter estimates for the model.

Code

as_draws_df(Model3_SBS$draws()) %>%
  select(b_SBS1_Intercept, 
         b_SBS2_Intercept, 
         b_SBS3_Intercept,
         `b_SBS2[1]`,
         `b_SBS3[1]`
         ) %>%
  mutate(Occ1_C0 = b_SBS1_Intercept,
         Occ2_C0 = b_SBS2_Intercept,
         Occ3_C0 = b_SBS3_Intercept,
         Occ1_C1 = b_SBS1_Intercept,
         Occ2_C1 = b_SBS2_Intercept + `b_SBS2[1]`,
         Occ3_C1 = b_SBS3_Intercept + `b_SBS3[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only","Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
   scale_colour_manual(
    values = Conditions_colours[c(1,2)]
  ) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 18: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Switches between sources)

Creating a plot that visualizes the comparison with the baseline study.

Code

Draws_baseline_SBS <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_SBS.RDS"
  )
)

Draws_baseline_SBS <- Draws_baseline_SBS %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_SBS <- as_draws_df(
  Model3_SBS$draws()
)

Fig6_Posterior_SBS <- Draw_Exp_SBS %>%
  cbind(Draws_baseline_SBS$b_Sess1) %>%
  cbind(Draws_baseline_SBS$b_Sess2) %>%
  cbind(Draws_baseline_SBS$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_SBS$b_Sess1`,
    Baseline2 = `Draws_baseline_SBS$b_Sess2`,
    Baseline3 = `Draws_baseline_SBS$b_Sess3`,
  ) %>%
  select(
    b_SBS1_Intercept, 
    b_SBS2_Intercept, 
    b_SBS3_Intercept,
    `b_SBS2[1]`,
    `b_SBS3[1]`,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_SBS1_Intercept,
    Occ2_C0 = b_SBS2_Intercept,
    Occ3_C0 = b_SBS3_Intercept,
    Occ1_C1 = b_SBS1_Intercept,
    Occ2_C1 = b_SBS2_Intercept + `b_SBS2[1]`,
    Occ3_C1 = b_SBS3_Intercept + `b_SBS3[1]`,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Feedback only",
      grepl("C1", name) ~ "Feedback + Observation",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  scale_fill_manual(
    values = Conditions_colours
  ) +
  ylab("") + 
  xlab("Switches between sources") +
  labs(fill = "Condition: ") 

Fig6_Posterior_SBS

ggsave(
  Fig6_Posterior_SBS,
  device = "jpeg",
  file = here("Figs", "Fig6_Posterior_SBS.jpg"), 
  dpi = 300
  )

Figure 19: Posterior probability distributions for the average number of switches between sources at each session, for the two conditions and the baseline students.

Switches between sources and text produced so far

Code

ProcesData %>%
  mutate(
    To_remove = case_when(
      is.na(Product_ID) ~ 1,
    )
  ) %>%
  filter(!is.na(Conditie) & is.na(To_remove)) %>%
  mutate(
    Condition = case_when(
      Conditie == "Feedback" ~ "Feedback_only",
      Conditie == "Instructie" ~ "Instruction"
    )
  ) %>%
  select(Product_ID, Session, Condition, Source_Text_Switch) %>%
  rename(KPM = Source_Text_Switch) %>%
  ggplot(aes(x = Session, y = KPM, group = Product_ID)) +
  geom_path() +
  facet_wrap(.~Condition) +
  scale_x_continuous(breaks = c(1,2,3)) +
  labs(x = "Session", y = "Switches sources - synthesis in Interval 1") 

Figure 20: Switches between sources and text in the first interval of the writing process for each session and both conditions

Following table summarizes the results of the model comparison.

Code

loo_M1_SST <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M1_SST.RDS"
             )
     )

loo_M2_SST <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M2_SST.RDS"
             )
     )

loo_M3_SST <- readRDS(
  file = here("Output", 
              "Loo", 
              "loo_M3_SST.RDS"
             )
     )

loo_comparison_SST <- loo_compare(loo_M1_SST, loo_M2_SST, loo_M3_SST)

kable(
  loo_comparison_SST, 
  digits = 2,
  caption = "Model comparisons based on Leave-one-out cross validation for the models on Switches sources-synthesis"
  ) %>%
  kable_classic(html_font = "Calibri", full_width = T)

Model comparisons based on Leave-one-out cross validation for the models on Switches sources-synthesis

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model2	0.00	0.00	-368.23	18.38	12.89	3.03	736.45	36.76
model3	-1.14	0.09	-369.37	18.41	13.89	3.21	738.74	36.82
model1	-1.88	1.99	-370.11	18.60	12.31	3.01	740.22	37.20

From this model comparison we learn that the best fitting model is Model 2. We will report on Model 2.

Code

Model2_SST <- readRDS(file = 
                       here(
                         "Output",
                         "Fitted_Models",
                         "Model2_SST.RDS"
                       )
                     )


Draws_SST_M2 <-  as_draws_df(Model2_SST$draws()) %>%
  select(
    starts_with("b_")
  ) %>%
  mutate(
    diff_M2_M1_C0 = b_SST2_Intercept - b_SST1_Intercept,
    diff_M3_M1_C0 = b_SST3_Intercept - b_SST1_Intercept,
    diff_M3_M2_C0 = b_SST3_Intercept - b_SST2_Intercept,
    diff_M2_M1_C1 = b_SST2_Intercept + `b_SST[1]` - b_SST1_Intercept,
    diff_M3_M1_C1 = b_SST3_Intercept + `b_SST[1]` - b_SST1_Intercept,
    diff_M3_M2_C1 = (b_SST3_Intercept + `b_SST[1]`) - (b_SST2_Intercept + `b_SST[1]`)
  )

P <- describe_posterior(
  Draws_SST_M2,
  ci = 0.89,
  ci_method = "hdi") %>%
  select(
    Parameter, 
    Median, 
    CI_low, 
    CI_high, 
    pd
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_SST1_Intercept" ~ "Intercept Session1",
      Parameter == "b_SST2_Intercept" ~ "Intercept Session2",
      Parameter == "b_SST3_Intercept" ~ "Intercept Session3",
      Parameter == "b_SST[1]" ~ "Effect Condtion (1 = Instr.)",
      Parameter == "diff_M2_M1_C0" ~ "Contrast: Session 2 - Session 1 (F)",
      Parameter == "diff_M2_M1_C1" ~ "Contrast: Session 2 - Session 1 (FO)",
      Parameter == "diff_M3_M1_C0" ~ "Contrast: Session 3 - Session 1 (F)",
      Parameter == "diff_M3_M1_C1" ~ "Contrast: Session 3 - Session 1 (FO)",
      Parameter == "diff_M3_M2_C0" ~ "Contrast: Session 3 - Session 2 (F)",
      Parameter == "diff_M3_M2_C1" ~ "Contrast: Session 3 - Session 2 (FO)"
    )
  )

kable(P, 
      digits = 3, 
      caption = "Parameter estimates information for Model 2 (Switches sources-synthesis) and contrasts",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction")

Parameter estimates information for Model 2 (Switches sources-synthesis) and contrasts

Parameter	Median	CI_low	CI_high	pd
Effect Condtion (1 = Instr.)	0.382	0.108	0.645	0.988
Intercept Session1	1.152	0.990	1.313	1.000
Intercept Session2	1.486	1.252	1.700	1.000
Intercept Session3	1.657	1.407	1.913	1.000
Contrast: Session 2 - Session 1 (F)	0.332	0.103	0.552	0.990
Contrast: Session 3 - Session 1 (F)	0.503	0.244	0.755	0.999
Contrast: Session 3 - Session 2 (F)	0.171	0.002	0.342	0.948
Contrast: Session 2 - Session 1 (FO)	0.715	0.515	0.910	1.000
Contrast: Session 3 - Session 1 (FO)	0.885	0.649	1.112	1.000
Contrast: Session 3 - Session 2 (FO)	0.171	0.002	0.342	0.948
Note:
CI_low = lower limit of a 89% highest density interval; CI_high = upper limit of a 89% highest density interval; pd = probability of direction

Creating the ROPE information table.

Code

Test <- Model2_SST$summary(
  "b_SST[1]", 
   pr_lt_half = ~ mean(. <= 0.106))%>%
  rename(
    Parameter = variable
  ) %>%
  mutate(
    ROPE_Percentage = pr_lt_half *100
  ) %>%
  select(
    Parameter,
    ROPE_Percentage
  ) %>%
  mutate(
    Parameter = case_when(
      Parameter == "b_SST[1]" ~ "Effect Condtion (1 = Instruction)",
    )
  )

kable(Test, 
      digits = 3, 
      caption = "Test for practical equivalence for the effect of condition (Model 2 for Switches sources-synthesis)",
      row.names = F) %>% 
  kable_classic(
    html_font = "Calibri", 
    full_width = T) %>%
  footnote(
    general = "ROPE = Region Of Practical Equivalence")

Test for practical equivalence for the effect of condition (Model 2 for Switches sources-synthesis)

Parameter	ROPE_Percentage
Effect Condtion (1 = Instruction)	5.027
Note:
ROPE = Region Of Practical Equivalence

Plotting the parameter estimates of the best fitting model.

Code

as_draws_df(Model2_SST$draws()) %>%
  select(b_SST1_Intercept, 
         b_SST2_Intercept, 
         b_SST3_Intercept,
         `b_SST[1]`
         ) %>%
  mutate(Occ1_C0 = b_SST1_Intercept,
         Occ2_C0 = b_SST2_Intercept,
         Occ3_C0 = b_SST3_Intercept,
         Occ1_C1 = b_SST1_Intercept,
         Occ2_C1 = b_SST2_Intercept + `b_SST[1]`,
         Occ3_C1 = b_SST3_Intercept + `b_SST[1]`) %>%
  select(Occ1_C0, Occ2_C0, Occ3_C0, Occ1_C1, Occ2_C1, Occ3_C1) %>%
  pivot_longer(everything()) %>%
  mutate(Occ1 = ifelse(grepl("Occ1",name), 1, 0),
         Occ2 = ifelse(grepl("Occ2",name), 1, 0),
         Occ3 = ifelse(grepl("Occ3",name), 1, 0),
         Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
         Condition = ifelse(grepl("C0", name), "Feedback only","Feedback + Observation")
         )  %>%
  ggplot(aes(x = value, 
             y = Occ,
             color = as.factor(Condition))) +
  stat_pointinterval(position = "dodgejust",
                     .width = c(.50, .89, .95)) +
  scale_fill_ramp_discrete(na.translate = FALSE) +
   scale_colour_manual(
    values = Conditions_colours[c(1,2)]
  ) +
  ylab("Occasion") + 
  xlab("Marginal posterior") +
  coord_flip() +
  labs(color = "Condition: ") +
  theme(legend.position = "top")

Figure 21: Credible intervals (50%, 89% and 95%) for the averages at each session for both conditions (Switches between sources and text)

Plotting the information of the baseline study together with the posterior draws from the best fitting model.

Code

Draws_baseline_SST <- readRDS(
  file = here(
    "Output",
    "Baseline_models",
    "Draws_baseline_SST.RDS"
  )
)

Draws_baseline_SST <- Draws_baseline_SST %>%
  select(
    b_Intercept,
    `r_Task_Code[A1,Intercept]`,
    `r_Task_Code[A2,Intercept]`,
    `r_Task_Code[A3,Intercept]`,
    `r_Task_Code[A4,Intercept]`,
    `r_Task_Code[B1,Intercept]`,
    `r_Task_Code[B2,Intercept]`,
    `r_Task_Code[B3,Intercept]`,
    `r_Task_Code[B4,Intercept]`,
    `r_Task_Code[D1,Intercept]`,
    `r_Task_Code[D2,Intercept]`,
    `r_Task_Code[D3,Intercept]`,
    `r_Task_Code[D4,Intercept]`
  ) %>%
  mutate(
    b_Sess1 = ((b_Intercept + `r_Task_Code[B1,Intercept]`) +
              (b_Intercept + `r_Task_Code[B2,Intercept]`) +
              (b_Intercept + `r_Task_Code[B3,Intercept]`) +
              (b_Intercept + `r_Task_Code[B4,Intercept]`))/4,
    b_Sess2 = ((b_Intercept + `r_Task_Code[A1,Intercept]`) +
              (b_Intercept + `r_Task_Code[A2,Intercept]`) +
              (b_Intercept + `r_Task_Code[A3,Intercept]`) +
              (b_Intercept + `r_Task_Code[A4,Intercept]`))/4,
    b_Sess3 = ((b_Intercept + `r_Task_Code[D1,Intercept]`) +
              (b_Intercept + `r_Task_Code[D2,Intercept]`) +
              (b_Intercept + `r_Task_Code[D3,Intercept]`) +
              (b_Intercept + `r_Task_Code[D4,Intercept]`))/4,
  ) %>%
  select(
    b_Intercept,
    b_Sess1,
    b_Sess2,
    b_Sess3
  ) 


Draw_Exp_SST <- as_draws_df(
  Model2_SST$draws()
)

Fig7_Posterior_SST <- Draw_Exp_SST %>%
  cbind(Draws_baseline_SST$b_Sess1) %>%
  cbind(Draws_baseline_SST$b_Sess2) %>%
  cbind(Draws_baseline_SST$b_Sess3) %>%
  rename(
    Baseline1 = `Draws_baseline_SST$b_Sess1`,
    Baseline2 = `Draws_baseline_SST$b_Sess2`,
    Baseline3 = `Draws_baseline_SST$b_Sess3`,
  ) %>%
  select(
    b_SST1_Intercept, 
    b_SST2_Intercept, 
    b_SST3_Intercept,
    `b_SST[1]`,
    Baseline1,
    Baseline2,
    Baseline3
    ) %>%
  mutate(
    Occ1_C0 = b_SST1_Intercept,
    Occ2_C0 = b_SST2_Intercept,
    Occ3_C0 = b_SST3_Intercept,
    Occ1_C1 = b_SST1_Intercept,
    Occ2_C1 = b_SST2_Intercept + `b_SST[1]`,
    Occ3_C1 = b_SST3_Intercept + `b_SST[1]`,
    Occ1_Base = Baseline1,
    Occ2_Base = Baseline2,
    Occ3_Base = Baseline3
    ) %>%
  select(
    Occ1_C0, 
    Occ2_C0, 
    Occ3_C0, 
    Occ1_C1, 
    Occ2_C1, 
    Occ3_C1,
    Occ1_Base,
    Occ2_Base,
    Occ3_Base) %>%
  pivot_longer(everything()) %>%
  mutate(
    Occ1 = ifelse(grepl("Occ1",name), 1, 0),
    Occ2 = ifelse(grepl("Occ2",name), 1, 0),
    Occ3 = ifelse(grepl("Occ3",name), 1, 0),
    Occ  = factor(Occ1 + 2*Occ2 + 3*Occ3),
    Condition = case_when(
      grepl("C0", name) ~ "Feedback only",
      grepl("C1", name) ~ "Feedback + Observation",
      grepl("Base", name) ~ "Baseline")
      ) %>%
  mutate(
    Occ = case_when(
      Occ == "1" ~ "Session 1",
      Occ == "2" ~ "Session 2",
      Occ == "3" ~ "Session 3"
    )
  ) %>%
  ggplot(aes(
    x = value, 
    fill = as.factor(Condition))) +
  stat_halfeye(
    .width = c(.50, .89, .95),
    alpha = 0.8
    ) +
  facet_wrap(.~ Occ) +
  scale_fill_manual(
    values = Conditions_colours
  ) +
  ylab("") + 
  xlab("Switches sources-synthesis") +
  labs(fill = "Condition: ") 

Fig7_Posterior_SST

ggsave(
  Fig7_Posterior_SST,
    device = "jpeg",
    file = here("Figs", "Fig7_Posterior_SST.jpg"), 
    dpi = 300
    )

Figure 22: Posterior probability distributions for the average number of switches between sources and text produced so far, at each session, for the two conditions and the baseline students.

Integrating figures

To create the integrated figure of the comparisons with the baseline study we used patchwork. Following code creates this figure and saves it in a publication-ready format.

Code

Fig3_Posterior_AWT + Fig4_Posterior_TIS + Fig5_Posterior_KPM + Fig6_Posterior_SBS + Fig7_Posterior_SST + guide_area( ) + plot_layout(ncol = 2, guides = "collect")

Code

ggsave(
  width = 14,
  height = 10,
    device = "jpeg",
    file = here("Figs", "Fig8_Process_variables.jpg"), 
    dpi = 300
    )

References

Bürkner, Paul-Christian. 2017. “Brms: An r Package for Bayesian Multilevel Models Using Stan” 80. https://doi.org/10.18637/jss.v080.i01.

Gelman, Andrew, Jennifer Hill, and Aki Vehtari. 2021. Regression and other stories. Analytical methods for social research. Cambridge New York, NY Port Melbourne, VIC New Delhi Singapore: Cambridge University Press. https://doi.org/10.1017/9781139161879.

Kay, Matthew. 2022. “Ggdist: Visualizations of Distributions and Uncertainty.” https://doi.org/10.5281/zenodo.3879620.

Makowski, Dominique, Mattan S. Ben-Shachar, and Daniel Lüdecke. 2019. “bayestestR: Describing Effects and Their Uncertainty, Existence and Significance Within the Bayesian Framework.” Journal of Open Source Software 4 (40): 1541. https://doi.org/10.21105/joss.01541.

Müller, Kirill. 2020. “Here: A Simpler Way to Find Your Files.” https://CRAN.R-project.org/package=here.

Pedersen, Thomas Lin. 2022. “Patchwork: The Composer of Plots.” https://CRAN.R-project.org/package=patchwork.

R Core Team. 2022. “Foreign: Read Data Stored by ’Minitab’, ’s’, ’SAS’, ’SPSS’, ’Stata’, ’Systat’, ’Weka’, ’dBase’, ...” https://CRAN.R-project.org/package=foreign.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse” 4: 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, and Jennifer Bryan. 2022. “Readxl: Read Excel Files.” https://CRAN.R-project.org/package=readxl.

Zhu, Hao. 2021. “kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax.” https://CRAN.R-project.org/package=kableExtra.

Footnotes

1. To define the ROPE (region of practical equivalence) we used the suggestion to call parameter estimates that correspond to a small effect size or larger (e.g., Cohen’s d > |0.1|) to be of ‘practical equivalence’.