---
title: "MRTool"
author: "Anne-Sophie Kieslinger"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
               

vignette: >
  %\VignetteIndexEntry{MRTool}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
Dieses Tool simuliert basierend auf bestehenden Genotyp-Daten $G$ Daten für MR-Szenarien (Confounder oder Störfaktor $U$, Exposure oder Risikofaktor $X$ und Outcome $Y$). Außerdem benötigt werden  Parameter, die die Beziehungsstärke beschreiben (siehe Bild). 

Beispiel: $U$ wird berechnet aus $\beta_{GX} * G + N(0,\sigma_U)$. Es kann entweder $\beta_{GX}$ und $\sigma_U$ oder die erklärte Varianz $pve$ angegeben werden. Das gleiche gilt für die anderen Beziehungen. 

 ```{r,engine='tikz', echo=FALSE, fig.cap= "Parameter für die Simulation"}
 \begin{tikzpicture}[scale=.3, transform shape]
 	\begin{scope}
 	\tikzstyle{every node} = [circle, fill=gray!30,minimum size=80]
 	\node (g1)[fill=green!30,rectangle] at (-3, -4) {Genvariante G};
 	\node (x)[fill=blue!30] at (-3,0) {Risikofaktor X};
 	\node (g2)[fill=green!30,rectangle] at (3,-4) {Genvariante i.G};
 	\node (y)[fill=red!30] at (3,0){Outcome Y};
 	\node(u) at (0,3) {Stoerfaktor U};
 	\end{scope}
 	\draw[ ->] (g1) -- (x) node[pos=.5,left] {$GX$};
 	\draw[ ->] (g1) -- (y)  node[pos=.1,right] {$GY$};

 	\draw[ ->] (x) -- (y)  node[pos=.5,above] {$XY$};
 	\draw[ ->] (g2) -- (y)  node[pos=.5,right] {$i.GY$};
 	\draw[ ->] (g2) -- (x)  node[pos=.1,left] {$i.GX$};

 	\draw[->] (u) -- (x)  node[pos=.5,left] {$UX$};
 	\draw[->] (u) -- (y)  node[pos=.5,right] {$UY$};
 	\draw[->] (g1) -- (u)  node[pos=.1,right] {$GU$};
 	\draw[->] (g2) -- (u)  node[pos=.1,left] {$i.GU$};
 	
 \end{tikzpicture}
 ```
 

## Reihenfolge der Funktionen

0. (optional) $get\_corrs$
1. $cal\_betas$
2. $sim\_and\_mr$
3. (optional) $transform\_results$

oder: 

0. (optional) $get\_corrs$
1. $cal\_betas$ oder $SetMRParams$
2. $SimulateMRData$
3. $CreateMRInputObject$
4. $MendelianRandomization::mr\_allmethods$
5. weitere Funktionen des MR Pakets


## Szenario mit 1 Genvariante 



