single-test-analysis.Rmd

---
title: "Floodlight practice effects"
author: "Tim Woelfle"
date: "29.10.2021"
output:
  html_document:
    code_folding: hide
    toc: yes
    toc_depth: 4
    toc_float:
      collapsed: no
      smooth_scroll: no
params:
  test_code: "ips"
  test_metric_code: "correct_responses"
  unit: "e-SDMT: Correct responses"
  unit_n: "patient"
  unit_time: "repetition"
  min_repetitions: 5
  min_weeks: 5
  up_to_date: "2021-08-01"
---

```{r}
rm(list=ls()[ls() != "params"]) # otherwise some variables might be "carried over" from the previous analysis
```

```{r}
suppressPackageStartupMessages(library(data.table)) # fread
suppressPackageStartupMessages(library(parsedate)) # parse_date
suppressPackageStartupMessages(library(dplyr)) # group_by
suppressPackageStartupMessages(library(tidyr)) # pivot_wider
suppressPackageStartupMessages(library(lme4)) # lmer
suppressPackageStartupMessages(library(mgcv)) # gamm
suppressPackageStartupMessages(library(quantreg)) # rq
suppressPackageStartupMessages(library(patchwork)) # plot_layout
suppressPackageStartupMessages(library(gridExtra)) # grid.arrange
suppressPackageStartupMessages(library(ggpubr)) # ggscatter
suppressPackageStartupMessages(library(ggtext)) # geom_text
suppressPackageStartupMessages(library(sjPlot)) # plot_model

# Download from: https://dataset.floodlightopen.com/public-blobs-prod/complete_dataset.csv
data = fread("complete_dataset.csv", data.table=F)

# Prepare dataset
if (params$test_code == "ips_over_ips_baseline") {
  data = data[data$testCode %in% c("ips", "ips_baseline") & !data$participantIsControl,]
} else {
  data = data[data$testCode == params$test_code & !data$participantIsControl,]
}
data$testStartedAt = parse_date(data$testStartedAt)
data = data[data$testStartedAt <= as.POSIXct(params$up_to_date, tz="UTC"),] # only analyse data up to (excluding) params$up_to_date
data = data[!duplicated(data),] # contains true duplicates for some reason (even with the same testResultMetricId)

data = data[order(data$floodlightOpenId, data$testStartedAt),]
rownames(data) = NULL

# Create composite test: ips (e-SDMT) over ips_baseline (digit tapping)
if (params$test_code == "ips_over_ips_baseline") {
  data = data[data$testMetricCode == "correct_responses",]
  
  data$testResultMetricCreatedOn = parse_date(data$testResultMetricCreatedOn)
  
  # 0 result values are discarded
  data = data[data$testResultMetricValue != 0,]
  
  for (i in 1:nrow(data)) {
    if (i == 1) {
      data[i, "block"] = 1
      next;
    }
    lastBlock = data[i-1,"block"]
    if (sum(data$block==lastBlock, na.rm=T) < 2 & data[i, "floodlightOpenId"] == data[i-1, "floodlightOpenId"] & !data[i, "testCode"] %in% data[data$block==lastBlock, "testCode"] & (difftime(data[i, "testStartedAt"], data[i-1, "testStartedAt"], units="secs")<=300 | difftime(data[i, "testResultMetricCreatedOn"], data[i-1, "testResultMetricCreatedOn"], units="secs")<=300)) {
      data[i, "block"] = lastBlock
      # sometimes testStartedAt has distinct values for IPS & IPS Baseline
      data[i, "testStartedAt"] = data[data$block==lastBlock, "testStartedAt"][1]
    } else {
      data[i, "block"] = lastBlock + 1
    }
  }
  
  # The mismatch between these two tables suggests that for some tests, either only IPS Baseline or IPS have been transmitted
  table(table(data$block)) # should be the same as table(table(data$testStartedAt))
  table(data$testCode)
  #View(data[data$block %in% names(table(data$block))[table(data$block)==1],])
  
  # Select only data with both complete IPS & IPS Baseline
  data = data[as.character(data$block) %in% names(table(data$block))[table(data$block)==2], ]
  
  data = as.data.frame(pivot_wider(data, id_cols=c("floodlightOpenId", "participantIsControl", "participantSex", "participantBirthYear", "testStartedAt"), names_from="testCode", values_from="testResultMetricValue"))
  
  colnames(data)[7] = "baseline"
  
  print(ggplot(data, aes(ips, ips/baseline, color=baseline)) +
    geom_point(alpha=0.2) +
    geom_smooth(method="lm", se=T, color="black") +
    ggtitle(paste0("n=", nrow(data), " (", sum(data$ips/data$baseline > 8), " outliers not shown), Pearson R=",round(cor(data$ips, data$ips/data$baseline), 2), ", Spearman ρ=", round(cor(data$ips, data$ips/data$baseline, method="spearman"), 2))) +
    ylim(0,8) + xlab("e-SDMT") + ylab("e-SDMT corrected"))
  
  data$ips = data$ips / data$baseline # becomes "value" later
  data$baseline = NA # becomes "hand_used" later
# For "Finger Pinching" and "Draw A Shape" hand_used has to be determined
} else if (params$test_code %in% c("pinching", "draw_a_shape")) {
  data = data[data$testMetricCode %in% c(params$test_metric_code, "hand_used"),]
  # just one means either "hand" or "successful_pinches" values are missing, remove those
  table(table(data$testStartedAt))
  #View(data[as.character(data$testStartedAt) %in% names(table(data$testStartedAt))[table(data$testStartedAt)==4],])
  data = data[as.character(data$testStartedAt) %in% names(which(table(data$testStartedAt) == 2 | table(data$testStartedAt) == 4)), ]
  data = as.data.frame(pivot_wider(data, id_cols=c("floodlightOpenId", "participantIsControl", "participantSex", "participantBirthYear", "testStartedAt"), names_from="testMetricCode", values_from="testResultMetricValue"))
} else {
  data = data[data$testMetricCode == params$test_metric_code,]
  data$hand_used = NA
  data = data[c("floodlightOpenId", "participantIsControl", "participantSex", "participantBirthYear", "testStartedAt", "testResultMetricValue", "hand_used")]
}

colnames(data) = c("id", "control", "sex", "birthyear", "time", "value", "hand_used")
data$age = year(data$time)-data$birthyear # Estimate age

# 0 result values are discarded
data = data[data$value != 0,]

# Consider ages supposedly younger than 18 (minimum study age) and older than 90 as NA
data$age[data$age < 18 | data$age > 90] = NA

data$id_original = data$id
data$id = paste0(data$id, "_hand", data$hand_used)
```

