/
simulation3_plotresults.Rmd
76 lines (58 loc) · 3.16 KB
/
simulation3_plotresults.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
---
title: "results"
author: "Alan Min"
date: "7/27/2020"
output: html_document
---
```{r, include = F}
get_legend<-function(myggplot){
tmp <- ggplot_gtable(ggplot_build(myggplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)
}
library(tidyverse)
library(gridExtra)
library(cowplot)
knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE)
```
```{r fig.width=12, fig.height = 4}
yaxis = c(.4, 1.2)
# Load results from cousins data
load("results_cousins.Rda")
results = df
# Rename some of the results
results$names = recode(results$names, elizabeth_no_diag = "HE", elizabeth_diag = "HE2")
results = results %>% filter(names != "R_maxlike" & names != "schwartzman" & names != "HE2")
results = results %>% mutate(names = factor(names, levels = c("HE", "dicker1", "dicker2", "gcta", "ldak", "gold_standard")))
# Take the real part of the results. Although this doesn't happen here
# taking eigenvalues can sometimes give an imaginary result
results$result_hsq = Re(results$result_hsq)
results$sigmae = Re(results$sigmae)
results$sigmag = Re(results$sigmag)
# Rename the levels of cousinship
to_plot = results
to_plot$structure = recode(results$structure, cousin0='(iv)', cousin1='(i)', cousin2='(ii)', cousin3='(iii)')
to_plot = to_plot %>% mutate(structure = factor(structure, levels = c("(i)", "(ii)", "(iii)", "(iv)")))
# Manually set colors for the estimators
colorscheme = c("#F8766D", "#BB9D00", "#00B81F", "#000000", "#00A5FF", "#E76BF3")
colorscheme = setNames(colorscheme, c( "HE", "dicker1", "dicker2", "gcta", "ldak", "gold_standard"))
# Plot data for cousins
ggplot(to_plot) + geom_boxplot(aes(x = as.factor(p), y = result_hsq, col = names)) + coord_cartesian(ylim = yaxis) + xlab("Number of Markers") + ylab("Heritability") + facet_grid(rows = vars(structure)) + theme(legend.position = "none", strip.background = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_color_manual(values = colorscheme, drop = T)
ggsave(file = "related_estimates.pdf", height = 8, width = 8)
```
```{r}
# Plot MSE of each of the number of markers for cousinships
plot_data = results %>%
group_by(names, p, structure) %>%
summarise(MSE = sqrt(sum((result_hsq - .8)^2)/length(result_hsq)), .groups = "keep")
names(plot_data)[1] = "estimator"
plot_data$structure = as.character(plot_data$structure)
plot_data$structure[plot_data$structure == "cousin0"] = "unrelated"
plot_data$structure[plot_data$structure == "cousin1"] = "1st Cousins"
plot_data$structure[plot_data$structure == "cousin2"] = "2nd Cousins"
plot_data$structure[plot_data$structure == "cousin3"] = "3rd Cousins"
names(plot_data)[3] = "relatedness"
print(ggplot(plot_data) + geom_line(aes(x = p, y = MSE, col = estimator, linetype = relatedness)) + coord_cartesian(ylim = c(0, .3)) + ggtitle("Number of Markers by Relatedness") + scale_linetype_manual(breaks=c("1st Cousins","2nd Cousins", "3rd Cousins", "unrelated"), values=c(3,2,6,1)) + xlab("Number of Markers") + scale_color_manual(values = colorscheme, drop = T))
ggsave("~/Documents/UW Documents/gcta_proj/Scripts/21_03_01_relatededits20_12_20/mse_by_relatedness.pdf", width = 7, height = 5)
```