```{r}
library(ggplot2)
# Daten vorbereiten --------------------------------------------------
# SNPs erzeugen
# alle SNP Daten müssen zwei Spalten Identifier haben
# alle SNP Namen müssen ':' haben
# V1 und V2 sind Platzhalter für Identifier bei echten Gendaten, hier ist der Wert egal
#TODO
genedata = data.table::data.table(
  V1 = c('id11', 'id21', 'id31'),
  V2 = c('id12', 'id22', 'id32'),
  'snp:1' = c(rep(0, 504 * 0.5 ^ 2),
              rep(2, 504 * 0.5 ^ 2),
              rep(1, 504 * 0.5 * 0.5 * 2)),
  'snp:2' = c(rep(0, 504 * 0.3 ^ 2),
              rep(2, 504 * 0.7 ^ 2),
              rep(1, 504 * 0.3 * 0.7 * 2))
)
genedata
# Minor Allele Frequency
#TODO
mafs = c('snp:1' = 0.5, 'snp:2' = 0.3)
# dieser SNP soll betrachtet werden
#TODO
snp = 'snp:1'
# die MAF von unserem SNP
maf = mafs[snp]

# Parameter (erklärte Varianzen), die für die Simulation verwendet werden sollen
# siehe Graphik oben
#TODO
pves <- list(
  G_X = c(0.005, 0.2),
  G_Y = c(0),
  G_U = 0,
  X_Y =  c(0, 0.045, 0.295),
  U_X = c(0),
  U_Y = c(0)
)
pves

# Parameter-Berechnung -----------------------------------------------
# calculate the simulation coefficients from explained variance
l <- MRTool::cal_betas(pves, maf)
# parameters for simulation
params <- l[[1]]
# adjustierte erklärte Varianzen
pve_grid <- l[[2]]

params[]
pve_grid

# so oft sollen die Daten neu simuliert und MR ausgeführt werden
#TODO
iterations = 10

# Simulation und Berechnung MR --------------------------------------------
#TODO wieviele cores sollen eingesetzt werden
results <-
  parallel::mclapply(
    1:iterations,
    FUN = function(i)
      MRTool::sim_and_mr(SNP_data = genedata[, .SD, .SDcols = c('V1', 'V2', snp)], SNP =
                           snp,  params),
    mc.cores = 1
  )

# Plotting -------------------------------------------------------
transformed = MRTool::transform_results(results, snp, iterations, pve_grid)

mr_res = transformed$res_list
mr_res_snp1 = mr_res[[1]]
# Spalte power: Anteil der p-Werte des MR-Schätzers <0.05 pro Szenario über die iterations
# caus_eff ist Lineare Regression zwischen X und Y in die Richtung der MR, also estimate_X_Y
# für die kausale Richtung und estimate_Y_X für die antikausale MR.
# Ergebnisse der MR sind 'Estimate', 'Std Error', 'P-Value', '95% CI upper' und '95% CI lower'. 
# Im datatable davor sind Ergebnisse der LR und Parameter der Simulation. Am Ende des datatable sind die pves.

```

Es folgen einige Plots, die genutzt werden können, um die Simulationsdaten und das Ergebnis der MR zu veranschaulichen. Sie dienen als Beispiel, wie mit den erzeugten Datensets umgegangen werden kann.

Der nächste Plot zeigt eine Art, das Ergebnis der MR zu visualisieren. In diesem Szenario ohne Confounder sollten MR-Schätzer und lineare Regression identisch sein. Abweichungen entstehen dadurch, dass 10 Iterations zu wenig sind. 
```{r}
ggplot(data = mr_res_snp1,
       aes(x = assigned.effect_X_Y,
           color = factor(G_X))) +
  geom_abline(slope = 1, intercept = 0) +
  geom_hline(yintercept = 0) +
  geom_line(aes(y = Estimate, linetype = 'MR')) +
  geom_point(aes(y = Estimate, shape = 'MR'), size = 3) +
  geom_line(aes(y = caus_eff, linetype = 'LR')) +
  geom_point(aes(y = caus_eff, shape = 'LR'), size = 3) +
  scale_shape_manual(values = c(1, 3, 5), name = 'Y') +
  scale_linetype(name = 'Y') +
  labs(
    title = 'MR-Estimate relative to X_Y ',
    color = 'Instrumenten-\nstärke',
    x = 'Betakoeffizient xy',
    y = 'Schätzer xy'
  )
```

Der nächste Plot zeigt, ob die Methode, die Betakoeffizienten über die erklärte Varianz zu bestimmen, funktioniert. Die beobachtete erklärte Varianz sollte der vorgegebenen entsprechen. Abweichungen entstehen, wenn der SNP nicht im HWE ist oder die Anzahl der Iterations zu gering.
```{r}
ggplot(data = mr_res_snp1,
       aes(x = G_X, color = factor(G_X))) +
  geom_abline(slope = 1, intercept = 0) +
  geom_line(aes(y = adjRsquared_G_X)) +
  geom_point(aes(y = adjRsquared_G_X)) +
  labs(
    title = 'Beobachtete erklärte Varianz GX vs vorgegebene',
    color = 'Instrumenten-\nstärke',
    y = 'GX real',
    x = 'GX vorgegeben'
  )
```

Der nächste Plot zeigt den Median der Standardabweichung des Schätzers der MR für verschiedene Instrumentenstärken.
```{r}
ggplot(data = mr_res_snp1,
       aes(assigned.effect_X_Y,
           `Std Error`,
           color = factor(G_X))) +
  geom_line() +
  geom_point() +
  labs(title = 'SE of MR estimate rel to beta_X_Y',
       color = 'Instrumenten-\nstärke')

```