```{r}
round = function(x, digits=0) sprintf(paste0("%.", digits, "f"), x)
 
signif_p = function(x, digits=1) {
  x = signif(x, digits)
  if (x < 0.001) return("< .001")
  return(paste0("= ", sub("0\\.", ".", as.character(x))))
}

no_x = theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
no_y = theme(axis.text.y=element_blank(), axis.ticks.y=element_blank())

linecolor = "#c71138"

unit_time_xlab = ifelse(params$unit_time == "repetition", "Repetitions", "Weeks")
```

### Participant selection

```{r}
# At least x weeks & repetitions
for (id in unique(data$id)) {
  subset = data$id == id
  n = sum(subset)
  data[subset, "repetition"] = (1:n)-1
  data[subset, "week"] = as.numeric(difftime(data[subset, "time"], data[subset, "time"][1], unit="weeks"))
}

n_orig = nrow(data)
n_patients_orig = length(unique(data$id_original))
n_hands_orig = length(unique(data$id))

participants_summary = data %>% group_by(id) %>% summarise(week=last(week), repetition=n(), .groups="keep")
```

Among the total n=`r nrow(participants_summary)` `r paste0(params$unit_n, "s")` with n=`r nrow(data)` repetitions, the median length of participation is `r round(median(participants_summary$week),1)` weeks (IQR `r round(quantile(participants_summary$week, 0.25),1)`-`r round(quantile(participants_summary$week, 0.75),1)`, range `r round(min(participants_summary$week),1)`-`r round(max(participants_summary$week),1)`) and the median number of repetitions is `r median(participants_summary$repetition)` (IQR `r quantile(participants_summary$repetition, 0.25)`-`r quantile(participants_summary$repetition, 0.75)`, range `r min(participants_summary$repetition)`-`r max(participants_summary$repetition)`).

