---
title: The proportional recovery rule redux
subtitle: Arguments for biological and predictive relevance
author: Jeff Goldsmith, Robinson Kundert, Andreas Luft, Cathy Stinear, Winston Byblow,
  Angel Garcia de la Garza, Geert Kwakkel, Tomoko Kitago, John W. Krakauer
output:
  word_document:
    reference_docx: my_docx_styles.docx
    fig_caption: yes
---

```{r setup, include=FALSE}
library(tidyverse)
library(viridis)

library(modelr)
library(patchwork)

knitr::opts_chunk$set(
	echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.width = 8,
  fig.height = 6
)

options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis"
)

scale_colour_discrete = scale_colour_viridis_d
scale_fill_discrete = scale_fill_viridis_d

theme_set(
	theme_minimal() + 
	theme(legend.position = "none") + 
	theme(text = element_text(size = 14)))

source(here::here("source", "bootstrap_helpers.R"))
source(here::here("source", "ceiling_effect.R"))
source(here::here("source", "cor_xd.R"))
source(here::here("source", "simulate_data.R"))
source(here::here("source", "fit_const.R"))
source(here::here("source", "clustering_helpers.R"))
source(here::here("source", "fit_models.R"))
source(here::here("source", "simulate_data.R"))

real_dfs = 
	read_csv(here::here("data", "combined_dfs.csv")) 
```

Abstract
--------

This document produces figures and tables included in the paper. The output is a Word document with all figures and their captions; high res versions of figures are exported to `results/`. 

This document also includes some code chunks that don't produce figures but are otherwise relevant to the paper. This includes empirical confirmations of some derived results, and the gap-statistic based analysis used to compare clustering results for Winters to a null of random recovery. 


CorPlot Base
------------

```{r, cor_plot_base}
cor_df = 
	expand.grid(
	  k = c(.01, seq(.1, 1.5, by = .1), 1.7, 2, 3, 100),
	  cor_xy = seq(-1, 1, by = 0.01)) %>% 
	mutate(cor_xd = map2_dbl(k, cor_xy, cor_xd))

cor_plot_base = 
	cor_df %>% 
	ggplot(aes(x = cor_xy, y = cor_xd)) + 
	geom_path(aes(group = k, color = log(k)), alpha = .3) +
	xlab("cor(x, y)") + 
	ylab("cor(x, \u03b4)") + 
  geom_path(
  	data = filter(cor_df, k == 1), 
  	aes(group = k, color = log(k)), alpha = 1, size = 1.2) +
	guides(color = guide_colorbar(barwidth = 0.5, barheight = 10)) 

ggsave(
	here::here("results", "cor_plot_base.png"),
	cor_plot_base, width = 4, height = 4, dpi = 300)

expand.grid(
	  k = 1,
	  cor_xy = seq(0, .9, by = 0.1)) %>% 
	mutate(
		cor_xd = map2_dbl(k, cor_xy, cor_xd),
	  cor_xd = round(cor_xd, 3)) %>% 
	write_csv(here::here("results", "null_vals.csv"))
```

Simulated datasets
------------------

```{r}
sim_dfs = 
	tibble(
		name = LETTERS[1:4],
		mu_x = c(30, 30, 30, 30),
		mu_y = c(30, 30, 53, 53),
		var_x = 196,
		var_ratio = c(1, 1, .16, .16), 
		cor_xy = c(0, .9, 0, .9)
	) %>% 
	mutate(., data = pmap(., simulate_data)) %>% 
	unnest

obs_values = 
	sim_dfs %>% 
	group_by(name) %>% 
  summarize(cor_xd = cor(x, delta), cor_xy = cor(x, y ), k = var(y) / var(x))
```

