generated from opensafely/research-template
/
plot_outcomes_byage.R
58 lines (44 loc) · 1.92 KB
/
plot_outcomes_byage.R
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
################################################################
# This script:
# - Calculates number of outcomes by week by age
################################################################
# For running locally only #
# setwd("C:/Users/aschaffer/OneDrive - Nexus365/Documents/GitHub/vax-fourth-dose-RD")
# getwd()
# Import libraries #
library('tidyverse')
library('lubridate')
library('arrow')
library('here')
library('reshape2')
library('dplyr')
library('fs')
library('ggplot2')
library('RColorBrewer')
library('data.table')
## Create directories
dir_create(here::here("output", "covid_outcomes", "by_start_date"), showWarnings = FALSE, recurse = TRUE)
dir_create(here::here("output", "covid_outcomes", "figures"), showWarnings = FALSE, recurse = TRUE)
# Load functions
source(here::here("analysis", "custom_functions.R"))
data <- read.csv(here::here("output", "covid_outcomes", "by_start_date", "outcomes_byage_6mon_2022-11-26_red.csv"))
data2 <- data %>% dplyr::select(c(contains("rate_"), age_mos6)) %>%
reshape2::melt(id = c("age_mos6"))
############################################################
### Plot event rate by age in months and index date
############################################################
ggplot(data2,
aes(x = age_mos6 / 12, y = value)) +
geom_vline(aes(xintercept = 50), linetype = "longdash") +
geom_point() +
scale_y_continuous(expand = expansion(mult = c(0, .2))) +
scale_x_continuous(breaks = seq(47,53,1)) +
facet_wrap(~ variable, ncol = 3, scales = "free_y") +
xlab("Age in months") + ylab("No. events per 100,000") +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.title = element_blank(), legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(here::here("output", "covid_outcomes", "figures", "plot_outcomes.png"),
dpi = 300, units = "in", width = 10, height = 7.5)