```{r}
data = data[data$id %in% participants_summary$id[participants_summary$week >= params$min_weeks & participants_summary$repetition >= params$min_repetitions],]

for (id in unique(data$id)) {
  subset = data$id == id
  n = sum(subset)
  data[subset, "daysSinceLast"] = as.numeric(difftime(data[subset, "time"], c(data[subset, "time"][1], data[subset, "time"][1:(n-1)]), unit="days"))
}

participants_summary = as.data.frame(data %>% group_by(id) %>% summarise(sex=first(sex), mean_age=mean(age), week=last(week), repetition=n(), median_intertest_interval=median(daysSinceLast), IQR_intertest_interval=IQR(daysSinceLast), first=first(value), fifth=nth(value, 5), last=last(value), .groups="keep") %>% mutate(diff_first_last=last-first, diff_first_fifth=fifth-first, diff_fifth_last=last-fifth))

participants_summary$unit_time = participants_summary[, params$unit_time]

data$unit_time = data[,params$unit_time]
```

Inclusion criteria: at least `r params$min_weeks` weeks between first and last repetition, at least `r params$min_repetitions` repetitions performed per test, leading to the analysis of n=`r length(unique(data$id_original))` / `r n_patients_orig` patients `r if (params$unit_n == "hand") paste0(", ", length(unique(data$id)), " / ", n_hands_orig, " hands ")` and n=`r nrow(data)` / `r n_orig` tests.

```{r}
t(data.frame(
  n_patients = paste0(length(unique(data$id_original)), "/", n_patients_orig, " (", round(length(unique(data$id_original))/n_patients_orig*100,1), "%)"),
  n_hands = paste0(length(unique(data$id)), "/", n_hands_orig, " (", round(length(unique(data$id))/n_hands_orig*100,1), "%)"),
  n_tests = paste0(nrow(data), "/", n_orig, " (", round(nrow(data)/n_orig*100,1), "%)"),
  percent_female = paste0(sum(participants_summary$sex == "female"), "/", nrow(participants_summary), " (", round(prop.table(table(participants_summary$sex == "female"))[[2]]*100, 1), "%)"),
  age = paste0(round(median(participants_summary$mean_age,na.rm=T),1), " (", round(quantile(participants_summary$mean_age, 0.25, na.rm=T),1), "-", round(quantile(participants_summary$mean_age, 0.75, na.rm=T),1), ", range ", round(min(participants_summary$mean_age, na.rm=T),1), "-", round(max(participants_summary$mean_age, na.rm=T),1), ")"),
  repetitions = paste0(median(participants_summary$repetition), " repetitions (IQR ", quantile(participants_summary$repetition, 0.25), "-", quantile(participants_summary$repetition, 0.75), ", range ", min(participants_summary$repetition), "-", max(participants_summary$repetition), ")"),
  median_intertest_interval = paste0(round(median(participants_summary$median_intertest_interval),1), " days (IQR ", round(quantile(participants_summary$median_intertest_interval, 0.25),1), "-", round(quantile(participants_summary$median_intertest_interval, 0.75),1), ", range ", round(min(participants_summary$median_intertest_interval),1), "-", round(max(participants_summary$median_intertest_interval),1), ")"),
  IQR_intertest_interval = paste0(round(median(participants_summary$IQR_intertest_interval),1), " days (IQR ", round(quantile(participants_summary$IQR_intertest_interval, 0.25),1), "-", round(quantile(participants_summary$IQR_intertest_interval, 0.75),1), ", range ", round(min(participants_summary$IQR_intertest_interval),1), "-", round(max(participants_summary$IQR_intertest_interval),1), ")"),
  weeks = paste0(round(median(participants_summary$week),1), " weeks (IQR ", round(quantile(participants_summary$week, 0.25),1), "-", round(quantile(participants_summary$week, 0.75),1), ", range ", round(min(participants_summary$week),1), "-", round(max(participants_summary$week),1), ")")
)) 
```

### Summary level analysis

#### Difference test