```{r sim_plots_create}
ggp_xy = 
	sim_dfs %>% 
	select(-delta) %>% 
	gather(key = visit, value = value, x:y) %>% 
	ggplot(aes(x = visit, y = value, group = subj)) + geom_point() + geom_path() + 
	ylab("Outcome value") + xlab("Visit") + 
	facet_grid(.~name)

ggp_xd = 
	sim_dfs %>% 
	ggplot(aes(x = 66 - x, y = delta)) + geom_point() +
	facet_grid(.~name) + 
	ylab("\u03b4") + xlab("Initial Impairment") +
	stat_smooth(method = "lm", se = FALSE)

ggp_c_xd_xy = 
	cor_plot_base + 
	geom_point(data = obs_values, color = "red", size = 2) + 
	geom_text(data = obs_values, aes(label = name), 
						nudge_x = -0.05, nudge_y = 0.05) + 
	theme(legend.position = "right")

ggp_xy_alt = 
	sim_dfs %>% 
	select(-delta) %>% 
	gather(key = visit, value = value, x:y) %>% 
	ggplot(aes(x = visit, y = value, group = subj)) + geom_point() + geom_path() + 
	facet_wrap(.~name, nrow = 2, ncol = 2)

ggp_ex_data = ({ggp_xy / ggp_xd } | ggp_c_xd_xy) + plot_layout(widths = c(3, 1.5))

ggsave(
	here::here("results", "example_data.png"), 
	ggp_ex_data, width = 10, height = 4, dpi = 300)
```

```{r sim_plots_show, fig.height = 4, fig.cap = "Figure 1: The left panels show four simulated datasets, labeled A, B, C, and D. In the top row, panels show outcome values and baseline (x) and follow-up (y). In the bottom row, panels show change (delta) against initial impairment (66 - x). The right panel shows a contour plot of Equation 1, with contours corresponding to values of the variance ratio k and the contour for k = 1 highlighted. Points on this surface show correlation values obtained for Dataset A through D."}
({ggp_xy / ggp_xd } | ggp_c_xd_xy) + plot_layout(widths = c(3, 1))
```

Not included: verification of computations
------------------------------------------

```{r, include = FALSE}
## the connections between correlation and PRR are verified empirically here

b_and_rsq = function(df) {
	
	df %>% 
		summarize(
		  k = var(y) / var(x),
		  b = 1 - sqrt(k) * cor(x, y),
		  r_sqar = ((cor(x, y) * sqrt(k) - 1) / 
						  sqrt(1  + (k) - 2 * cor(x,y) * sqrt(k )    ))^2,
		  r_sqar_alt = ((-b) / sqrt((k) + 2 * b - 1 ))^2, 
		  r_sqar_alt2 = (cor(x, y - x))^2) %>% 
		select(k, b, r_sqar, r_sqar_alt, r_sqar_alt2)
	
}

sim_dfs %>% 
	mutate(ii = 66 - x) %>% 
	group_by(name) %>% 
	nest() %>% 
  mutate(
  	model = map(data, ~lm(delta ~ 1 + ii, data = .x)),
  	res = map(model, broom::tidy),
  	r_sq = map2_dbl(.x = model, .y = data, ~modelr::rsquare(.x, .y)),
  	t = map(data, b_and_rsq)) %>% 
	unnest(t) %>% 
	select(-data, -model) %>% 
	unnest(res) %>% 
	filter(term == "ii") %>% 
	select(name, k, estimate, b, r_sq, r_sqar, r_sqar_alt, r_sqar_alt2) %>% 
	knitr::kable(digits = 3)
```


Simulated data
--------------

