In askieslinger/MRTool: Simulating Data for Mendelian Randomization Scenarios Based on Genotype Data
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
\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
- (optional) $get_corrs$
- $cal_betas$
- $sim_and_mr$
- (optional) $transform_results$
oder:
- (optional) $get_corrs$
- $cal_betas$ oder $SetMRParams$
- $SimulateMRData$
- $CreateMRInputObject$
- $MendelianRandomization::mr_allmethods$
- weitere Funktionen des MR Pakets
Szenario mit 1 Genvariante
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.
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.
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.
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.
# 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.
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.
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')
askieslinger/MRTool documentation built on April 29, 2020, 12:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 \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
- (optional) $get_corrs$
- $cal_betas$
- $sim_and_mr$
- (optional) $transform_results$
oder:
- (optional) $get_corrs$
- $cal_betas$ oder $SetMRParams$
- $SimulateMRData$
- $CreateMRInputObject$
- $MendelianRandomization::mr_allmethods$
- weitere Funktionen des MR Pakets
Szenario mit 1 Genvariante
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.
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.
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.
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.
# 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.
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.
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')
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.