```{r}
summary(aov(value~name,pivot_longer(participants_summary[c("first", "fifth", "last")], c("first", "fifth", "last"))))

test_first_last = t.test(participants_summary$last, participants_summary$first, paired=T)
test_first_last

print("Average observed difference (last-first) / first: ")
improvement_first_last = paste0(round(test_first_last$estimate/mean(participants_summary$first)*100, 1), "% (", round(test_first_last$conf[1]/mean(participants_summary$first)*100, 1), "-", round(test_first_last$conf[2]/mean(participants_summary$first)*100, 1), ")")
print(improvement_first_last)

test_first_fifth = t.test(participants_summary$fifth, participants_summary$first, paired=T)
test_first_fifth

print("Average observed improvement (fifth-first) / first: ")
improvement_first_fifth = paste0(round(test_first_fifth$estimate/mean(participants_summary$first)*100, 1), "% (", round(test_first_fifth$conf[1]/mean(participants_summary$first)*100, 1), "-", round(test_first_fifth$conf[2]/mean(participants_summary$first)*100, 1), ")")
print(improvement_first_fifth)

test_fifth_last = t.test(participants_summary$last, participants_summary$fifth, paired=T)
test_fifth_last

print("Average observed improvement (last-fifth) / first: ")
improvement_fifth_last = paste0(round(test_fifth_last$estimate/mean(participants_summary$first)*100, 1), "% (", round(test_fifth_last$conf[1]/mean(participants_summary$first)*100, 1), "-", round(test_fifth_last$conf[2]/mean(participants_summary$first)*100, 1), ")")
improvement_fifth_last

# mod_first_last = lm(diff_first_last ~ I(mean_age/10) + I(first/10) + log10(unit_time), participants_summary)
# confint(mod_first_last)
# summ_first_last = summary(mod_first_last)
# summ_first_last
# 
# mod_first_fifth = lm(diff_first_fifth ~ I(mean_age/10) + I(first/10) + log10(unit_time), participants_summary)
# confint(mod_first_fifth)
# summ_first_fifth = summary(mod_first_fifth)
# summ_first_fifth

mod_fifth_last = lm(diff_fifth_last ~ I(mean_age/10) + I(first/10) + I(fifth/10) + log10(unit_time), participants_summary)
summ_fifth_last = summary(mod_fifth_last)
summ_fifth_last
confint(mod_fifth_last)
```

#### Difference plot

```{r, fig.width=9, fig.height=3}
#lab.y = 1.1*max(mean(participants_summary$fifth), mean(participants_summary$last))

p1 = ggbarplot(data.frame(Timepoint=rep(c("First","Fifth","Last"), each=nrow(participants_summary)), value=c(participants_summary$first,participants_summary$fifth,participants_summary$last)), "Timepoint", "value", add="mean_se", label=T, lab.nb.digits=1, lab.vjust=1.9, ylab=params$unit) + xlab("trial") #+ stat_pvalue_manual(data.frame(group1="Fifth", group2="Last", label=paste0("p ", signif_p(test_fifth_last$p.value)), y.position=lab.y), label="label") + scale_y_continuous(expand=expansion(mult=c(0,0.1))) #+ ggtitle(paste0("ANOVA p=", signif(summary(aov(value~name,pivot_longer(participants_summary[c("first", "fifth", "last")], c("first", "fifth", "last"))))[[1]]["Pr(>F)"][1,],1)))

#p2 = plot_model(summ_first_last, show.values=T, show.intercept=F, colors=linecolor, title=paste0("Difference from First to Last Score, R²=", round(summ_first_last$r.squared, 2)), axis.labels=rev(c("Intercept", "Age (per 10 years)", "First score (per 10)", paste0(unit_time_xlab, " (log 10)"))), value.offset=0.3, show.p=F) + ylab("β estimates") + ylab(NULL) + geom_hline(yintercept=0, alpha=0.5)

p3 = plot_model(summ_fifth_last, show.values=T, show.intercept=F, colors=linecolor, title=paste0("Difference from fifth to last score, R²=", round(summ_fifth_last$r.squared, 2)), axis.labels=rev(c("Age (per 10 years)", "First score (per 10)", "Fifth score (per 10)", paste0(unit_time_xlab, " (log 10)"))), value.offset=0.3, show.p=F) + ylab("β estimates") + geom_hline(yintercept=0, alpha=0.5)

#(p1 + (p2/p3)) + plot_layout(widths=c(2,5)) & theme_pubr(base_family="Serif")
(p1 + p3) + plot_layout(widths=c(2,5)) & theme_pubr(base_family="Serif")
```

#### Confounders