```{r, cv_sim, fig.height = 6, fig.cap = "Figure 2: The top panels show three simulated datasets with true outcome values and baseline (x) and follow-up (y); a horizontal line indicates a ceiling on observed values. In the bottom row, panels show the observed (ceiled) value at follow-up against the baseline value."}
sim_dfs = 
	tibble(
		name = c('C[200]','D[200]','E[200]'),
		n = c(200, 200, 200),
		mu_x = c(30, 30, 30),
		mu_y = c(53, 53, 55),
		var_x = 169,
		var_ratio = c(.16, .16, 1), 
		cor_xy = c(0, .9, .9),
		seed = c(2, 2, 14)
	) %>% 
	mutate(., data = pmap(., simulate_data)) %>% 
	unnest %>% 
	mutate(
    fmii = 66 - x,
		y_obs = map_dbl(y, ceiling_effect),
		delta_obs = y_obs - x
	)

model_fits = 
	sim_dfs %>% 
	nest(data = -name) %>% 
	mutate(
		fit = map2(data, data, ~fit_models(train = .x, test = .y, resid = FALSE))
	) %>% 
  select(-data) %>% 
	unnest(fit) %>% 
	select(name, x, starts_with("pred_y_")) %>% 
	pivot_longer(
		starts_with("pred_y"),
		names_to = "model",
		values_to = "y",
		names_prefix = "pred_y_"
	) %>% 
	mutate(model = fct_relevel(model, c("Int", "Const", "GAM", "Exp", "PRR")))
	

data_labs <- c("C~200", "D~200", "E~200")
names(data_labs) <- c("C", "D", "E")

ggp_xy = 
	sim_dfs %>% 
	select(-delta) %>% 
	gather(key = visit, value = value, x:y) %>% 
	ggplot(aes(x = visit, y = value, group = subj)) + 
	geom_point(size = .2) + geom_path(alpha = .2) + 
	geom_hline(aes(yintercept = 66), size = .2, color = "blue") + 
	ylab("True outcome value") + xlab("Visit") + 
	facet_grid(.~name, labeller = label_parsed) 

ggp_yx = 
	sim_dfs %>% 
	ggplot(aes(x = x, y = y_obs)) + geom_point(size = .2) +
	geom_line(data = model_fits, aes(y = y, color = model), size = 1.2) +
	ylab("Observed y") + 
	facet_grid(.~name, labeller = label_parsed) +
	theme(legend.position = "bottom")

ggp_xd = 
	sim_dfs %>% 
	mutate(delta = y_obs - x) %>% 
	ggplot(aes(x = 66 - x, y = delta)) + geom_point(size = .2) +
	ylab("Observed y") + 
	facet_grid(.~name, labeller = label_parsed)

ggp_compare_model = {ggp_xy / ggp_yx }

ggp_compare_model

ggsave(
	here::here("results", "compare_model.png"), 
	ggp_compare_model, width = 9, height = 6, dpi = 300)
```

Cross Validation: Sim data
--------------------------

```{r, pred_cvs, fig.height = 3, cache = TRUE, dependson = "cv_sim", fig.cap = "Figure 3: Each panel shows the distribution of median absolute prediction errors (MAPE) obtained using cross validation for each of three models. Panels correspond to the simulated datasets shown in Figure 2. Models compared are an intercept-only model, a generalized additive model, and the PRR."}
pred_df = 
	sim_dfs %>% 
	mutate(fmii = 66 - x) 

cv_results = 
	pred_df %>% 
	nest(data = -name) %>% 
	mutate(cv_dfs = map(data, crossv_mc, n = 1000)) %>% 
	select(-data) %>% 
	unnest(cv_dfs) %>% 
	mutate(
		test_df = map2(train, test, ~fit_models(train = .x, test = .y)),
		mad = map(test_df, ~summarize_all(.x, list(MAD = ~median(abs(.)))))
		) %>% 
	select(name, .id, mad) %>% 
	unnest(cols = mad) 

cv_results = 
	cv_results %>% 
	pivot_longer(
		contains("MAD"),
		names_to = "method", 
		values_to = "MAD") %>% 
	mutate(
		method = str_remove(method, "resid_"),
		method = str_remove(method, "_MAD"),
		method = fct_relevel(method, c("Int", "Const", "GAM", "Exp", "PRR"))
	)

ggp_compare_model_cv = 
	cv_results %>% 
	ggplot(aes(x = method, y = MAD)) + 
	  geom_violin(adjust = 2) + 
		stat_summary(fun = "median", colour = "red", size = 2, geom = "point") + 
	ylab("MAPE") + xlab("Model") + 
	facet_wrap(~name, labeller = label_parsed)

ggp_compare_model_cv

ggsave(
	here::here("results", "compare_model_cv.png"), 
	ggp_compare_model_cv, width = 9, height = 3, dpi = 300)
```

