Simulate under the liability threshold model (single phenotype).
Source:R/Simulate_under_LTM.R
simulate_under_LTM_single.Rdsimulate_under_LTM_single simulates families and thresholds under
the liability threshold model for a given family structure and a single
phenotype. Please note that it is not possible to simulate different
family structures.
Usage
simulate_under_LTM_single(
fam_vec = c("m", "f", "s1", "mgm", "mgf", "pgm", "pgf"),
n_fam = NULL,
add_ind = TRUE,
h2 = 0.5,
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
A number representing the liability-scale heritability for a single phenotype. Must be non-negative. Note that under the liability threshold model, the heritability must also be at most 1. Defaults to 0.5.
- n_sim
A positive number representing the number of simulations. Defaults to 1000.
- pop_prev
A positive number representing the population prevalence, i.e. the overall prevalence in the population. Must be 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 holding 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.
In addition, 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.
Examples
simulate_under_LTM_single()
#> $sim_obs
#> # A tibble: 1,000 × 26
#> fid g o m f s1 mgm mgf pgm pgf
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 fid_1 -0.695 0.0137 0.371 -0.809 -1.07 -0.144 -1.37 -0.135 1.68
#> 2 fid_2 -1.22 -1.22 -0.374 1.41 0.414 0.00726 -1.46 1.55 -0.810
#> 3 fid_3 1.47 1.52 0.281 1.35 -0.379 -1.10 1.34 -1.87 0.268
#> 4 fid_4 0.630 0.581 -0.132 0.289 0.147 -0.956 -0.214 0.195 1.17
#> 5 fid_5 -0.465 -0.0850 0.344 -2.17 -1.67 1.30 0.199 -0.295 -1.96
#> 6 fid_6 0.249 1.86 0.199 -1.07 1.05 -0.0487 0.527 -1.17 1.29
#> 7 fid_7 0.298 0.994 1.87 -0.0417 1.72 0.778 0.482 0.134 1.75
#> 8 fid_8 1.01 0.350 -0.368 0.739 -0.593 -0.932 -0.479 -0.484 0.720
#> 9 fid_9 0.420 1.27 -0.555 -0.414 -0.262 0.684 0.282 -0.634 -0.668
#> 10 fid_10 0.549 -0.0648 1.95 0.255 0.288 -0.374 0.0944 1.07 0.0511
#> # ℹ 990 more rows
#> # ℹ 16 more variables: 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 2.72
#> 2 fid_2 fid_2_1 o -Inf 3.10
#> 3 fid_3 fid_3_1 o 1.51 1.51
#> 4 fid_4 fid_4_1 o -Inf 2.91
#> 5 fid_5 fid_5_1 o -Inf 3.28
#> 6 fid_6 fid_6_1 o 1.85 1.85
#> 7 fid_7 fid_7_1 o -Inf 3.35
#> 8 fid_8 fid_8_1 o -Inf 3.24
#> 9 fid_9 fid_9_1 o -Inf 3.03
#> 10 fid_10 fid_10_1 o -Inf 2.87
#> # ℹ 7,990 more rows
#>
simulate_under_LTM_single(fam_vec = NULL, n_fam = stats::setNames(c(1,1,1,2),
c("m","mgm","mgf","mhs")))
#> $sim_obs
#> # A tibble: 1,000 × 20
#> fid g o m mgm mgf mhs1 mhs2 o_status m_status
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl> <lgl>
#> 1 fid_1 -0.875 -1.80 -0.455 -0.122 0.264 -1.83 0.571 FALSE FALSE
#> 2 fid_2 -0.0196 -0.343 0.0647 -0.0320 0.580 -0.132 1.47 FALSE FALSE
#> 3 fid_3 0.218 0.197 0.0768 -0.938 -0.0490 1.93 0.822 FALSE FALSE
#> 4 fid_4 0.202 1.18 0.0658 -0.786 -0.0139 0.634 0.924 FALSE FALSE
#> 5 fid_5 0.0784 -0.195 -1.82 -0.158 0.419 -0.596 -1.22 FALSE FALSE
#> 6 fid_6 0.130 -0.241 0.447 0.356 0.958 -2.05 0.130 FALSE FALSE
#> 7 fid_7 -0.0486 1.02 0.966 0.336 -0.514 -0.407 0.0377 FALSE FALSE
#> 8 fid_8 -0.943 -0.501 0.772 1.93 0.720 0.910 0.452 FALSE FALSE
#> 9 fid_9 -0.488 -1.18 0.693 -0.385 -0.0115 1.08 0.0553 FALSE FALSE
#> 10 fid_… -0.801 -0.697 -0.566 0.572 -1.21 -1.32 1.13 FALSE FALSE
#> # ℹ 990 more rows
#> # ℹ 10 more variables: mgm_status <lgl>, mgf_status <lgl>, mhs1_status <lgl>,
#> # mhs2_status <lgl>, o_aoo <dbl>, m_aoo <dbl>, mgm_aoo <dbl>, mgf_aoo <dbl>,
#> # mhs1_aoo <dbl>, mhs2_aoo <dbl>
#>
#> $thresholds
#> # A tibble: 6,000 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 o -Inf 3.31
#> 2 fid_2 fid_2_1 o -Inf 3.52
#> 3 fid_3 fid_3_1 o -Inf 2.76
#> 4 fid_4 fid_4_1 o -Inf 2.68
#> 5 fid_5 fid_5_1 o -Inf 3.24
#> 6 fid_6 fid_6_1 o -Inf 3.52
#> 7 fid_7 fid_7_1 o -Inf 2.87
#> 8 fid_8 fid_8_1 o -Inf 3.06
#> 9 fid_9 fid_9_1 o -Inf 2.76
#> 10 fid_10 fid_10_1 o -Inf 2.68
#> # ℹ 5,990 more rows
#>
simulate_under_LTM_single(fam_vec = c("m","f","s1"), n_fam = NULL, add_ind = FALSE,
h2 = 0.5, n_sim = 500, pop_prev = .05)
#> $sim_obs
#> # A tibble: 500 × 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.457 -0.146 0.190 FALSE FALSE FALSE 43 40 17
#> 2 fid_2 -1.12 -0.109 -1.28 FALSE FALSE FALSE 45 51 25
#> 3 fid_3 1.70 1.22 0.704 TRUE FALSE FALSE 77 38 12
#> 4 fid_4 -1.45 0.957 -0.819 FALSE FALSE FALSE 52 43 25
#> 5 fid_5 -1.74 -2.18 -1.32 FALSE FALSE FALSE 65 56 35
#> 6 fid_6 -0.410 -0.555 -0.774 FALSE FALSE FALSE 51 39 21
#> 7 fid_7 0.0720 -0.573 0.236 FALSE FALSE FALSE 58 60 32
#> 8 fid_8 -1.03 0.0949 -0.621 FALSE FALSE FALSE 41 38 11
#> 9 fid_9 0.517 -0.766 1.44 FALSE FALSE FALSE 44 39 15
#> 10 fid_10 -0.190 -0.337 -0.423 FALSE FALSE FALSE 51 49 22
#> # ℹ 490 more rows
#>
#> $thresholds
#> # A tibble: 1,500 × 5
#> fid indiv_ID role lower upper
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 fid_1 fid_1_1 m -Inf 2.55
#> 2 fid_2 fid_2_1 m -Inf 2.48
#> 3 fid_3 fid_3_1 m 1.70 1.70
#> 4 fid_4 fid_4_1 m -Inf 2.21
#> 5 fid_5 fid_5_1 m -Inf 1.84
#> 6 fid_6 fid_6_1 m -Inf 2.25
#> 7 fid_7 fid_7_1 m -Inf 2.02
#> 8 fid_8 fid_8_1 m -Inf 2.63
#> 9 fid_9 fid_9_1 m -Inf 2.51
#> 10 fid_10 fid_10_1 m -Inf 2.25
#> # ℹ 1,490 more rows
#>
simulate_under_LTM_single(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.357 -1.41 FALSE 40
#> 2 fid_2 1.36 1.57 FALSE 32
#> 3 fid_3 -1.04 -1.14 FALSE 22
#> 4 fid_4 -0.367 -0.150 FALSE 35
#> 5 fid_5 1.20 1.13 FALSE 15
#> 6 fid_6 -0.814 0.189 FALSE 11
#> 7 fid_7 1.06 1.99 TRUE 59
#> 8 fid_8 0.806 0.0122 FALSE 17
#> 9 fid_9 -0.298 -0.0825 FALSE 38
#> 10 fid_10 0.656 -0.114 FALSE 30
#> # ℹ 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 2.67
#> 2 fid_2 fid_2_1 o -Inf 2.97
#> 3 fid_3 fid_3_1 o -Inf 3.33
#> 4 fid_4 fid_4_1 o -Inf 2.86
#> 5 fid_5 fid_5_1 o -Inf 3.57
#> 6 fid_6 fid_6_1 o -Inf 3.70
#> 7 fid_7 fid_7_1 o 1.99 1.99
#> 8 fid_8 fid_8_1 o -Inf 3.50
#> 9 fid_9 fid_9_1 o -Inf 2.75
#> 10 fid_10 fid_10_1 o -Inf 3.05
#> # ℹ 190 more rows
#>