```{r, fig.width=15, fig.height=15, warning=FALSE, message=FALSE}
p_age_first = ggscatter(participants_summary, "mean_age", "first", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Mean age") + ylab("First score") + theme_pubr(base_family="Serif")

p_age_fifth = ggscatter(participants_summary, "mean_age", "fifth", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Mean age") + ylab("Fifth score") + theme_pubr(base_family="Serif")

p_age_pred = ggscatter(participants_summary, "mean_age", "unit_time", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_y_log10() + xlab("Mean age") + ylab(paste0(unit_time_xlab, " (log10)")) + theme_pubr(base_family="Serif")

p_age_last = ggscatter(participants_summary, "mean_age", "last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Mean age") + ylab("Last score") + theme_pubr(base_family="Serif")

p_age_diff = ggscatter(participants_summary, "mean_age", "diff_fifth_last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Mean age") + ylab("Difference fifth to last score") + theme_pubr(base_family="Serif")

p_first_fifth = ggscatter(participants_summary, "first", "fifth", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("First score") + ylab("Fifth score") + theme_pubr(base_family="Serif") #+ geom_abline(intercept=0,slope=1)

p_first_pred = ggscatter(participants_summary, "first", "unit_time",  add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_y_log10() + xlab("First score") + ylab(paste0(unit_time_xlab, " (log10)")) + theme_pubr(base_family="Serif")

p_first_last = ggscatter(participants_summary, "first", "last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("First score") + ylab("Last score") + theme_pubr(base_family="Serif") #+ geom_abline(intercept=0,slope=1)

p_first_diff = ggscatter(participants_summary, "first", "diff_fifth_last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("First score") + ylab("Difference fifth to last score") + theme_pubr(base_family="Serif") #+ geom_abline(intercept=0,slope=1)

p_fifth_pred = ggscatter(participants_summary, "fifth", "unit_time", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_y_log10() + xlab("Fifth score") + ylab(paste0(unit_time_xlab, " (log10)")) + theme_pubr(base_family="Serif")

p_fifth_last = ggscatter(participants_summary, "fifth", "last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Fifth score") + ylab("Last score") + theme_pubr(base_family="Serif")

p_fifth_diff = ggscatter(participants_summary, "fifth", "diff_fifth_last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Fifth score") + ylab("Difference fifth to last score") + theme_pubr(base_family="Serif")

p_pred_last = ggscatter(participants_summary, "unit_time", "last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_x_log10(expand = expansion(mult = c(.05, .15))) + xlab(paste0(unit_time_xlab, " (log10)")) + ylab("Last score") + theme_pubr(base_family="Serif")

p_pred_diff = ggscatter(participants_summary, "unit_time", "diff_fifth_last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_x_log10(expand = expansion(mult = c(.05, .15))) + xlab(paste0(unit_time_xlab, " (log10)")) + ylab("Difference fifth to last score") + theme_pubr(base_family="Serif")

p_last_diff = ggscatter(participants_summary, "last", "diff_fifth_last", add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + xlab("Last score") + ylab("Difference fifth to last score") + theme_pubr(base_family="Serif") #+ geom_abline(intercept=0,slope=1)

p_last_pred = ggscatter(participants_summary, "last", "unit_time",  add="reg.line", alpha=0.2, cor.coef=T, cor.coeff.args=list(color=linecolor), conf.int=T, add.params=list(color=linecolor)) + scale_y_log10() + xlab("Last score") + ylab(paste0(unit_time_xlab, " (log10)")) + theme_pubr(base_family="Serif")


p_age = gghistogram(participants_summary, "mean_age", bins=15) + xlab(NULL) + theme_pubr(base_family="Serif")
p_first = gghistogram(participants_summary, "first", bins=15) + xlab(NULL) + theme_pubr(base_family="Serif")
p_fifth = gghistogram(participants_summary, "fifth", bins=15) + xlab(NULL) + theme_pubr(base_family="Serif")
p_last = gghistogram(participants_summary, "last", bins=15) + xlab(NULL) + theme_pubr(base_family="Serif")
p_pred = gghistogram(participants_summary, "unit_time", bins=15) + scale_x_log10() + xlab(NULL) + theme_pubr(base_family="Serif")
p_diff = gghistogram(participants_summary, "diff_fifth_last", bins=15) + xlab("Difference fifth to last") + theme_pubr(base_family="Serif")

m <- matrix(NA, 6, 6)
m[lower.tri(m, diag = T)] <- 1:21
grid.arrange(grobs=list(
  p_age, p_age_first+xlab(NULL), p_age_fifth+xlab(NULL), p_age_pred+xlab(NULL), p_age_last+xlab(NULL), p_age_diff,
  p_first, p_first_fifth+xlab(NULL)+ylab(""), p_first_pred+xlab(NULL)+ylab(""), p_first_last+xlab(NULL)+ylab(""), p_first_diff+ylab(""),
  p_fifth, p_fifth_pred+xlab(NULL)+ylab(""), p_fifth_last+xlab(NULL)+ylab(""), p_fifth_diff+ylab(""),
  p_pred, p_pred_last+xlab(NULL)+ylab(""), p_pred_diff+ylab(""),
  p_last, p_last_diff+ylab(""), 
  p_diff
), layout_matrix=m, heights=c(1,1,1,1,1,1.1))


#GGally::ggpairs(participants_summary[c("mean_age", "first", "fifth", "unit_time", "last", "diff_fifth_last")])
#pairs(participants_summary[c("mean_age", "first", "fifth", "unit_time", "last", "diff_fifth_last")], upper.panel=NULL)
#corrplot::corrplot(cor(participants_summary[c("mean_age", "first", "fifth", "unit_time", "last", "diff_fifth_last")], use="complete.obs"))
```

