CTC_MR_combined
Claire Morrison
2023-02-13
to plot:
thoughts…
- plot correlations (see corr plot script)
- plot between/within associations - separate for ABCD and LTS
- plot sig MR associations - if a significant IVW, plot other methods AND related traits
- however we finalize the co-twin and MR plot together
results<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/ctc_results.csv", header=T, data.table=F)
MR_out<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sleep_exposure_MR.csv", header=T, data.table=F)
results$OR<- exp(results$beta)
results$OR_lo<- NA
results$OR_hi<- NA
head(MR_out)## id.exposure id.outcome outcome
## 1 lVP89E ieu-b-102 Major depression || id:ieu-b-102
## 2 lVP89E ieu-b-102 Major depression || id:ieu-b-102
## 3 lVP89E ieu-b-102 Major depression || id:ieu-b-102
## 4 lVP89E ieu-b-102 Major depression || id:ieu-b-102
## 5 lVP89E ieu-b-102 Major depression || id:ieu-b-102
## 6 lVP89E ieu-b-42 schizophrenia || id:ieu-b-42
## exposure method nsnp b
## 1 Sleep duration (unit increase) Maximum likelihood 44 -0.003790686
## 2 Sleep duration (unit increase) MR Egger 44 -0.003997578
## 3 Sleep duration (unit increase) Weighted mode 44 -0.002878508
## 4 Sleep duration (unit increase) Weighted median 44 -0.003027430
## 5 Sleep duration (unit increase) Inverse variance weighted 44 -0.003700711
## 6 Sleep duration (unit increase) Maximum likelihood 44 0.006757675
## se pval lo_ci up_ci or or_lci95
## 1 0.0008569348 9.709624e-06 -0.005470278 -0.0021110938 0.9962165 0.9945447
## 2 0.0043827653 3.669151e-01 -0.012587798 0.0045926424 0.9960104 0.9874911
## 3 0.0020044580 1.582228e-01 -0.006807246 0.0010502292 0.9971256 0.9932159
## 4 0.0013214918 2.196815e-02 -0.005617554 -0.0004373057 0.9969771 0.9943982
## 5 0.0013220044 5.121062e-03 -0.006291840 -0.0011095826 0.9963061 0.9937279
## 6 0.0021850579 1.983609e-03 0.002474961 0.0110403886 1.0067806 1.0024780
## or_uci95
## 1 0.9978911
## 2 1.0046032
## 3 1.0010508
## 4 0.9995628
## 5 0.9988910
## 6 1.0111016table(MR_out$outcome)##
## ADHD || id:ieu-a-1183 Anti Social Behavior
## 55 50
## Anxiety Bipolar Disorder
## 55 55
## Cigarettes per Day Drinks per Week
## 55 55
## Major depression || id:ieu-b-102 schizophrenia || id:ieu-b-42
## 55 55table(MR_out$exposure)##
## Accelerometer sleep duration Alertness Factor
## 40 40
## Chronotype Circadian Preference Factor
## 40 40
## Daytime alertness (sleepiness GWAS) Efficiency
## 40 35
## Efficiency Factor No napping (napping GWAS)
## 40 40
## Non-Insomnia Sleep duration (unit increase)
## 40 40
## Sleep Episodes
## 40table(MR_out$method)##
## Inverse variance weighted Maximum likelihood MR Egger
## 87 87 87
## Weighted median Weighted mode
## 87 87MR_out<- MR_out %>% filter(!method=="Outlier-corrected")
MR_out$exposure[MR_out$exposure=="Insomnia symptoms (never/rarely vs. sometimes/usually)"] <- "Insomnia"
MR_out$exposure[MR_out$exposure=="Daytime nap || id:ebi-a-GCST011494"] <- "Napping"
MR_out$exposure[MR_out$exposure=="Sleep duration (unit increase)"] <- "Self-reported sleep duration"
MR_out$outcome[MR_out$outcome=="Efficiency"] <- "Accelerometer Efficiency"
MR_out$outcome[MR_out$outcome=="ADHD || id:ieu-a-1183"] <- "ADHD"
MR_out$outcome[MR_out$outcome=="schizophrenia || id:ieu-b-42"] <- "Schizophrenia"
MR_out$outcome[MR_out$outcome=="bipolar disorder"] <- "Bipolar Disorder"
MR_out$outcome[MR_out$outcome=="cigarettes per day"] <- "Cigarettes per Day"
MR_out$outcome[MR_out$outcome=="Major depression || id:ieu-b-102"] <- "Depression"
MR_out$id.exposure<- NULL
MR_out$id.outcome<- NULL
MR_out$lo_ci<- NULL
MR_out$up_ci<- NULL
colnames(MR_out)[1:10]<- c("psychiatric_trait_specific", "sleep_trait_specific", "model", "n_SNPs", "beta", "beta_se", "beta_p", "OR", "OR_lo", "OR_hi")
MR_out$sample<- "MR"
MR_out$sleep[MR_out$sleep_trait_specific=="Accelerometer sleep duration" | MR_out$sleep_trait_specific=="Self-reported sleep duration"] <- "Duration"
MR_out$sleep[MR_out$sleep_trait_specific=="No napping (napping GWAS)" | MR_out$sleep_trait_specific=="Alertness Factor" | MR_out$sleep_trait_specific=="Daytime alertness (sleepiness GWAS)"] <- "Alertness"
MR_out$sleep[MR_out$sleep_trait_specific=="Chronotype"| MR_out$sleep_trait_specific=="Circadian Preference Factor"] <- "Chronotype"
MR_out$sleep[MR_out$sleep_trait_specific=="Efficiency" | MR_out$sleep_trait_specific=="Efficiency Factor" | MR_out$sleep_trait_specific=="Sleep Episodes"] <- "Efficiency"
MR_out$psychiatric[MR_out$psychiatric_trait_specific=="Depression" | MR_out$psychiatric_trait_specific=="Anxiety"] <- "Internalizing"
MR_out$psychiatric[MR_out$psychiatric_trait_specific=="Schizophrenia" | MR_out$psychiatric_trait_specific=="Bipolar Disorder"] <- "Psychosis"
MR_out$psychiatric[MR_out$psychiatric_trait_specific=="ADHD"] <- "Attention Problems"
MR_out$psychiatric[MR_out$psychiatric_trait_specific=="Anti Social Behavior" | MR_out$psychiatric_trait_specific=="Cigarettes per Day" | MR_out$psychiatric_trait_specific=="Drinks per Week"] <- "Externalizing"
MR_out<- MR_out %>% select(colnames(results))
results<- rbind(results, MR_out)
results$sample_psych<- ifelse(is.na(results$psychiatric_trait_specific), paste(results$sample, results$psychiatric, sep="_"), paste(results$sleep_trait_specific, results$psychiatric_trait_specific, sep="_"))co twin within results
within<- results[grep("Within", results$model),]
ctc_p<- .05/4
# chrono
within %>% select(sleep, sleep_trait_specific,psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Chronotype") %>%
filter(beta_p<ctc_p )## sleep sleep_trait_specific psychiatric beta beta_p model
## 1 Chronotype Chronotype Externalizing -0.16354 0.00751 Within
## 2 Chronotype Chronotype Externalizing -0.42456 0.01120 Within sib
## 3 Chronotype Social Jetlag Externalizing 0.17500 0.00704 Within
## 4 Chronotype Social Jetlag Externalizing 0.47200 0.00560 Within sib
## sample
## 1 ABCD
## 2 ABCD
## 3 ABCD
## 4 ABCD# alertness
within %>% select(sleep, psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Alertness") %>%
filter(beta_p<ctc_p )## [1] sleep psychiatric beta beta_p model sample
## <0 rows> (or 0-length row.names)# effic
within %>% select(sleep, sleep_trait_specific,psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Efficiency") %>%
filter(beta_p<ctc_p )## sleep sleep_trait_specific psychiatric beta beta_p
## 1 Efficiency Accelerometer Weekday Efficiency Attention Problems -0.24 0.002000
## 2 Efficiency Accelerometer Weekday Efficiency Externalizing -0.24 0.000279
## model sample
## 1 Within DZ ABCD
## 2 Within DZ ABCD# duration
within %>% select(sleep, sleep_trait_specific,psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Duration") %>%
filter(beta_p<ctc_p )## sleep sleep_trait_specific psychiatric beta
## 1 Duration Weekday Duration Psychosis -0.158859
## 2 Duration Accelerometer Weekday Duration Attention Problems -0.251620
## 3 Duration Weekday Duration Psychosis -0.282930
## 4 Duration Weekday Duration Psychosis -0.277147
## beta_p model sample
## 1 1.100e-02 Within MZ ABCD
## 2 5.247e-03 Within DZ ABCD
## 3 2.090e-05 Within ABCD
## 4 1.030e-07 Within DZ ABCD# satisfaction
within %>% select(sleep, psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Satisfaction") %>%
filter(beta_p<ctc_p )## [1] sleep psychiatric beta beta_p model sample
## <0 rows> (or 0-length row.names)# insom
within %>% select(sleep, psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Insomnia") %>%
filter(beta_p<ctc_p )## sleep psychiatric beta beta_p model sample
## 1 Insomnia Internalizing 0.12 0.012 Within Colorado# var
within %>% select(sleep, psychiatric, beta, beta_p, model, sample) %>%
filter(sleep=="Variability") %>%
filter(beta_p<ctc_p )## [1] sleep psychiatric beta beta_p model sample
## <0 rows> (or 0-length row.names)MR IVW results
ivw<- results[grep("Inverse variance weighted", results$model),]
ivw_p<- .05/8
# chrono
ivw %>% select(sleep, sleep_trait_specific, psychiatric_trait_specific, OR, beta_p) %>%
filter(sleep=="Chronotype") %>%
filter(beta_p<ivw_p )## sleep sleep_trait_specific psychiatric_trait_specific OR
## 1 Chronotype Circadian Preference Factor Drinks per Week 1.121713
## beta_p
## 1 0.0001020805results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Circadian Preference Factor") %>%
filter(psychiatric_trait_specific=="Drinks per Week")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific
## 1 Circadian Preference Factor Drinks per Week
## 2 Circadian Preference Factor Drinks per Week
## 3 Circadian Preference Factor Drinks per Week
## model OR beta_p
## 1 Maximum likelihood 1.123158 2.818292e-08
## 2 Weighted median 1.103489 2.394274e-03
## 3 Inverse variance weighted 1.121713 1.020805e-04# alertness
ivw %>% select(sleep, sleep_trait_specific, psychiatric_trait_specific, OR, beta_p) %>%
filter(sleep=="Alertness") %>%
filter(beta_p<ivw_p)## sleep sleep_trait_specific psychiatric_trait_specific OR
## 1 Alertness No napping (napping GWAS) Depression 0.6815710
## 2 Alertness No napping (napping GWAS) Schizophrenia 0.4540639
## 3 Alertness No napping (napping GWAS) Bipolar Disorder 0.5329535
## beta_p
## 1 5.280796e-07
## 2 3.198915e-04
## 3 1.511661e-03results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="No napping (napping GWAS)") %>%
filter(psychiatric_trait_specific=="Depression")## sleep_trait_specific psychiatric_trait_specific
## 1 No napping (napping GWAS) Depression
## 2 No napping (napping GWAS) Depression
## 3 No napping (napping GWAS) Depression
## 4 No napping (napping GWAS) Depression
## 5 No napping (napping GWAS) Depression
## model OR beta_p
## 1 Maximum likelihood 0.6738814 3.595328e-12
## 2 MR Egger 1.1054844 6.999369e-01
## 3 Weighted mode 1.0076273 9.617286e-01
## 4 Weighted median 0.8096428 1.965047e-02
## 5 Inverse variance weighted 0.6815710 5.280796e-07results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="No napping (napping GWAS)") %>%
filter(psychiatric_trait_specific=="Schizophrenia") %>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific
## 1 No napping (napping GWAS) Schizophrenia
## 2 No napping (napping GWAS) Schizophrenia
## 3 No napping (napping GWAS) Schizophrenia
## 4 No napping (napping GWAS) Schizophrenia
## model OR beta_p
## 1 Maximum likelihood 0.4536748 3.025806e-08
## 2 Weighted mode 0.2003513 4.482934e-03
## 3 Weighted median 0.4285470 3.347649e-04
## 4 Inverse variance weighted 0.4540639 3.198915e-04results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="No napping (napping GWAS)") %>%
filter(psychiatric_trait_specific=="Bipolar Disorder") %>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific
## 1 No napping (napping GWAS) Bipolar Disorder
## 2 No napping (napping GWAS) Bipolar Disorder
## model OR beta_p
## 1 Maximum likelihood 0.5282921 3.288510e-06
## 2 Inverse variance weighted 0.5329535 1.511661e-03# effic
ivw %>% select(sleep, sleep_trait_specific, psychiatric_trait_specific, OR, beta_p) %>%
filter(sleep=="Efficiency") %>%
filter(beta_p<ivw_p ) ### remove ASB --- unrealistic## sleep sleep_trait_specific psychiatric_trait_specific OR
## 1 Efficiency Efficiency Factor Depression 0.5779504
## beta_p
## 1 0.0002680735results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Efficiency Factor") %>%
filter(psychiatric_trait_specific=="Depression")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific model
## 1 Efficiency Factor Depression Maximum likelihood
## 2 Efficiency Factor Depression Weighted median
## 3 Efficiency Factor Depression Inverse variance weighted
## OR beta_p
## 1 0.5674757 5.514911e-06
## 2 0.6115202 4.701350e-03
## 3 0.5779504 2.680735e-04# dur
ivw %>% select(sleep, sleep_trait_specific, psychiatric_trait_specific, OR,beta, beta_p) %>%
filter(sleep=="Duration") %>%
filter(beta_p<ivw_p ) ### remove ASB## sleep sleep_trait_specific psychiatric_trait_specific OR
## 1 Duration Self-reported sleep duration Depression 0.9963061
## 2 Duration Accelerometer sleep duration Schizophrenia 1.6416609
## 3 Duration Accelerometer sleep duration Anti Social Behavior 1006.7755661
## beta beta_p
## 1 -0.003700711 5.121062e-03
## 2 0.495708484 8.110242e-05
## 3 6.914507994 5.208813e-03results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Self-reported sleep duration") %>%
filter(psychiatric_trait_specific=="Depression")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific
## 1 Self-reported sleep duration Depression
## 2 Self-reported sleep duration Depression
## model OR beta_p
## 1 Maximum likelihood 0.9962165 9.709624e-06
## 2 Inverse variance weighted 0.9963061 5.121062e-03results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Accelerometer sleep duration") %>%
filter(psychiatric_trait_specific=="Schizophrenia")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific
## 1 Accelerometer sleep duration Schizophrenia
## 2 Accelerometer sleep duration Schizophrenia
## 3 Accelerometer sleep duration Schizophrenia
## model OR beta_p
## 1 Maximum likelihood 1.646551 1.093846e-04
## 2 Weighted median 1.550797 5.303916e-03
## 3 Inverse variance weighted 1.641661 8.110242e-05# insom
ivw %>% select(sleep, sleep_trait_specific, psychiatric_trait_specific, OR, beta_p) %>%
filter(sleep=="Non-Insomnia") %>%
filter(beta_p<ivw_p )## sleep sleep_trait_specific psychiatric_trait_specific OR
## 1 Non-Insomnia Non-Insomnia Depression 0.3850484
## 2 Non-Insomnia Non-Insomnia Anxiety 5.6429822
## 3 Non-Insomnia Non-Insomnia ADHD 0.3167471
## 4 Non-Insomnia Non-Insomnia Drinks per Week 1.1283073
## beta_p
## 1 5.366771e-17
## 2 1.964025e-04
## 3 6.654396e-05
## 4 9.708955e-04results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Non-Insomnia") %>%
filter(psychiatric_trait_specific=="Depression")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific model
## 1 Non-Insomnia Depression Maximum likelihood
## 2 Non-Insomnia Depression Weighted mode
## 3 Non-Insomnia Depression Weighted median
## 4 Non-Insomnia Depression Inverse variance weighted
## OR beta_p
## 1 0.3763089 2.571577e-26
## 2 0.3274129 5.073069e-03
## 3 0.4205445 1.163517e-09
## 4 0.3850484 5.366771e-17results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Non-Insomnia") %>%
filter(psychiatric_trait_specific=="Anxiety")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific model
## 1 Non-Insomnia Anxiety Maximum likelihood
## 2 Non-Insomnia Anxiety Inverse variance weighted
## OR beta_p
## 1 6.204112 1.546891e-06
## 2 5.642982 1.964025e-04results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Non-Insomnia") %>%
filter(psychiatric_trait_specific=="ADHD")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific model
## 1 Non-Insomnia ADHD Maximum likelihood
## 2 Non-Insomnia ADHD Inverse variance weighted
## OR beta_p
## 1 0.3328918 2.312973e-05
## 2 0.3167471 6.654396e-05results %>% select(sleep_trait_specific, psychiatric_trait_specific, model, OR, beta_p) %>%
filter(sleep_trait_specific=="Non-Insomnia") %>%
filter(psychiatric_trait_specific=="Drinks per Week")%>%
filter(beta_p<ivw_p )## sleep_trait_specific psychiatric_trait_specific model
## 1 Non-Insomnia Drinks per Week Maximum likelihood
## 2 Non-Insomnia Drinks per Week Inverse variance weighted
## OR beta_p
## 1 1.130406 0.0009705446
## 2 1.128307 0.0009708955full across methods results
converging across: - chrono & ext - effic & ext - dur & ext - dur & psych - insom & int - insom & ext
both<- results %>% filter(model=="Inverse variance weighted" | model=="Within" | model=="Within MZ" | model=="Within DZ")
# chrono
both %>% select(sleep, psychiatric, sample,model, beta, beta_p) %>%
filter(sleep=="Chronotype") %>%
filter(beta_p<ctc_p ) %>%
arrange(psychiatric) ### both, but chrono factor## sleep psychiatric sample model beta
## 1 Chronotype Externalizing ABCD Within -0.1635400
## 2 Chronotype Externalizing ABCD Within 0.1750000
## 3 Chronotype Externalizing MR Inverse variance weighted 0.1148572
## beta_p
## 1 0.0075100000
## 2 0.0070400000
## 3 0.0001020805# alertness
both %>% select(sleep, psychiatric, sample, beta, beta_p) %>%
filter(sleep=="Alertness") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric) ### only MR## sleep psychiatric sample beta beta_p
## 1 Alertness Internalizing MR -0.3833548 5.280796e-07
## 2 Alertness Internalizing MR -0.6832640 1.236765e-02
## 3 Alertness Psychosis MR -0.7895173 3.198915e-04
## 4 Alertness Psychosis MR -0.6293212 1.511661e-03# effic
both %>% select(sleep, sleep_trait_specific, psychiatric, psychiatric_trait_specific, model, sample, beta, beta_p) %>%
filter(sleep=="Efficiency") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric) ## only co twin## sleep sleep_trait_specific psychiatric
## 1 Efficiency Accelerometer Weekday Efficiency Attention Problems
## 2 Efficiency Accelerometer Weekday Efficiency Externalizing
## 3 Efficiency Efficiency Factor Internalizing
## psychiatric_trait_specific model sample beta
## 1 <NA> Within DZ ABCD -0.2400000
## 2 <NA> Within DZ ABCD -0.2400000
## 3 Depression Inverse variance weighted MR -0.5482673
## beta_p
## 1 0.0020000000
## 2 0.0002790000
## 3 0.0002680735# dur
both %>% select(sleep, sleep_trait_specific, psychiatric, psychiatric_trait_specific, model, sample, beta, beta_p) %>%
filter(sleep=="Duration") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric) ### all three samples## sleep sleep_trait_specific psychiatric
## 1 Duration Accelerometer Weekday Duration Attention Problems
## 2 Duration Accelerometer sleep duration Externalizing
## 3 Duration Self-reported sleep duration Internalizing
## 4 Duration Weekday Duration Psychosis
## 5 Duration Weekday Duration Psychosis
## 6 Duration Weekday Duration Psychosis
## 7 Duration Accelerometer sleep duration Psychosis
## psychiatric_trait_specific model sample beta
## 1 <NA> Within DZ ABCD -0.251620000
## 2 Anti Social Behavior Inverse variance weighted MR 6.914507994
## 3 Depression Inverse variance weighted MR -0.003700711
## 4 <NA> Within MZ ABCD -0.158859000
## 5 <NA> Within ABCD -0.282930000
## 6 <NA> Within DZ ABCD -0.277147000
## 7 Schizophrenia Inverse variance weighted MR 0.495708484
## beta_p
## 1 5.247000e-03
## 2 5.208813e-03
## 3 5.121062e-03
## 4 1.100000e-02
## 5 2.090000e-05
## 6 1.030000e-07
## 7 8.110242e-05# insom
both %>% select(sleep, sleep_trait_specific, psychiatric, psychiatric_trait_specific, model, sample, beta, beta_p) %>%
filter(sleep=="Insomnia") %>%
filter(beta_p<ctc_p ) ### MR and co twin colorado## sleep sleep_trait_specific psychiatric psychiatric_trait_specific model
## 1 Insomnia Insomnia Internalizing <NA> Within
## sample beta beta_p
## 1 Colorado 0.12 0.012across measurement types
pheno<- results %>% filter(model=="Phenotypic")
# effic
pheno %>% select(sleep, sleep_trait_specific, psychiatric, sample, beta, beta_p) %>%
filter(sleep=="Efficiency") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric)## sleep sleep_trait_specific psychiatric sample beta beta_p
## 1 Efficiency Accelerometer Weekday Efficiency Externalizing ABCD -0.1 0.005# dur
pheno %>% select(sleep, sleep_trait_specific, psychiatric, sample, beta, beta_p) %>%
filter(sleep=="Duration") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric)## sleep sleep_trait_specific psychiatric sample beta
## 1 Duration Weekday Duration Attention Problems ABCD -0.08558
## 2 Duration Weekday Duration Externalizing ABCD -0.08000
## 3 Duration Accelerometer Weekend Duration Psychosis ABCD -0.09000
## 4 Duration Accelerometer Weekday Duration Psychosis ABCD -0.10000
## 5 Duration Weekend Duration Psychosis ABCD -0.08090
## 6 Duration Weekday Duration Psychosis ABCD -0.22699
## beta_p
## 1 3.58e-03
## 2 3.00e-03
## 3 8.00e-03
## 4 3.00e-03
## 5 3.92e-03
## 6 2.00e-16# effic
within %>% select(sleep, sleep_trait_specific, psychiatric, sample, beta, beta_p) %>%
filter(sleep=="Efficiency") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric)## sleep sleep_trait_specific psychiatric sample beta
## 1 Efficiency Accelerometer Weekday Efficiency Attention Problems ABCD -0.24
## 2 Efficiency Accelerometer Weekday Efficiency Externalizing ABCD -0.24
## beta_p
## 1 0.002000
## 2 0.000279# dur
within %>% select(sleep, sleep_trait_specific, psychiatric, sample, beta, beta_p) %>%
filter(sleep=="Duration") %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric)## sleep sleep_trait_specific psychiatric sample beta
## 1 Duration Accelerometer Weekday Duration Attention Problems ABCD -0.251620
## 2 Duration Weekday Duration Psychosis ABCD -0.158859
## 3 Duration Weekday Duration Psychosis ABCD -0.282930
## 4 Duration Weekday Duration Psychosis ABCD -0.277147
## beta_p
## 1 5.247e-03
## 2 1.100e-02
## 3 2.090e-05
## 4 1.030e-07# accel
accel<- both[grep("Accelerometer", both$sleep_trait_specific),]
accel %>% select(sleep, sleep_trait_specific, psychiatric, sample, beta, beta_p) %>%
filter(beta_p<ctc_p )%>%
arrange(psychiatric)## sleep sleep_trait_specific psychiatric sample
## 1 Duration Accelerometer Weekday Duration Attention Problems ABCD
## 2 Efficiency Accelerometer Weekday Efficiency Attention Problems ABCD
## 3 Efficiency Accelerometer Weekday Efficiency Externalizing ABCD
## 4 Duration Accelerometer sleep duration Externalizing MR
## 5 Duration Accelerometer sleep duration Psychosis MR
## beta beta_p
## 1 -0.2516200 5.247000e-03
## 2 -0.2400000 2.000000e-03
## 3 -0.2400000 2.790000e-04
## 4 6.9145080 5.208813e-03
## 5 0.4957085 8.110242e-05plots
set up model colors and order
table(results$sleep) ### need to combine as insomnia... remember MR is non insomnia, CTC is insomnia!##
## Alertness Chronotype Duration Efficiency Insomnia Non-Insomnia
## 130 138 196 163 10 40
## Satisfaction Variability
## 10 24table(results$psychiatric)##
## Attention Problems Externalizing Internalizing Psychosis
## 109 244 194 164results$sleep[results$sleep=="Non-Insomnia"] <- "Insomnia"
phenotypic_col<- "burlywood"
between_col<- "peachpuff"
within_col<- "dodgerblue"
within_MZ_col<- "steelblue2"
within_DZ_col<- "deepskyblue"
within_sib_col<- "skyblue"
IVW_col<- "tan3"
ML_col<- "darkorange"
wmode_col<- "darkgoldenrod2"
wmed_col<- "goldenrod1"
mregger_col<- "gold1"
model_colors<- c("Phenotypic"=phenotypic_col, "Between"=between_col, "Within"=within_col, "Within MZ"=within_MZ_col, "Within DZ"=within_DZ_col, "Within sib"=within_sib_col, "Maximum likelihood"=ML_col, "Inverse variance weighted"=IVW_col, "Weighted mode"= wmode_col, "Weighted median"=wmed_col, "MR Egger" = mregger_col)
results<- results %>%
mutate(model=fct_relevel(model,"Within sib", "Within DZ", "Within MZ", "Within", "Between", "Phenotypic", "Inverse variance weighted", "Maximum likelihood", "Weighted mode", "Weighted median"),
psychiatric= fct_relevel(psychiatric, "Internalizing", "Externalizing", "Attention Problems", "Psychosis"))co twin plots
## ABCD
abcd<- results %>% filter(sample=="ABCD")
abcd_sleep<- unique(abcd$sleep)
abcd_plt<- list()
for (s in abcd_sleep) {
dat<- results %>% filter(sample=="ABCD" & sleep==s) %>%
mutate(CI=1.96*beta_se)
plt<- ggplot(dat, aes(colour=model, shape=sleep_trait_specific,x=beta, y=model, xmin = beta-CI, xmax = beta+CI)) +
geom_point(aes(colour=model), size=7, position=position_dodge(width=0.9)) +
geom_errorbarh(aes(x = CI), height=1, position=position_dodge(width=0.7))+
facet_wrap(~psychiatric) +
xlab("Beta") +
ylab("Coefficient")+
geom_vline(xintercept = 0, linetype="dashed")+
ggtitle(paste0(s, " in ABCD"))+
scale_color_manual(values = model_colors)+
theme(legend.title = element_blank(),
#legend.text=element_text(size=10), couldnt fit all the duration legend items with a bigger text
legend.position="bottom",
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 17),
plot.title = element_text(size=17))
plt<- plt+guides(colour = "none")
abcd_plt[[s]]<- plt
}## Warning in geom_errorbarh(aes(x = CI), height = 1, position = position_dodge(width = 0.7)): Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x## LTS
lts<- results %>% filter(sample=="Colorado")
lts_sleep<- unique(lts$sleep)
lts_plt<- list()
for (s in lts_sleep) {
dat<- results %>% filter(sample=="Colorado" & sleep==s)%>%
mutate(CI=1.96*beta_se)
plt<- ggplot(dat, aes(colour=model, shape=sleep_trait_specific,x=beta, y=model, xmin = beta-CI, xmax = beta+CI)) +
geom_point(size=7,position=position_dodge(width=0.7)) +
geom_errorbarh(aes(x = CI), height=1, position=position_dodge(width=0.7))+
facet_wrap(~psychiatric) +
geom_vline(xintercept = 0, linetype="dashed")+
xlab("Beta") +
ylab("Coefficient")+
ggtitle(paste0(s, " in Colorado Sample"))+
scale_color_manual(values = model_colors)+
theme(legend.title = element_blank(),
legend.position="bottom",
legend.text=element_text(size=15),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 17),
plot.title = element_text(size=17))
plt<- plt+guides(colour = "none")
lts_plt[[s]]<- plt
}## Warning in geom_errorbarh(aes(x = CI), height = 1, position = position_dodge(width = 0.7)): Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: x
## Ignoring unknown aesthetics: xlts<- ggarrange(lts_plt[[3]], lts_plt[[4]], lts_plt[[2]],
lts_plt[[5]], lts_plt[[1]], nrow=3, ncol=2, common.legend = F)
ggsave(plot = lts,filename = "/Users/claire/Desktop/dissertation//figs/LTS_CTC.PNG", width = 15, height = 14, units = "in")
abcd<- ggarrange(abcd_plt[[3]], abcd_plt[[4]], abcd_plt[[2]],
abcd_plt[[1]], nrow=2, ncol=2, common.legend = F)
ggsave(plot = abcd,filename = "/Users/claire/Desktop/dissertation//figs/ABCD_CTC.PNG", width = 18, height = 18, units = "in")MR plots
MR<- results %>% filter(sample=="MR") %>%
mutate(model=fct_relevel(model, "Inverse variance weighted", "Maximum likelihood", "Weighted mode", "Weighted median", "MR Egger"))
MR<- MR %>% filter(OR<10) %>%
filter(psychiatric_trait_specific!="Anti Social Behavior")
head(MR)## sleep_trait_specific sleep psychiatric
## 1 Self-reported sleep duration Duration Internalizing
## 2 Self-reported sleep duration Duration Internalizing
## 3 Self-reported sleep duration Duration Internalizing
## 4 Self-reported sleep duration Duration Internalizing
## 5 Self-reported sleep duration Duration Internalizing
## 6 Self-reported sleep duration Duration Psychosis
## psychiatric_trait_specific beta beta_se beta_p
## 1 Depression -0.003790686 0.0008569348 9.709624e-06
## 2 Depression -0.003997578 0.0043827653 3.669151e-01
## 3 Depression -0.002878508 0.0020044580 1.582228e-01
## 4 Depression -0.003027430 0.0013214918 2.196815e-02
## 5 Depression -0.003700711 0.0013220044 5.121062e-03
## 6 Schizophrenia 0.006757675 0.0021850579 1.983609e-03
## model sample n_SNPs OR OR_lo OR_hi
## 1 Maximum likelihood MR 44 0.9962165 0.9945447 0.9978911
## 2 MR Egger MR 44 0.9960104 0.9874911 1.0046032
## 3 Weighted mode MR 44 0.9971256 0.9932159 1.0010508
## 4 Weighted median MR 44 0.9969771 0.9943982 0.9995628
## 5 Inverse variance weighted MR 44 0.9963061 0.9937279 0.9988910
## 6 Maximum likelihood MR 44 1.0067806 1.0024780 1.0111016
## sample_psych
## 1 Self-reported sleep duration_Depression
## 2 Self-reported sleep duration_Depression
## 3 Self-reported sleep duration_Depression
## 4 Self-reported sleep duration_Depression
## 5 Self-reported sleep duration_Depression
## 6 Self-reported sleep duration_Schizophreniadat<- list()
MR_plt<- list()
plt<- list()
for (s in unique(MR$sleep)) {
dat[[s]]<- MR %>% filter(sleep==s)%>%
mutate(CI=1.96*beta_se)
for (p in unique(MR$psychiatric)) {
dat[[paste0(s,p)]]<- dat[[s]] %>% filter(psychiatric==p)
plt[[p]]<- ggplot(dat[[paste0(s,p)]], aes(shape=sleep_trait_specific, colour=model, x=OR, y=model, xmin = OR, xmax=OR)) +
geom_point(size=8,position=position_dodge(width=0.5)) +
geom_errorbarh(aes(xmin = OR_lo, xmax = OR_hi), height=1.2, position=position_dodge(width=0.5))+
facet_wrap(~psychiatric_trait_specific,scales = "free_x") +
geom_vline(xintercept = 1, linetype="dashed")+
xlab("Coefficient") +
ylab("Beta")+
ggtitle(paste(s, "Predicting", p, "MR", sep=" "))+
scale_color_manual(values = model_colors)+
theme(legend.title = element_blank(),
legend.position="bottom",
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 10),
plot.title = element_text(size=17))
plt[[p]]<- plt[[p]]+guides(colour = "none")
MR_plt[[paste0(s,p)]]<- plt[[p]]
}
}
### dur
dur<- ggarrange(MR_plt[[1]], MR_plt[[4]], MR_plt[[3]], MR_plt[[2]], nrow=2, ncol=2, common.legend = F)## Warning: `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervalsggsave(plot = dur,filename = "/Users/claire/Desktop/dissertation/figs/Duration_MR.png", width = 20, height = 12, units = "in")
### alert
alert<- ggarrange(MR_plt[[5]], MR_plt[[8]], MR_plt[[7]], MR_plt[[6]], nrow=2, ncol=2, common.legend = F)## Warning: `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervalsggsave(plot = alert,filename = "/Users/claire/Desktop/dissertation/figs/Alertness_MR.png", width = 20, height = 12, units = "in")
### chrono
chrono<- ggarrange(MR_plt[[9]], MR_plt[[12]], MR_plt[[11]], MR_plt[[10]], nrow=2, ncol=2, common.legend = F)## Warning: `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervalsggsave(plot = chrono,filename = "/Users/claire/Desktop/dissertation/figs/Chronotype_MR.png", width = 20, height = 12, units = "in")
### insom
insom<- ggarrange(MR_plt[[13]], MR_plt[[16]], MR_plt[[15]], MR_plt[[14]], nrow=2, ncol=2, common.legend = F)## Warning: `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervalsggsave(plot = insom,filename = "/Users/claire/Desktop/dissertation/figs/Insomnia_MR.png", width = 20, height = 12, units = "in")
### effic
effic<- ggarrange(MR_plt[[17]], MR_plt[[20]], MR_plt[[19]], MR_plt[[18]], nrow=2, ncol=2, common.legend = F)## Warning: `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervals
## `position_dodge()` requires non-overlapping x intervalsggsave(plot = effic,filename = "/Users/claire/Desktop/dissertation/figs/Efficiency_MR.png", width = 20, height = 12, units = "in")co twin and MR combined plot
issues with this are the categories aren’t super neat.
technically the MR for satisfaction is napping, and sleep duration is average habitual sleep duration for MR
had to add specific “sample_psych” on X axis because otherwise it was pooling results if more than one significant finding in the same sleep-psych domain
### only sig
MR$beta_sig<- ifelse(MR$beta_p<0.05, MR$beta, NA)
### had already gotten rid of those with huge ORs above ^^
within$beta_sig<- ifelse(within$beta_p<0.05, within$beta, NA)
results_sig<- rbind(MR, within)
table(results_sig$sleep)##
## Alertness Chronotype Duration Efficiency Insomnia Satisfaction
## 111 108 146 136 41 6
## Variability
## 16unique(results_sig$sleep)## [1] "Duration" "Alertness" "Chronotype" "Insomnia" "Efficiency"
## [6] "Satisfaction" "Variability"combo_plt<- list()
for (s in unique(results_sig$sleep)) {
dat<- results_sig %>% filter(sleep==s) %>%
filter(model!="Phenotypic") %>%
filter(model!="Between")
plt<- ggplot(dat, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
xlab(" ") +
ylab("Beta")+
# geom_errorbarh(aes(x = CI), height=.4)+
ggtitle(paste0(s))+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
# axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
axis.text.x = element_blank(),
plot.title = element_text(size=20))
combo_plt[[s]]<- plt
}
combo<- ggarrange(combo_plt[[3]], combo_plt[[2]],
combo_plt[[5]], combo_plt[[1]], combo_plt[[4]],
combo_plt[[6]], combo_plt[[7]], nrow=7, ncol=1, common.legend = T)## Warning: Removed 88 rows containing missing values (`geom_bar()`).
## Removed 88 rows containing missing values (`geom_bar()`).## Warning: Removed 87 rows containing missing values (`geom_bar()`).## Warning: Removed 119 rows containing missing values (`geom_bar()`).## Warning: Removed 128 rows containing missing values (`geom_bar()`).## Warning: Removed 21 rows containing missing values (`geom_bar()`).## Warning: Removed 6 rows containing missing values (`geom_bar()`).## Warning: Removed 16 rows containing missing values (`geom_bar()`).ggsave(plot = combo,filename = "/Users/claire/Desktop/dissertation/figs/MR_CTC.PNG", width = 11, height = 20, units = "in")plotting specific findings for defense
### combo only sig
### only sig
MR$beta_sig<- ifelse(MR$beta_p<ivw_p, MR$beta, NA)
### had already gotten rid of those with huge ORs above ^^
within$beta_sig<- ifelse(within$beta_p<ctc_p, within$beta, NA)
results_sig<- rbind(MR, within)
dat<- results_sig %>% filter(model!="Phenotypic") %>%
filter(model!="Between") %>%
filter(!is.na(beta_sig)) %>%
filter(OR<10)
chrono<- dat %>% filter(sleep=="Chronotype")
c<- ggplot(chrono, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity", width=.5)+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
xlab(" ") +
ylab("Beta")+
ggtitle("Chronotype")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
alert<- dat %>% filter(sleep=="Alertness")
a<- ggplot(alert, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity", width=.5)+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
scale_x_discrete(labels = function(x) str_wrap(x, width = 2))+
xlab(" ") +
ylab("Beta")+
ggtitle("Alertness")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
dur<- dat %>% filter(sleep=="Duration")
d<- ggplot(dur, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity", width=.5)+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
xlab(" ") +
ylab("Beta")+
ggtitle("Sleep Duration")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
effic<- dat %>% filter(sleep=="Efficiency")
e<- ggplot(effic, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
xlab(" ") +
ylab("Beta")+
ggtitle("Efficiency")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
insom<- dat %>% filter(sleep=="Insomnia")
i<- ggplot(insom, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
facet_wrap(~psychiatric, ncol = 4, scales = "free")+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
xlab(" ") +
ylab("Beta")+
ggtitle("Insomnia")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
combo<- ggarrange(c, a, e, d, i, nrow = 5, common.legend = T)
ggsave(plot = combo,filename = "/Users/claire/Desktop/dissertation/figs/defense_sig_CTC_MR.PNG", width = 10, height = 17, units = "in")
### specific assoc.
effic<- effic %>% filter(psychiatric=="Externalizing")
ggplot(effic, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
xlab(" ") +
ylab("Beta")+
ggtitle("Efficiency")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_fill_manual(values = model_colors)+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
dur<- dur %>% filter(psychiatric=="Psychosis")
ggplot(dur, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
xlab(" ") +
ylab("Beta")+
ggtitle("Duration predicting psychosis")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))
insom<- insom %>% filter(psychiatric=="Internalizing")
ggplot(insom, aes(fill=model, y=beta_sig, x=sample_psych)) +
geom_bar(position="dodge", stat="identity")+
xlab(" ") +
ylab("Beta")+
ggtitle("Insomnia predicting Internalizing")+
geom_hline(yintercept = 0, linetype="dashed", linewidth=1)+
scale_x_discrete(labels = function(x) str_wrap(x, width = 5))+
scale_fill_manual(values = model_colors)+
theme(legend.title = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_text(size=10, angle = 90, vjust = 0.5, hjust=1),
plot.title = element_text(size=20))