## Szenario mit 2 Genvarianten

Der Vorgang bei mehreren Genvarianten funktioniert genau so wie bei einer Genvariante. In der Variable pves wird festgelegt, welcher SNP wie auf X und Y wirkt. 

Durch die Variable rev kann festgelegt werden, ob die MR in die kausale oder antikausale Richtung durchgeführt werden soll. 

In diesem Beispiel werden aus dem oben verwendeten Beispielset genedata die zwei SNPs mit der geringsten Korrelation ausgewählt und die Ergebnisse der MR angeschaut für das Szenario, das ein SNP auf X wirkt und der andere auf Y, für verschiedene Instrumentenstärken. 

```{r,fig.width=5,fig.height=5}
# aus dem Set genedata korrelierte SNPs raussuchen
snp_cors <- MRTool::get_corrs(genedata)
snp_cors
snp_cors[cor=='lowest']

#TODO Richtung der MR, reverse=FALSE für kausale Richtung, reverse=TRUE für antikausale Richtung
rev = FALSE

#TODO
pves <- list(
  # list of snps, each pos in vector is one scenario
  G_X = list(c(0.2, 0.005, 0.2, 0.005), c(0, 0, 0, 0)),
  G_Y = list(c(0, 0, 0, 0), c(0.2, 0.005, 0.005, 0.2)),
  G_U = list(c(0, 0, 0, 0), c(0, 0, 0, 0)),
  X_Y =  c(0, 0.045, 0.295),
  U_X = c(0),
  U_Y = c(0)
)

# für jede Korrelationsstufe
t <- lapply(snp_cors[, cor], function(cor_scen) {
  SNP <- snp_cors[cor == cor_scen, .(snp1, snp2)]
  SNP <- unlist(SNP)
  
  mafs <- mafs[order(names(mafs))]
  SNP <- SNP[order(SNP)]
  
  #TODO
  iterations = 10
  
  # alle snp daten müssen zwei Spalten Identifier haben
  SNP_data =  cbind(genedata[, 1:2, with = F], genedata[, ..SNP])
  maf <- mafs[SNP]
  
  # calculate the simulation coefficients from explained variance
  l <- MRTool::cal_betas(pves, maf)
  # parameters for simulation
  params <- l[[1]]
  
  #TODO falls negative Betakoeffizienten genutzt werden sollen
  #params[,U_Y:=-U_Y]
  #params[,U_X:=-U_X]
  
  # grid of explained variances after calculation of coefficients,
  # pves maybe needed adjustment if impossible
  pve_grid <- l[[2]]
  params
  # remember the true IVs for X and Y respectively
  xsnp = paste0('`', SNP[1], '`')
  ysnp = paste0('`', SNP[2], '`')
  # simulate and MR in parallel for number of iterations
  results <-
    parallel::mclapply(
      1:iterations,
      FUN = function(i)
        MRTool::sim_and_mr(SNP_data, SNP,  params, reverse = rev),
      #TODO
      mc.cores = 1
    )
  
  transformed = MRTool::transform_results(
    results = results,
    SNP = SNP,
    iterations = iterations,
    pve_grid = pve_grid
  )
  res_list = transformed$res_list
  
  # SNP Kategorie hinzufügen: Quelle für X oder Y
  lapply(res_list, function(x)
    x[, snp_type := ifelse(unique(snp) == xsnp, 'xsnp', 'ysnp')])
  
  # return
  data.table::rbindlist(lapply(res_list, function(x) {
    test <- split(x, by = 'snp_id')
    total_res_3 <- test[[xsnp]][test[[ysnp]], on = 'Scenario']
  }), use.names = T)
  
})


#Verarbeitung für Plotting --------------------------------------------------
names(t) <- snp_cors[, cor]
res <- data.table::rbindlist(t, use.names = T, idcol = 'corr_scen')
res[, corr_scen := factor(
  corr_scen,
  levels = c('neg_high', 'neg_low', 'lowest', 'pos_low', 'pos_high'),
  ordered = T
)] 

# je nach Richtung der MR soll die richtige Richtung der LR mit der MR verglichen werden
lr_est = ''
if (rev) {
  res[, caus_eff := estimate_Y_X]
  lr_est = 'Lineare Regression YX'
} else{
  res[, caus_eff := estimate_X_Y]
  lr_est = 'Lineare Regression XY'
}

res_2 <- split(res, by = 'snp_type')
```

