setup.Rmd

---
title: "Setup & Data Preparation"
author: Hampus Broman & William Levén
date: 2021-05
output: 
  html_document: 
    pandoc_args: [ "-o", "docs/setup.html" ]
---
## Libraries

```{r load-libraries, class.source = 'fold-show', message=FALSE}
# Data manipulation
library(plyr)
library(dplyr)
library(tidyr)
library(reshape2)
library(tibble)

# Pretty printing
library(scales)

# Fetching data
library(curl)

# Computing hashes, used for efficient model caching
library(hashr)

# Easier plotting
library(ggplot2); theme_set(theme_minimal())
library(likert)
library(bayesplot)
library(ggpubr)
library(HDInterval)
library(ggridges)

# Baysian modeling
library(brms)
```

## Data fetching

We begin by importing the csv data from the [data repository](https://github.com/BrokenWindowsInvestigation/Data).

```{r fetch-data}
d.orig <- read.csv(curl("https://raw.githubusercontent.com/BrokenWindowsInvestigation/Data/master/data.csv"))
d.orig
```

## Encoding the data with correct types

```{r data-encoding}
### Utility functions for encoding ###
encode.categorical <- function(column, categories) {
  factor(column, level = categories)
}

encode.bool <- function(column) {
  encode.categorical(column, c("true", "false"))
}
encode.logic <- function(column) {
  encode.categorical(column, c(TRUE, FALSE))
}

encode.orderedcategorical <- function(column, categories) {
  as.ordered(encode.categorical(column, categories))
}

encode.likert <- function(column) {
  encode.orderedcategorical(column, c(-3, -2, -1, 0, 1, 2, 3))
}

### Encode the original data ###
d <- data.frame(
  session = factor(d.orig$session), 
  
  time = d.orig$time,
  
  reused_logic_constructor = encode.bool(d.orig$reused_logic_constructor),
  reused_logic_validation  = encode.bool(d.orig$reused_logic_validation),
  
  equals.state   = encode.orderedcategorical(
    d.orig$equals_state, 
    c("Not implemented", "Duplicated", "Good")
  ),
  hashcode.state = encode.orderedcategorical(
    d.orig$hashcode_state, 
    c("Not implemented", "Duplicated", "Good")
  ),
  
  documentation = factor(d.orig$documentation),
  
  var_names_copied_all        = d.orig$var_names_copied_all,
  var_names_copied_good       = d.orig$var_names_copied_good,
  var_names_copied_good.ratio = d.orig$var_names_copied_good / d.orig$var_names_copied_all,
  var_names_new_all           = d.orig$var_names_new_all,
  var_names_new_good          = d.orig$var_names_new_good,
  var_names_new_good.ratio    = d.orig$var_names_new_good / d.orig$var_names_new_all,
  var_names_edited_all        = d.orig$var_names_edited_all,
  var_names_edited_good       = d.orig$var_names_edited_good,
  var_names_edited_good.ratio = d.orig$var_names_edited_good / d.orig$var_names_edited_all,
  
  sonarqube_issues          = 
    d.orig$sonarqube_issues_major + 
    d.orig$sonarqube_issues_minor + 
    d.orig$sonarqube_issues_info + 
    d.orig$sonarqube_issues_critical,
  sonarqube_issues.major    = d.orig$sonarqube_issues_major,
  sonarqube_issues.minor    = d.orig$sonarqube_issues_minor,
  sonarqube_issues.info     = d.orig$sonarqube_issues_info,
  sonarqube_issues.critical = d.orig$sonarqube_issues_critical,
  
  group = factor(d.orig$group),
  
  education_level = encode.orderedcategorical(d.orig$education_level, c(
    "None", 
    "Some bachelor studies", 
    "Bachelor degree", 
    "Some master studies", 
    "Master degree", 
    "Some Ph.D. studies", 
    "Ph. D."
  )),
  education_field = factor(d.orig$education_field),
  
  work_domain                 = factor(d.orig$work_domain),
  work_experience_programming = d.orig$work_experience_programming,
  work_experience_java        = d.orig$work_experience_java,
  
  workplace_pair_programming = encode.bool(d.orig$workplace_pair_programming),
  workplace_peer_review      = encode.bool(d.orig$workplace_peer_review),
  workplace_td_tracking      = encode.bool(d.orig$workplace_td_tracking),
  workplace_coding_standards = encode.bool(d.orig$workplace_coding_standards),
  
  task_completion = encode.orderedcategorical(d.orig$task_completion, c(
    "Not submitted", 
    "Does not compile", 
    "Invalid solution", 
    "Completed"
  )),
  
  quality_pre_task  = encode.likert(d.orig$quality_pre_task),
  quality_post_task = encode.likert(d.orig$quality_post_task),
  
  high_debt_version = encode.bool(d.orig$high_debt_version),
  scenario          = encode.categorical(d.orig$scenario, c("booking", "tickets")),
  order             = encode.orderedcategorical(d.orig$order, c(0, 1)),
  
  modified_lines         = d.orig$modified_lines,
  large_structure_change = encode.bool(d.orig$large_structure_change)
)

d$equals.exists <- encode.logic(d$equals.state != "Not implemented")
d$hashcode.exists <-encode.logic(d$hashcode.state != "Not implemented")

str(d)
```

## Partial data sets and aggregates {.tabset}

For some models partial data sets and aggregates are needed.

### Sessions as rows

```{r partial-sessions}
d.sessions <- d %>% group_by(session) %>% dplyr::summarise(
  across(task_completion, min),
  across(c(
    education_level, 
    education_field, 
    work_domain, 
    group,
    work_experience_java, 
    work_experience_programming, 
    workplace_coding_standards, 
    workplace_pair_programming, 
    workplace_peer_review, 
    workplace_td_tracking
  ), first)
  ) 

d$work_experience_programming.s = scale(d$work_experience_programming)
d$work_experience_java.s = scale(d$work_experience_java)

d.sessions
```

### Sessions as rows (only completed)

```{r partial-sessions-completed}
d.sessions.completed <- d.sessions %>% filter(task_completion == "Completed")

d.sessions.completed
```

### Only completed

```{r partial-submitted}
d.completed <- d %>% filter(task_completion == "Completed")

d.completed$work_experience_programming.s = scale(d.completed$work_experience_programming)
d.completed$work_experience_java.s = scale(d.completed$work_experience_java)
d.completed$time.s = scale(d.completed$time)
d.completed$sonarqube_issues.s = scale(d.completed$sonarqube_issues)

d.completed
```

### Only both completed
```{r partial-submitted-both}
d.both_completed <- d %>% semi_join(d.sessions.completed, by = "session")

d.both_completed$work_experience_programming.s = scale(d.both_completed$work_experience_programming)
d.both_completed$work_experience_java.s = scale(d.both_completed$work_experience_java)
d.both_completed$time.s = scale(d.both_completed$time)
d.both_completed$sonarqube_issues.s = scale(d.both_completed$sonarqube_issues)

d.both_completed
```

## Model expansion utility function

The function `extendable_model` takes some basic arguments for creating brms models and returns a function that can be called with additioanl parameters to combine with those passed to `extendable_model`. The `extendable_model` takes the following arguments:

* `base_name` is a name that is used to identify this extendable model while caching.
* `base_formula` the formula that will be extended and passed to `brms::brm`, represented as string.
* `data` the data frame to be passed to `brms::brm`.
* `base_priors` (`NULL`) is a vector of priors to be passed to `brms::brm`.
* `base_control` (`NULL`) is a vector of control options to be passed to `brms::brm`.

The returned function takes the following arguments:

* `additional_variables` (`NULL`) a vector of aditional variables (predictors) to pass to pass to `brms::brm` in adition to `base_formula`.
* `additional_priors` (`NULL`) a vector of additioanl priors to pass to `brms::brm` in adition to `base_priors`.
* `only_priors` (`FALSE`) indicates if the model should be epty and not compiled, usefull to extract default priors of a model.
* `sample_prior` (`"no"`) is passed to the `sample_prior` of `brms::brm`.
* `control_override` (`NULL`) takes a vector of `control` arguments for `brms::brm` that will override `base_control`.

```{r extendable-model}
extendable_model <- function(
  base_name, 
  base_formula, 
  family, 
  data, 
  base_priors = NULL, 
  base_control = NULL
) {
  function(
    additional_variables = NULL, 
    additional_priors = NULL, 
    only_priors = FALSE, 
    sample_prior = "no", 
    control_override = NULL
  ) {
    # Sort variable names for consistent caching and naming
    additional_variables.sorted <- sort(additional_variables)
    
    # Build priors
    priors <- base_priors
    if (!is.null(additional_priors)) {
      priors <- c(base_priors, additional_priors)
    }
    if (only_priors) {
      priors <- NULL
    }
    
    # Build formula
    additional_variables.formula <- paste(additional_variables.sorted, collapse = " + ")
    formula <- base_formula
    if (!is.null(additional_variables)) {
      formula <- paste(base_formula, additional_variables.formula, sep = " + ")
    }
    
    # Build cache file name
    additional_variables.name <- paste(additional_variables.sorted, collapse = ".")
    name <- base_name
    if (!is.null(additional_variables)) {
      name <- paste(base_name, hash(additional_variables.name), sep = ".")
    }
    name <- paste(name, paste("sample_priors-", sample_prior, sep = ""), sep = ".")
    name <- paste(name, paste("priors_hash-", hash(priors), sep = ""), sep = ".")
    name <- paste(name, paste("formula_hash-", hash(formula), sep = ""), sep = ".")
    
    # Get control options
    control <- base_control
    if (!is.null(control_override)) {
      control <- control_override
    }
    
    # Create and return the brms model
    brm(
      formula = as.formula(formula),
      family = family,
      data = as.data.frame(data),
      prior = priors,
      empty = only_priors,
      sample_prior = sample_prior,
      file = paste("fits", name, sep = "/"),
      file_refit = "on_change",
      seed = 20210421,
      control = control
    )
  }
}
```

Example usage:

```{r, eval=FALSE, class.source = 'fold-show'}
## Not run
m.with <- extendable_model(
  base_name = "m", 
  base_formula = "time ~ 1",
  family = negbinomial(), 
  data = d.both_completed, 
  base_priors = c(
    prior(normal(0, 1), class = "Intercept")
  )
)

prior_summary(m.with(only_priors = TRUE))

pp_check(m.with(sample_prior = "only"), nsamples = 200)

summary(m.with())

pp_check(m.with(), nsamples = 200)

loo(
  m.with(),
  m.with("high_debt_version"),
  m.with(c("high_debt_version", "scenario"))
)
## End(Not run)
```