### Quantile regression

```{r fig.width=8.2, fig.height=4}
data_censored = data[data$repetition < params$min_repetitions,]

percentiles = c(0.05,0.25,0.5,0.75,0.95)

QR = rq(value ~ repetition, tau=percentiles, data_censored)

QR_summ = summary(QR, se="ker")

p_vals = sapply(1:length(QR_summ), function(i) {
  summ = coef(QR_summ[[i]])
  print(summ)
  print(paste0("Intercept: ", round(summ[1,1],1), " (", round(summ[1,1]-1.96*summ[1,2],1), "-", round(summ[1,1]+1.96*summ[1,2],1), "), beta: ", round(summ[2,1],2), " (", round(summ[2,1]-1.96*summ[2,2],2), "-", round(summ[2,1]+1.96*summ[2,2],2), ")"))
  summ[2,4]
})

p_vals = p.adjust(p_vals, method="bonferroni")

try({
  ANOVA = anova(QR)
})

quant = quantile(table(data_censored$id))
print(paste0("n tests = ", nrow(data_censored), " (median tests per ", params$unit_n, ": ", quant["50%"], ", IQR ", quant["25%"], "-", quant["75%"], ")"))

#data_censored$repetition = data_censored$repetition+1
p = ggplot() + geom_line(aes(repetition, value, group=id), data_censored, alpha=0.2, color="darkgrey") +
  theme_pubr(base_family="Serif") + theme(legend.position = "none") +
  scale_x_continuous(expand = expansion(mult = c(0, 0)), labels=function(breaks) breaks+1) + scale_y_continuous(limits=c(0, quantile(data$value, probs=0.975)), expand = expansion(mult = c(0, 0))) +
  xlab("Trial") + ylab(params$unit) +
  geom_abline(intercept=coef(QR)[1,], slope=coef(QR)[2,], color=linecolor) +
  geom_richtext(data=data.frame(intercept=coef(QR)[1,], label=paste0(percentiles*100, "th percentile: β = ", round(coef(QR)[2,],2), ", ", "adjusted *P* ", sapply(p_vals,signif_p))),
            mapping=aes(x=min(data_censored$repetition)+0.2,y=intercept, label=label), color=linecolor, hjust="left", vjust=1, family="Serif", fill = NA, label.color = NA, label.padding = grid::unit(rep(0, 4), "pt"))

if (exists("ANOVA")) {
  p + geom_richtext(aes(x=x,y=y,label=label), data.frame(x=0.7*params$min_repetitions, y=0, label=paste0("ANOVA *P* ", signif_p(ANOVA$table$pvalue, 1))), vjust=-1.5, family="Serif", fill = NA, label.color = NA, label.padding = grid::unit(rep(0, 4), "pt"))
} else {
  p
}
```

```{r}
if (test_fifth_last$p.value > 0.05) {  
  knitr::knit_exit()
}
```

### Learning curve: Model selection

