Init
library(kirkegaard)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
##
##
## Attaching package: 'magrittr'
##
##
## The following object is masked from 'package:purrr':
##
## set_names
##
##
## The following object is masked from 'package:tidyr':
##
## extract
##
##
## Loading required package: weights
##
## Loading required package: Hmisc
##
##
## Attaching package: 'Hmisc'
##
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
##
## The following objects are masked from 'package:base':
##
## format.pval, units
##
##
## Loading required package: assertthat
##
##
## Attaching package: 'assertthat'
##
##
## The following object is masked from 'package:tibble':
##
## has_name
##
##
## Loading required package: psych
##
##
## Attaching package: 'psych'
##
##
## The following object is masked from 'package:Hmisc':
##
## describe
##
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
##
##
## Attaching package: 'kirkegaard'
##
##
## The following object is masked from 'package:psych':
##
## rescale
##
##
## The following object is masked from 'package:assertthat':
##
## are_equal
##
##
## The following object is masked from 'package:purrr':
##
## is_logical
##
##
## The following object is masked from 'package:base':
##
## +
load_packages(
)
theme_set(theme_bw())
options(
digits = 3
)
Analysis
#make a predictor that is more valid for some groups than others
set.seed(1)
error_sds = seq(1, 3, length.out = 10)
group_means = rnorm(10, sd = 1)
group_means[1] = 0
group_n = 1000
#simulate data
d_sim = map_df(1:10, function(group_i) {
tibble(
group = group_i,
true_PGS = rnorm(group_n, mean = group_means[group_i]),
phenotype = true_PGS + rnorm(group_n, sd = 1),
obs_PGS = (true_PGS + rnorm(group_n, sd = error_sds[group_i]))
)
}) %>% mutate(
group = factor(group)
)
#scatter by group
d_sim %>%
ggplot(aes(obs_PGS, phenotype, color = group)) +
geom_point(alpha = 0.2) +
geom_smooth() +
xlab("Observed polygenic score")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
GG_save("PGS_validity_by_group.png")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
#correlations by group
d_sim_means = d_sim %>%
group_by(group) %>%
summarise(
cor_obs_PGS_phenotype = cor(obs_PGS, phenotype),
true_PGS = mean(true_PGS),
phenotype = mean(phenotype),
obs_PGS = mean(obs_PGS)
)
d_sim_means %>%
GG_scatter("obs_PGS", "phenotype", case_names = "group") +
xlab("Observed polygenic score mean") +
ylab("Phenotypic mean")
## `geom_smooth()` using formula = 'y ~ x'
GG_save("group means.png")
## `geom_smooth()` using formula = 'y ~ x'