Der nächste Plot zeigt das Ergebnis der MR in der kausalen Richtung für verschiedene Instrumentenstärken, wenn der richtige SNP als Instrument für X verwendet wird. Die Facetten zeigen die unterschiedlichen Szenarios, wenn die SNPs auf X und Y unterschiedlich stark einwirken.
Die SNPs sind relativ stark korreliert, daher findet die MR nicht den kausalen Effekt, außer wenn der ysnp nur schwach auf Y wirkt. 

```{r,fig.width=5,fig.height=5}

ggplot(data = res_2$xsnp, aes(x = assigned.effect_X_Y,
                              color = corr_scen)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_hline(yintercept = 0) +
  geom_line(aes(y = Estimate, linetype = 'Kausalschätzer der MR')) +
  geom_point(aes(y = Estimate, shape = 'Kausalschätzer der MR'), size =
               3) +
  geom_line(aes(y = caus_eff, linetype = lr_est)) +
  geom_point(aes(y = caus_eff, shape = lr_est), size = 3) +
  scale_shape_manual(values = c(1, 3, 5), name = 'Y') +
  scale_linetype(name = 'Y') +
  scale_color_discrete(
    breaks = c('neg_high', 'neg_low', 'lowest', 'pos_low', 'pos_high'),
    labels = paste0(paste0(paste0(
      c(
        'negativ hoch',
        'negativ niedrig',
        'am niedrigsten',
        'positiv niedrig',
        'positiv hoch'
      )
    ), ' ~ '),
    c(-0.3,-0.1, snp_cors[cor=='lowest',round(cor_val,digits=2)], 0.1, 0.3))
  ) +
  facet_grid(i.G_Y ~ G_X, labeller = label_both) +
  labs(
    title = 'MR-Estimate relative to X_Y with confounders',
    color = 'Korrelationsszenario ',
    caption = paste0('IV: xsnp \nreverse: ', ifelse(rev, yes = 'yes', no =
                                                      'no')),
    x = 'Betakoeffizient XY',
    y = 'Kausalschätzer'
  )+
  theme(legend.position='bottom',legend.box='vertical')
```

Der nächste Plot zeigt das Ergebnis der MR, wenn der falsche SNP als Instrument für X verwendet wird, nämlich der, der auf Y einwirkt. Die Facetten zeigen die Stärken der Beziehungen der SNPs auf X und Y. Der kausale Effekt wird nicht erkannt.

```{r,fig.width=5,fig.height=5}
ggplot(data = res_2$ysnp,
       aes(x = assigned.effect_X_Y,
           color = corr_scen)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_hline(yintercept = 0) +
  geom_line(aes(y = Estimate,
                linetype = 'Kausalschätzer der MR')) +
  geom_point(aes(y = Estimate,
                 shape = 'Kausalschätzer der MR'),
             size = 3) +
  geom_line(aes(y = caus_eff,
                linetype = lr_est)) +
  geom_point(aes(y = caus_eff,
                 shape = lr_est),
             size = 3) +
  scale_shape_manual(values = c(1, 3, 5),
                     name = 'Y') +
  scale_linetype(name = 'Y') +
  scale_color_discrete(
    breaks = c('neg_high', 'neg_low', 'lowest', 'pos_low', 'pos_high'),
    labels = paste0(paste0(paste0(
      c(
        'negativ hoch',
        'negativ niedrig',
        'am niedrigsten',
        'positiv niedrig',
        'positiv hoch'
      )
    ), ' ~ '),
    c(-0.3,-0.1, snp_cors[cor=='lowest',round(cor_val,digits=2)], 0.1, 0.3))
  ) +
  facet_grid(i.G_Y ~ G_X,
             labeller = label_both) +
  labs(
    title = 'MR-Estimate relative to X_Y with confounders',
    color = 'Korrelationsszenario',
    caption = paste0('IV: ysnp \nreverse: ', ifelse(rev, yes = 'yes', no =
                                                      'no')),
    x = 'Betakoeffizient XY',
    y = 'Kausalschätzer'
  )+
  theme(legend.position='bottom',legend.box='vertical')

```