Reported datasets
-----------------

```{r, real_plot, fig.height = 4, fig.cap = "Figure 4: The left panels show three real datasets. In the top row, panels show outcome values and baseline (x) and follow-up (y); points are colored to indicate recoverers (purple) and non-recoverers (yellow) using the definitions from each paper describing the data. In the bottom row, panels show change (delta) against initial impairment (66 - x), again separating recoverers and non-recoverers. The right panel shows a contour plot of Equation 1, with contours corresponding to values of the variance ratio k and the contour for k = 1 highlighted. Points on this surface show correlation values obtained for the real datasets."}
obs_values = 
	real_dfs %>% 
	filter(severe == FALSE) %>% 
	group_by(name) %>% 
  summarize(cor_xd = cor(x, delta), cor_xy = cor(x, y ), k = var(y) / var(x))

ggp_xy = 
	real_dfs %>% 
	select(-delta) %>% 
	gather(key = visit, value = value, x:y) %>% 
	ggplot(aes(x = visit, y = value, group = subj, color = severe)) + geom_point() + geom_path() + 
	facet_grid(.~name) + 
	ylab("Outcome value") + xlab("Visit")
	theme(legend.position = "none")

ggp_xd = 
	real_dfs %>% 
	ggplot(aes(x = 66 - x, y = delta, color = severe)) + geom_point() +
	facet_grid(.~name) + 
	stat_smooth(
		data = real_dfs %>% filter(severe == FALSE), 
		method = "lm", se = FALSE) + 
	ylab("\u03b4") + xlab("Initial Impairment (66 - x)") +
	theme(legend.position = "none")

ggp_c_xd_xy = 
	cor_plot_base + 
	geom_point(data = obs_values, color = "red", size = 2) + 
	ggrepel::geom_text_repel(
		data = obs_values, aes(label = name)) + 
	theme(legend.position = "right")

({ggp_xy / ggp_xd } | ggp_c_xd_xy) + plot_layout(widths = c(2, 1))

ggp_xy_alt = 
	real_dfs %>% 
	select(-delta) %>% 
	filter(severe == FALSE) %>% 
	gather(key = visit, value = value, x:y) %>% 
	ggplot(aes(x = visit, y = value, group = subj)) + geom_point() + geom_path() + 
	facet_wrap(.~name, nrow = 2, ncol = 2) + 
	theme(legend.position = "none")

ggp_rep_data = ({ggp_xy / ggp_xd } | ggp_c_xd_xy) + plot_layout(widths = c(2, 1))

ggsave(
	here::here("results", "reported_data.png"), 
	ggp_rep_data, width = 9, height = 4, dpi = 300)

ggp_reg_plot = ggp_xy / ggp_xd

ggsave(here::here("results", "regression.png"),
			 ggp_reg_plot, width = 8, height = 4, dpi = 300)
```

Bootstrap: real data
--------------------