```{r}
# Linear
REM_linear = lmer(value ~ (1|id) + unit_time, data)
equ_linear = function(t) fixef(REM_linear)[1] + fixef(REM_linear)[2]*t
summary(REM_linear)

# Quadratic
REM_quadratic = lmer(value ~ (1|id) + unit_time + I(unit_time^2), data)
equ_quadratic = function(t) fixef(REM_quadratic)[1] + fixef(REM_quadratic)[2]*t + fixef(REM_quadratic)[3]*t^2
summary(REM_quadratic)

# Smoothing spline
smoothing_spline = gamm(value ~ s(unit_time, bs="ps"), random=list(id=~1), data=data)
summary(smoothing_spline$lme)
summary(smoothing_spline$gam)

# Bounded growth
try({
  #summary(nls(value ~ SSasymp(unit_time, yf, y0, log_alpha), data = data, start=c(yf=mean(participants_summary$last), y0=mean(participants_summary$first), log_alpha=-3)))
  REM_bounded = nlmer(value ~ SSasymp(unit_time, yf, y0, log_alpha) ~ (y0|id) + (yf|id), data = data, start=c(yf=mean(participants_summary$last), y0=mean(participants_summary$first), log_alpha=-2), control=nlmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=200000)), verbose=0)
  y0=fixef(REM_bounded)["y0"]
  yf=fixef(REM_bounded)["yf"]
  log_alpha=fixef(REM_bounded)["log_alpha"]
  equ_bounded = function(t, yf=fixef(REM_bounded)[["yf"]], y0=fixef(REM_bounded)[["y0"]], log_alpha=fixef(REM_bounded)[["log_alpha"]]) yf+(y0-yf)*exp(-exp(log_alpha)*t)
  print(summary(REM_bounded))
  print(exp(log_alpha))
  
  print(paste0("Average boundary improvement over baseline: ", (yf-y0)/y0*100))
})

# RMSE vs eDF
models = list(REM_linear, REM_quadratic, smoothing_spline$lme)
if (exists("REM_bounded")) models = c(models, REM_bounded)
RMSE = round(sapply(models, function(mod) sqrt(mean(resid(mod)^2))), 2)
models[[3]] = smoothing_spline$gam
eDF = round(sapply(models, function(mod) { nrow(data)-df.residual(mod) }), 1)
```

#### Comparison plot

```{r fig.width=8.2}
remaining_participants_bins = seq(0,350,by=10)
remaining_participants = data.frame(x=remaining_participants_bins, text=sapply(remaining_participants_bins, function(x) sum(participants_summary$unit_time >= x)))

xmax = remaining_participants$x[which(remaining_participants$text<=20)[1]]

if (xmax > 200) remaining_participants = remaining_participants[seq(1,36,by=2),]

range_95p =quantile(data$value, probs=c(0.025,0.975))

labeltext = paste0("Model RMSE (eDF):<br><span style='color:#0000ff'>····· Linear: ", RMSE[1], " (", eDF[1], ") </span><br><span style='color:#008b00'>- - - Quadratic: ", RMSE[2], " (", eDF[2], ") </span><br><span style='color:#000000'>—— Smoothing spline: ", RMSE[3], " (", eDF[3], ") </span>")
if (exists("REM_bounded")) labeltext = paste0(labeltext, "<br><span style='color:", linecolor, "'>— Bounded growth: ", RMSE[4], " (", eDF[4], ") </span>")

p1 = ggplot() +
  geom_line(aes_string("unit_time", "value", group="id"), data, color="darkgrey", alpha=0.1) +
  geom_line(aes(x,y), data.frame(x=0:xmax, y=predict(smoothing_spline$gam, newdata=data.frame(unit_time=0:xmax))), linetype="F1", size=1) +
  coord_cartesian(xlim=c(0,xmax), ylim=c(range_95p[1], range_95p[2])) + # +xlim(0,xmax) + ylim(range_95p[1], range_95p[2]) +
  stat_function(fun=equ_linear, color="blue", linetype="dotted", size=1) + 
  stat_function(fun=equ_quadratic, color="green4", linetype="dashed", size=1) + 
  theme_pubr(base_family="Serif") + no_x + xlab(NULL) + ylab(params$unit) +
  geom_richtext(aes(x,y,label=label,hjust=1), data.frame(x=0.8*xmax, y=range_95p[2]-(range_95p[2]-range_95p[1])/1.3, label=labeltext), family="Serif")

if (exists("REM_bounded")) p1 = p1 + stat_function(fun=equ_bounded, color=linecolor, size=1)

p2 = ggplot(remaining_participants[remaining_participants$x<=xmax,]) +
  geom_text(aes(x=x,y="A",label=text), family="Serif") +
  theme_pubr(base_family="Serif") + no_y + xlab(unit_time_xlab) + ylab(paste0("Remaining \n ", params$unit_n, "s")) +
  scale_x_continuous(breaks=seq(0,xmax,by=ifelse(xmax > 200, 20, 10)))

(p1 / p2) + plot_layout(heights=c(0.9,0.1))
```

