---
title: "Protoype for metaanalysis"
author: "Daniel Falster & Alistair Senior"
date: "03/07/2018"
output: 
  html_document:
    fig_height: 6
    fig_width: 10
    df_print: paged
    rows.print: 10
    code_folding: show
    theme: yeti
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE, echo=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, cache=FALSE)
knitr::opts_knit$set(root.dir =  rprojroot::find_root("README.md"))
```

```{r}
library(readr)
library(dplyr)
library(metafor)

source("R/meta_analysis.R")
source("R/data_load_clean.R")
source("R/calc_pop_stats.R")
```

# Setup 


Load data & clean data

```{r, eval = FALSE}
# Load raw data - save cleaned dataset as RDS for reuse
data_raw <- load_raw("data/dr7.0_all_control_data.csv") %>%
    clean_raw_data()
dir.create("export", F, F)
saveRDS(data_raw, "export/data_clean.rds")
```

```{r}
data <- readRDS("export/data_clean.rds")
```

Define population variable
```{r}
data <- data %>%
  mutate(population = sprintf("%s-%s", production_center, strain_name))
```

Create a copy of data for variance analysis
```{r}
data_par_age <- data_subset_parameter_individual_by_age(data, "body weight", age_min = 98, age_center = 100)
```

Now calculate statistics by population:

```{r}
population_stats <- calculate_population_stats(data_par_age)
```

Now analyse by population data:


Sanity checks

```{r}
# i <- seq(1, nrow(population_stats), by = 2)
# all(population_stats$sex[i] == "female")
# all(population_stats$sex[i + 1] == "male")
# all(population_stats$strain_name[i + 1] == population_stats$strain_name[i])
# all(population_stats$production_center[i + 1] == population_stats$production_center[i])
```
Check plot of mean vs variance

```{r}
plot(population_stats$x_bar, population_stats$x_sd, pch=16, log="xy", col = as.factor(population_stats$sex))
```

```{r}
results <- create_meta_analysis_effect_sizes(population_stats)
```

## Compare CVR

```{r}
plot(results$effect_size_CVR, 1/sqrt(results$sample_variance_CVR), ylab="sample precision", xlab="effect size", pch=16, xlim=c(-0.7, 0.7), col=as.factor(results$production_center))
abline(v=0)
```

```{r}
metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1| strain_name, ~1|production_center, ~1|err), data = results)
```

## Compare VR

```{r}
plot(results$effect_size_VR, 1/sqrt(results$sample_variance_VR), ylab="sample precision", xlab="effect size", pch=16, xlim=c(-0.8, 0.8), col=as.factor(results$production_center))
abline(v=0)
```

```{r}
metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1| strain_name, ~1|production_center, ~1|err), data = results)
```


## Compare Means

```{r}
plot(results$effect_size_RR, 1/sqrt(results$sample_variance_RR), ylab="sample precision", xlab="effect size", pch=16, xlim=c(-0.5, 0.5), col=as.factor(results$production_center))
abline(v=0)
```

```{r}
metafor::rma.mv(yi = effect_size_RR, V = sample_variance_RR, random = list(~1| strain_name, ~1|production_center, ~1|err), data = results)
```
