simulate_under_LTM simulates families and thresholds under
the liability threshold model for a given family structure and a
variable number of phenotypes.Please note that it is not possible
to simulate different family structures.
Usage
simulate_under_LTM(
fam_vec = c("m", "f", "s1", "mgm", "mgf", "pgm", "pgf"),
n_fam = NULL,
add_ind = TRUE,
h2 = 0.5,
genetic_corrmat = NULL,
full_corrmat = NULL,
phen_names = NULL,
n_sim = 1000,
pop_prev = 0.1
)Arguments
- fam_vec
A vector of strings holding the different family members. All family members must be represented by strings from the following list: -
m(Mother) -f(Father) -c[0-9]*.[0-9]*(Children) -mgm(Maternal grandmother) -mgf(Maternal grandfather) -pgm(Paternal grandmother) -pgf(Paternal grandfather) -s[0-9]*(Full siblings) -mhs[0-9]*(Half-siblings - maternal side) -phs[0-9]*(Half-siblings - paternal side) -mau[0-9]*(Aunts/Uncles - maternal side) -pau[0-9]*(Aunts/Uncles - paternal side). Defaults toc("m","f","s1","mgm","mgf","pgm","pgf").- n_fam
A named vector holding the desired number of family members. See
setNames. All names must be picked from the list mentioned above. Defaults toNULL.- add_ind
A logical scalar indicating whether the genetic component of the full liability as well as the full liability for the underlying target individual should be included in the covariance matrix. Defaults to
TRUE.- h2
Either a number or a numeric vector holding the liability-scale heritability(ies) for one or more phenotypes. All entries in
h2must be non-negative. Note that under the liability threshold model, the heritabilities must also be at most 1. Defaults to 0.5.- genetic_corrmat
Either
NULLor a numeric matrix holding the genetic correlations between the desired phenotypes. Must be specified, iflength(h2)\(>0\), and will be ignored ifh2is a number. All diagonal entries ingenetic_corrmatmust be equal to one, while all off-diagonal entries must be between -1 and 1. In addition, the matrix must be symmetric. Defaults toNULL.- full_corrmat
Either
NULLor a numeric matrix holding the full correlations between the desired phenotypes. Must be specified, iflength(h2)\(>0\), and will be ignored ifh2is a number. All diagonal entries infull_corrmatmust be equal to one, while all off-diagonal entries must be between -1 and 1. In addition, the matrix must be symmetric. Defaults toNULL.- phen_names
Either
NULLor character vector holding the phenotype names. These names will be used to create the row and column names for the covariance matrix. Must be specified, iflength(h2)\(> 0\), and will be ignored ifh2is a number. If it is not specified, the names will default to phenotype1, phenotype2, etc. Defaults toNULL.- n_sim
A positive number representing the number of simulations. Defaults to 1000.
- pop_prev
Either a number or a numeric vector holding the population prevalence(s), i.e. the overall prevalence(s) in the population. All entries in
pop_prevmust be positive and smaller than 1. Defaults to 0.1.
Value
If either fam_vec or n_fam is used as the argument,
if it is of the required format, if the liability-scale heritability h2
is a number satisfying \(0 \leq h^2\), n_sim is a strictly positive number,
and pop_prev is a positive number that is at most one,
then the output will be a list containing two tibbles.
The first tibble, sim_obs, holds the simulated liabilities, the disease
status and the current age/age-of-onset for all family members in each of the
n_sim families.
The second tibble, thresholds, holds the family identifier, the personal
identifier, the role (specified in fam_vec or n_fam) as well as the lower and
upper thresholds for all individuals in all families. Note that this tibble has
the format required in estimate_liability.
If either fam_vec or n_fam is used as the argument and if it is of the
required format, if genetic_corrmat and full_corrmat are two numeric
and symmetric matrices satisfying that all diagonal entries are one and that all
off-diagonal entries are between -1 and 1, if the liability-scale heritabilities in
h2_vec are numbers satisfying \(0 \leq h^2_i\) for all \(i \in \{1,...,n_pheno\}\),
n_sim is a strictly positive number, and pop_prev is a positive numeric
vector such that all entries are at most one, then the output will be a list containing
the following lists.
The first outer list, which is named after the first phenotype in phen_names,
holds the tibble sim_obs, which holds the simulated liabilities, the
disease status and the current age/age-of-onset for all family members in each of
the n_sim families for the first phenotype.
As the first outer list, the second outer list, which is named after the second
phenotype in phen_names, holds the tibble sim_obs, which holds
the simulated liabilities, the disease status and the current age/age-of-onset
for all family members in each of the n_sim families for the second phenotype.
There is a list containing sim_obs for each phenotype in phen_names.
The last list entry, thresholds, holds the family identifier, the personal
identifier, the role (specified in fam_vec or n_fam) as well as the lower and
upper thresholds for all individuals in all families and all phenotypes.
Note that this tibble has the format required in estimate_liability.
Finally, note that if neither fam_vec nor n_fam are specified, the function
returns the disease status, the current age/age-of-onset, the lower and upper
thresholds, as well as the personal identifier for a single individual, namely
the individual under consideration (called o).
If both fam_vec and n_fam are defined, the user is asked to '
decide on which of the two vectors to use.
Details
This function can be used to simulate the case-control status, the current
age and age-of-onset as well as the lower and upper thresholds for
a variable number of phenotypes for all family members in each of
the n_sim families.
If h2 is a number, simulate_under_LTM simulates the case-
control status, the current age and age-of-onset as well as thresholds
for a single phenotype.
However, if h2 is a numeric vector, if genetic_corrmat and
full_corrmat are two symmetric correlation matrices, and if
phen_names and pop_prev are to numeric vectors holding
the phenotype names and the population prevalences, respectively, then
simulate_under_LTM simulates the case-control status, the current
age and age-of-onset as well as thresholds for two or more (correlated)
phenotypes.
The family members can be specified using one of two possible formats.
Examples
simulate_under_LTM()
#> $sim_obs
#> # A tibble: 1,000 × 26
#> fid g o m f s1 mgm mgf pgm
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 fid_1 0.824 0.375 -1.02 0.00903 1.32 -0.268 2.07 0.172
#> 2 fid_2 0.230 0.495 0.107 0.0394 1.73 0.290 -0.681 -0.772
#> 3 fid_3 -0.115 0.223 0.866 -0.0989 0.0404 0.186 -0.0920 -0.453
#> 4 fid_4 -0.556 0.0704 -1.58 -1.57 0.457 -0.673 0.875 -0.0749
#> 5 fid_5 0.341 0.978 -0.0870 -1.06 -0.877 -0.0642 1.35 -1.22
#> 6 fid_6 0.569 -0.478 0.364 -0.822 -0.801 -1.46 1.22 -0.937
#> 7 fid_7 0.310 -0.772 0.989 0.691 1.02 -0.917 0.979 -0.452
#> 8 fid_8 -0.0480 -0.224 -0.000185 2.00 2.36 0.514 0.323 -0.180
#> 9 fid_9 0.325 0.467 0.380 0.330 1.17 0.279 0.0237 0.250
#> 10 fid_10 0.691 -0.699 0.273 -0.223 -1.66 -0.863 0.857 -0.137
#> # ℹ 990 more rows
#> # ℹ 17 more variables: pgf <dbl>, o_status <lgl>, m_status <lgl>,
#> # f_status <lgl>, s1_status <lgl>, mgm_status <lgl>, mgf_status <lgl>,
#> # pgm_status <lgl>, pgf_status <lgl>, o_aoo <dbl>, m_aoo <dbl>, f_aoo <dbl>,
#> # s1_aoo <dbl>, mgm_aoo <dbl>, mgf_aoo <dbl>, pgm_aoo <dbl>, pgf_aoo <dbl>
#>
#> $thresholds
#> # A tibble: 8,000 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 o -Inf 3.45
#> 2 fid_2 fid_2_1 o -Inf 3.17
#> 3 fid_3 fid_3_1 o -Inf 2.63
#> 4 fid_4 fid_4_1 o -Inf 2.95
#> 5 fid_5 fid_5_1 o -Inf 2.55
#> 6 fid_6 fid_6_1 o -Inf 2.63
#> 7 fid_7 fid_7_1 o -Inf 3.06
#> 8 fid_8 fid_8_1 o -Inf 3.03
#> 9 fid_9 fid_9_1 o -Inf 3.17
#> 10 fid_10 fid_10_1 o -Inf 2.72
#> # ℹ 7,990 more rows
#>
genetic_corrmat <- matrix(0.4, 3, 3)
diag(genetic_corrmat) <- 1
full_corrmat <- matrix(0.6, 3, 3)
diag(full_corrmat) <- 1
simulate_under_LTM(fam_vec = NULL, n_fam = stats::setNames(c(1,1,1,2,2),
c("m","mgm","mgf","s","mhs")))
#> $sim_obs
#> # A tibble: 1,000 × 26
#> fid g o m mgm mgf s1 s2 mhs1 mhs2
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 fid_1 -0.252 0.00373 0.737 0.994 1.08 -0.589 -1.22 -1.20 -0.123
#> 2 fid_2 -0.232 -1.09 0.610 -0.281 -0.115 -0.450 0.107 1.71 0.554
#> 3 fid_3 -0.201 -0.665 -0.270 0.985 -0.607 0.314 -1.76 0.112 0.613
#> 4 fid_4 0.682 0.698 0.547 0.261 -0.225 -0.0136 0.637 0.420 -0.448
#> 5 fid_5 0.147 0.349 -1.99 -0.248 0.587 -0.892 -0.772 -0.694 -2.16
#> 6 fid_6 0.300 0.993 0.267 0.131 1.16 1.50 -0.0664 0.998 -0.856
#> 7 fid_7 -0.835 -0.342 -0.0870 -0.0111 0.117 0.336 0.864 1.21 -0.565
#> 8 fid_8 -0.0510 0.882 1.63 -0.341 1.35 -1.14 -0.444 -0.712 -1.51
#> 9 fid_9 -0.941 0.0908 0.947 1.21 -1.67 -1.08 0.0891 0.721 1.36
#> 10 fid_10 0.642 0.962 -0.127 -0.449 1.45 -0.881 -0.402 -0.567 -0.960
#> # ℹ 990 more rows
#> # ℹ 16 more variables: o_status <lgl>, m_status <lgl>, mgm_status <lgl>,
#> # mgf_status <lgl>, s1_status <lgl>, s2_status <lgl>, mhs1_status <lgl>,
#> # mhs2_status <lgl>, o_aoo <dbl>, m_aoo <dbl>, mgm_aoo <dbl>, mgf_aoo <dbl>,
#> # s1_aoo <dbl>, s2_aoo <dbl>, mhs1_aoo <dbl>, mhs2_aoo <dbl>
#>
#> $thresholds
#> # A tibble: 8,000 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 o -Inf 3.03
#> 2 fid_2 fid_2_1 o -Inf 3.52
#> 3 fid_3 fid_3_1 o -Inf 3.48
#> 4 fid_4 fid_4_1 o -Inf 2.59
#> 5 fid_5 fid_5_1 o -Inf 2.47
#> 6 fid_6 fid_6_1 o -Inf 3.38
#> 7 fid_7 fid_7_1 o -Inf 2.68
#> 8 fid_8 fid_8_1 o -Inf 2.68
#> 9 fid_9 fid_9_1 o -Inf 2.51
#> 10 fid_10 fid_10_1 o -Inf 3.06
#> # ℹ 7,990 more rows
#>
simulate_under_LTM(fam_vec = c("m","f","s1"), n_fam = NULL, add_ind = FALSE,
genetic_corrmat = genetic_corrmat, full_corrmat = full_corrmat, n_sim = 200)
#> $sim_obs
#> # A tibble: 200 × 10
#> fid m f s1 m_status f_status s1_status m_aoo f_aoo s1_aoo
#> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl> <dbl> <dbl> <dbl>
#> 1 fid_1 0.703 1.05 1.03 FALSE FALSE FALSE 61 53 32
#> 2 fid_2 0.597 0.732 -0.700 FALSE FALSE FALSE 46 35 17
#> 3 fid_3 0.174 0.889 0.937 FALSE FALSE FALSE 59 48 29
#> 4 fid_4 1.86 0.328 -0.129 TRUE FALSE FALSE 54 52 33
#> 5 fid_5 -0.470 1.92 2.38 FALSE TRUE TRUE 42 52 41
#> 6 fid_6 -0.0129 0.182 -0.464 FALSE FALSE FALSE 60 58 33
#> 7 fid_7 -0.520 0.519 -0.672 FALSE FALSE FALSE 42 44 18
#> 8 fid_8 -1.05 0.0422 0.169 FALSE FALSE FALSE 51 57 27
#> 9 fid_9 0.638 0.197 1.25 FALSE FALSE FALSE 61 61 38
#> 10 fid_10 0.300 -1.22 -0.404 FALSE FALSE FALSE 60 70 40
#> # ℹ 190 more rows
#>
#> $thresholds
#> # A tibble: 600 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 m -Inf 1.62
#> 2 fid_2 fid_2_1 m -Inf 2.18
#> 3 fid_3 fid_3_1 m -Inf 1.68
#> 4 fid_4 fid_4_1 m 1.85 1.85
#> 5 fid_5 fid_5_1 m -Inf 2.34
#> 6 fid_6 fid_6_1 m -Inf 1.64
#> 7 fid_7 fid_7_1 m -Inf 2.34
#> 8 fid_8 fid_8_1 m -Inf 1.97
#> 9 fid_9 fid_9_1 m -Inf 1.62
#> 10 fid_10 fid_10_1 m -Inf 1.64
#> # ℹ 590 more rows
#>
simulate_under_LTM(fam_vec = c(), n_fam = NULL, add_ind = TRUE, h2 = 0.5,
n_sim = 200, pop_prev = 0.05)
#> Warning: Neither fam_vec nor n_fam is specified...
#> $sim_obs
#> # A tibble: 200 × 5
#> fid g o o_status o_aoo
#> <chr> <dbl> <dbl> <lgl> <dbl>
#> 1 fid_1 -0.143 0.0463 FALSE 10
#> 2 fid_2 -0.219 -1.18 FALSE 33
#> 3 fid_3 -0.826 -1.26 FALSE 30
#> 4 fid_4 -2.25 -2.09 FALSE 26
#> 5 fid_5 -0.678 0.426 FALSE 37
#> 6 fid_6 1.11 1.64 FALSE 30
#> 7 fid_7 0.585 0.360 FALSE 28
#> 8 fid_8 0.286 -0.0754 FALSE 35
#> 9 fid_9 -0.785 -0.875 FALSE 21
#> 10 fid_10 0.254 1.00 FALSE 37
#> # ℹ 190 more rows
#>
#> $thresholds
#> # A tibble: 200 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 o -Inf 3.73
#> 2 fid_2 fid_2_1 o -Inf 2.94
#> 3 fid_3 fid_3_1 o -Inf 3.05
#> 4 fid_4 fid_4_1 o -Inf 3.19
#> 5 fid_5 fid_5_1 o -Inf 2.79
#> 6 fid_6 fid_6_1 o -Inf 3.05
#> 7 fid_7 fid_7_1 o -Inf 3.12
#> 8 fid_8 fid_8_1 o -Inf 2.86
#> 9 fid_9 fid_9_1 o -Inf 3.37
#> 10 fid_10 fid_10_1 o -Inf 2.79
#> # ℹ 190 more rows
#>