```{r, fig.height = 3, fig.cap = "Figure 5: Each panel shows the results of the bootstrap procedure used to obtain inferences about the value of correlations and variance ratio. Red points are the values obtained for the full dataset, and blue points are values obtained in each of 1000 bootstrap samples; points are overlaid on the contour plot of Equation 1."}
boot_analysis = 
	real_dfs %>% 
	filter(severe == FALSE) %>% 
	nest(-name) %>% 
	mutate(boot = map(data, boot_dist, n_boot = 1000)) %>% 
	select(-data)

ggp_boot_analysis = 
	cor_plot_base + 
	geom_point(
		data = boot_analysis %>% select(name, boot) %>% unnest, 
		aes(x = c_xy, y = c_xd), alpha = .4, size = .4, color = "blue") +
	geom_point(data = obs_values, color = "red", size = 2) + 
	facet_grid(~name) + 
	theme(legend.position = "right")

ggp_boot_analysis 

ggsave(here::here("results", "boot_analysis.png"),
			 ggp_boot_analysis, width = 8, height = 3, dpi = 300)
```

```{r}
ci_boot = 
	boot_analysis %>% 
	mutate(ci = map(boot, ci_from_resamp)) %>% 
	select(name, ci) %>% 
	unnest(cols = ci) 

left_join(obs_values, ci_boot) %>% 
	select(-contains("_xd")) %>% 
	mutate_if(is.numeric, round, digits = 2) %>% 
	mutate(
		cor_xy = str_c(cor_xy, " [", c_xy_lb, ", ", c_xy_ub, "]"),
		k      = str_c(k     , " [", k_lb,    ", ", k_ub,    "]")) %>% 
	select(name, k, cor_xy) %>% 
	knitr::kable()
```            


Cross validation: real data
---------------------------

```{r, real_cvs, fig.height = 3, cache = TRUE, fig.cap = "Figure 6: Each panel shows the distribution of median absolute prediction errors (MAPE) obtained using cross validation for each of three models. Panels correspond to the datasets shown in Figure 4. Models compared are an intercept-only model, a generalized additive model, and the PRR."}
pred_df = 
	real_dfs %>% 
	mutate(fmii = 66 - x) %>% 
	filter(severe == FALSE) %>% 
	rename(y_obs = y, delta_obs = delta)

set.seed(11)

cv_results = 
	pred_df %>% 
	nest(data = -name) %>% 
	mutate(cv_dfs = map(data, crossv_mc, n = 1000)) %>% 
	select(-data) %>% 
	unnest(cv_dfs) %>% 
	mutate(
		test_df = map2(train, test, ~fit_models(train = .x, test = .y)),
		mad = map(test_df, ~summarize_at(.x, vars(starts_with("resid")), list(MAD = ~median(abs(.)))))
		) 

ggp_cv_res = 
	cv_results %>% 
	select(name, .id, mad) %>% 
	unnest(cols = mad) %>% 
	pivot_longer(
		contains("MAD"),
		names_to = "method", 
		values_to = "MAD") %>% 
	mutate(
		method = str_remove(method, "resid_"),
		method = str_remove(method, "_MAD"),
		method = fct_relevel(method, c("Int", "Const", "GAM", "Exp", "PRR"))
	) %>%  
	ggplot(aes(x = method, y = MAD)) + 
	  geom_violin(adjust = 2) + 
		stat_summary(fun.y = "median", colour = "red", size = 2, geom = "point") + 
	ylab("MAPE") + xlab("Model") +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.25, hjust = .5)) + 
	facet_wrap(~name)

ggp_cv_res

ggsave(here::here("results", "cross_valid_mape.png"),
			 ggp_cv_res, width = 8, height = 3, dpi = 300)


ggp_cv_fits = 
	cv_results %>% 
	select(.id, name, test_df) %>% 
	unnest(test_df) %>% 
	select(.id, name, x, starts_with("pred")) %>% 
	pivot_longer(
		cols = starts_with("pred"),
		names_to = "method", 
		values_to = "pred", 
		names_prefix = "pred_y_"
	) %>% 
	filter(method %in% c("Const", "Exp", "GAM", "PRR")) %>% 
	mutate(
		unique_id = str_c(.id, method),
		method = fct_relevel(method, c("Const", "GAM", "Exp", "PRR"))) %>% 
	ggplot(aes(x = x, y = pred, group = unique_id, color = method)) + 
	geom_point(data = pred_df, aes(y = y_obs, group = NULL, color = NULL), size = .3, alpha = .3) +
	geom_line(alpha = .05) + 
	facet_grid(name ~ method)

ggsave(here::here("results", "cross_valid_fits.png"),
			 ggp_cv_fits, width = 8, height = 6, dpi = 300)
```