```{r}
if (exists("REM_bounded")) {
  conf = confint(REM_bounded, c("y0","yf","log_alpha"), method="Wald", level=0.95)
  print(conf)
  print(exp(conf[3,]))
  print(paste0("95% CI (", round((conf[1,1]-conf[2,1])/conf[2,1]*100, 2), "-", round((conf[1,2]-conf[2,2])/conf[2,2]*100, 2), ")"))
} 
```

### Bounded growth model

```{r fig.width=8.2, fig.height=4}
if (exists("REM_bounded")) { 
  xmax = remaining_participants$x[which(remaining_participants$text<20)[1]]
  
  equ_diff_REM_bounded = function(t) exp(log_alpha)*(yf-y0)*exp(exp(log_alpha)*-t)
  equ_diff_get_time_REM_bounded = function(target_slope) log(exp(log_alpha)*(yf-y0)/target_slope)/exp(log_alpha)
  
  equ_bounded_get_x = function(target_value) exp(-log_alpha)*log((yf-y0)/(yf-target_value))
  
  #print(equ_bounded(4))
  
  growth_percentiles = c(0, 0.5, 0.9)
  names_percentiles = c("baseline", "50% practice", "90% practice")
  selected_timepoints = equ_bounded_get_x(y0+(yf-y0)*growth_percentiles)
  labelled_points = data.frame(
    x=selected_timepoints,
    y=equ_bounded(selected_timepoints),
    label=paste0("y=", round(equ_bounded(selected_timepoints),1), ", m=", signif(equ_diff_REM_bounded(selected_timepoints),2), " at ", params$unit_time, " ", round(selected_timepoints,0), ",", ifelse(selected_timepoints>=xmax/2, "\n", " "), names_percentiles),
    vjust=1.3,
    color=linecolor
  )
  labelled_points = rbind(labelled_points, list(x=0.87*xmax, y=yf, label=paste0("boundary: ", round(yf, 1)), vjust=-1.0, color=linecolor))
  
  quant = quantile(table(data$id))
  print(paste0("n tests = ", nrow(data), " (n ", params$unit_n, "s = ", length(unique(data$id)), ", median tests per ", params$unit_n, ": ", quant["50%"], ", IQR ", quant["25%"], "-", quant["75%"], ")"))
  
  p1 = ggplot() + geom_line(aes_string("unit_time", "value", group="id"), data, color="darkgrey", alpha=0.2) +
    theme_pubr(base_family="Serif") +
    coord_cartesian(xlim=c(0,xmax), ylim=c(range_95p[1], range_95p[2])) +
    scale_x_continuous(expand = expansion(mult = c(0,0.05))) + scale_y_continuous(expand = expansion(mult = 0)) +
    xlab(unit_time_xlab) + ylab(params$unit) +
    stat_function(fun=equ_bounded, color=linecolor, size=1) +
    geom_point(data=labelled_points[1:(nrow(labelled_points)-1),], aes(x,y,color=I(color)), size=5) +
    geom_text(data=labelled_points, aes(x,y,label=label, vjust=vjust, color=I(color)), hjust=-0.01, family="Serif")
  
  if (exists("conf")) {
    ymin = equ_bounded(seq(0,xmax*1.05,0.05), conf["yf",1], conf["y0",1], conf["log_alpha",1])
    ymax = equ_bounded(seq(0,xmax*1.05,0.05), conf["yf",2], conf["y0",2], conf["log_alpha",2])
    ribbon = data.frame(x=seq(0,xmax*1.05,0.05), ymin=ymin, ymax=ymax)
    p1 = p1 + geom_ribbon(aes(x=x, ymin=ymax, ymax=ymin), ribbon, fill=linecolor, alpha=0.2)
  }
  
  p1 
} 
```