MAPE: real data
---------------

```{r}
real_dfs %>% 
	mutate(fmii = 66 - x) %>% 
	filter(severe == FALSE) %>% 
	rename(y_obs = y, delta_obs = delta) %>% 
	nest(data = -name) %>% 
	mutate(
		test_df = map(data, ~fit_models(train = .x, test = .x)),
	) %>% 
	unnest() %>% 
	group_by(name) %>% 
	summarize(
		v_y = var(y_obs),
		v_resid = var(resid_PRR),
		R2 = (1 - v_resid / v_y) * 100
	) %>% 
	select(name, R2) %>% 
	knitr::kable(digits = 1)
```

Clustering
----------

```{r}
winters_df = 
	real_dfs %>% 
	filter(name == "Winters", x != 66) %>% 
	mutate(
		ii = 66 - x,
		delta = y  - x,
		prop = delta / (66 - x))

winters_df$cluster = hclust_fun(winters_df %>% select(x, y), k = 2)$cluster

random_df = 
	random_recov_data(1, winters_df) %>% 
	mutate(name = "Random Recovery")

random_df$cluster = 3 - hclust_fun(random_df %>% select(x, y), k = 2)$cluster

clust_df = 
	bind_rows(winters_df, random_df) %>% 
	mutate(cluster = str_c("Group ", cluster))
```

```{r, fig.height = 4, fig.cap = "Figure 7: Both panels show change between baseline and follow-up against initial impairment (66 - baseline). Left panel shows data simulated under a 'random recovery' process, in which outcome values are drawn from a uniform distribution over the baseline value and ceiling. Right panels shows data from Winters et al."}
ggp_postclust = 
	clust_df %>% 
	ggplot(aes(x = 66 - x, y = delta, color = cluster)) + 
	geom_point() + 
	geom_smooth(data = clust_df %>% filter(cluster == "Group 1"), method = "lm", se = FALSE) +
	facet_grid(~name) + 
	ylab("\u03b4") + xlab("Initial Impairment") + 
	theme(legend.position = "none")

ggp_postclust 

ggsave(here::here("results", "postcluster.png"),
	ggp_postclust, width = 8, height = 4, dpi = 300)
```


```{r, eval = FALSE}
random_recov_gaps = 
	tibble(seed = 1:1000) %>% 
	mutate(
		random_df = map(seed, random_recov_data, start_df = winters_df),
		clust_df = map(random_df, ~select(.x, x, y)),
		gap_fit = map(
			.x = clust_df, ~cluster::clusGap(.x, FUN = hclust_fun, K.max = 5, B = 50, verbose = FALSE, spaceH0 = "original")),
		gap_stat = map(gap_fit, gap_proc))

winters_gap = 
	select(winters_df, x, y) %>% 
	cluster::clusGap(FUN = hclust_fun, K.max = 5, B = 50, verbose = FALSE, spaceH0 = "original") %>% 
	gap_proc() %>% 
	mutate(seed = -100)

random_recov_gaps %>%
	select(seed, gap_stat) %>% 
	unnest(gap_stat) %>% 
	ggplot(aes(x = k, y = gap, group = seed)) + 
  geom_path() +
	geom_path(data = winters_gap, color = "red")

random_recov_gaps %>%
	select(seed, gap_stat) %>% 
	unnest(gap_stat) %>% 
	filter(k == 2) %>% 
	ggplot(aes(x = gap)) + geom_histogram() + 
  geom_histogram(data = winters_gap %>% filter(k == 2), fill = "red")
```

