Defines a log normal model for peak height variability

Usage

log_normal_model(
  template,
  degradation = rep(0, length(template)),
  LSAE = stats::setNames(rep(1, length(model_settings$locus_names)),
    model_settings$locus_names),
  c2,
  k2,
  model_settings
)

Arguments

template

Numeric vector

degradation

Numeric vector of same length as template. Degradation parameters for each contributor.

LSAE

Numeric vector (named) with Locus Specific Amplification Efficiencies. See sample_LSAE. Defaults to 1 for each locus.

c2

Numeric. Allele variance parameter.

k2

Optionally a numeric vector with stutter variance parameters. See sample_log_normal_stutter_variance.

model_settings

List. Possible parameters:

locus_names. Character vector.
degradation_parameter_cap. Numeric.
c2_prior. Numeric of length two with shape and scale.
LSAE_variance_prior. Numeric of length one.
detection_threshold. Numeric vector (named) with Detection Thresholds. Defaults to 50 for each locus.
size_regression. Function, see read_size_regression.
stutter_model. Optionally a stutter_model object that gives expected stutter heights. See global_stutter_model.
stutter_variability. Optionally peak height variability parameters for stutters. Required when stutter_model is supplied.

Value

Object of class pg_model.

Details

Define a log normal model for peak height variability with the parametrisation as described by Bright et al. The model may then be used to sample DNA profiles using the sample_mixture_from_genotypes function. Alternatively, to sample many models and profiles in one go with parameters according to a specified distribution, the sample_mixtures function can be used.

References

Bright, J.A. et al. (2016). Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles. Forensic Science International: Genetics, 23, 226-239. doi:10.1016/j.fsigen.2016.05.007

Examples

gf <- gf_configuration()
freqs <- read_allele_freqs(system.file("extdata","FBI_extended_Cauc_022024.csv",
                           package = "simDNAmixtures"))

k2 <- sample_log_normal_stutter_variance(gf$log_normal_settings$stutter_variability)

model <- log_normal_model(template = 1e3, c2 = 15, k2 = k2,
                          model_settings = gf$log_normal_settings)
model
#> $locus_names
#>  [1] "D3S1358"  "vWA"      "D16S539"  "CSF1PO"   "TPOX"     "AMEL"    
#>  [7] "D8S1179"  "D21S11"   "D18S51"   "D2S441"   "D19S433"  "TH01"    
#> [13] "FGA"      "D22S1045" "D5S818"   "D13S317"  "D7S820"   "SE33"    
#> [19] "D10S1248" "D1S1656"  "D12S391"  "D2S1338" 
#> 
#> $detection_threshold
#>  D3S1358      vWA  D16S539   CSF1PO     TPOX     AMEL  D8S1179   D21S11 
#>       75       75       75       75       75       75      100      100 
#>   D18S51   D2S441  D19S433     TH01      FGA D22S1045   D5S818  D13S317 
#>      100       60       60       60       60       80       80       80 
#>   D7S820     SE33 D10S1248  D1S1656  D12S391  D2S1338 
#>       80       80      100      100      100      100 
#> 
#> $parameters
#> $parameters$model
#> [1] "log_normal_model"
#> 
#> $parameters$template
#> [1] 1000
#> 
#> $parameters$degradation
#> [1] 0
#> 
#> $parameters$c2
#> [1] 15
#> 
#> $parameters$k2
#>       k2BackStutter    k2ForwardStutter    k22bpBackStutter k2DoubleBackStutter 
#>            8.904872            6.825219            1.497035            5.542487 
#> 
#> $parameters$LSAE
#>  D3S1358      vWA  D16S539   CSF1PO     TPOX     AMEL  D8S1179   D21S11 
#>        1        1        1        1        1        1        1        1 
#>   D18S51   D2S441  D19S433     TH01      FGA D22S1045   D5S818  D13S317 
#>        1        1        1        1        1        1        1        1 
#>   D7S820     SE33 D10S1248  D1S1656  D12S391  D2S1338 
#>        1        1        1        1        1        1 
#> 
#> 
#> $size_regression
#> function (locus, allele) 
#> {
#>     if (has_exceptions) {
#>         locus_exceptions <- exceptions[[locus]]
#>         if (!is.null(locus_exceptions)) {
#>             if (allele %in% names(locus_exceptions)) {
#>                 size <- locus_exceptions[[allele]]
#>                 return(size)
#>             }
#>         }
#>     }
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No size regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     allele_numeric <- as.numeric(allele)
#>     if (has_repeat_length_by_marker) {
#>         repeat_length <- repeat_length_by_marker[[locus]]
#>         if (!is.null(repeat_length)) {
#>             allele_numeric <- repeats_to_decimals(allele_numeric, 
#>                 repeat_length)
#>         }
#>     }
#>     intercept + slope * allele_numeric
#> }
#> <bytecode: 0x55ce0eb8eb40>
#> <environment: 0x55ce113d0360>
#> 
#> $build_expected_profile_and_sample_peak_heights
#> function (genotypes) 
#> {
#>     expected_profile <- log_normal_model_build_expected_profile(model, 
#>         genotypes)
#>     x <- log_normal_model_sample_peak_heights(model, expected_profile, 
#>         model$stutter_variability)
#>     x$LSAE <- LSAE[x$Marker]
#>     x
#> }
#> <bytecode: 0x55ce0ffca2b0>
#> <environment: 0x55ce08e2e830>
#> 
#> $sample_name_suffix
#> [1] "N1_1000_a"
#> 
#> $stutter_model
#> $stutter_types
#> $stutter_types$BackStutter
#> $name
#> [1] "BackStutter"
#> 
#> $delta
#> [1] -1
#> 
#> $applies_to_all_loci
#> [1] TRUE
#> 
#> $regression
#> function (locus, allele) 
#> {
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No stutter regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     max(min_stutter_ratio, intercept + slope * as.numeric(allele))
#> }
#> <bytecode: 0x55ce0edc0388>
#> <environment: 0x55ce098178d8>
#> 
#> $exceptions
#> $exceptions$D8S1179
#>       8      10      11      12      13      14      15      16      18 
#> 0.03152 0.04864 0.05584 0.06524 0.05913 0.06564 0.07019 0.07451 0.08698 
#> 
#> $exceptions$D21S11
#>      27      28      29      30    30.2      31    31.2      32    32.2    33.2 
#> 0.04241 0.05957 0.05909 0.07504 0.05746 0.07995 0.06253 0.08734 0.07836 0.08319 
#> 
#> $exceptions$D2S441
#>      10      11      12      13      14      15      16 
#> 0.03342 0.05468 0.06129 0.04058 0.03950 0.04500 0.05982 
#> 
#> $exceptions$D19S433
#>     6.2       8       9      10      11    11.1    11.2      12    12.1    12.2 
#> 0.00005 0.00876 0.00876 0.02618 0.03489 0.02618 0.04360 0.04360 0.00005 0.05231 
#>      13    13.2      14    14.2      15    15.2      16    16.2      17    17.2 
#> 0.05231 0.06102 0.06102 0.06973 0.06973 0.07844 0.07844 0.08715 0.08715 0.09151 
#>      18    18.2    19.2 
#> 0.09586 0.10457 0.11328 
#> 
#> $exceptions$TH01
#>       4       5       6       7       8     8.3       9     9.3      10    10.3 
#> 0.00223 0.00799 0.01375 0.01951 0.02527 0.00799 0.03103 0.01375 0.03679 0.01375 
#>      11      12    13.3 
#> 0.04255 0.04831 0.02527 
#> 
#> $exceptions$SE33
#>      13      14      15      16      17      18      19      20      21    22.2 
#> 0.05466 0.05775 0.06570 0.07151 0.07849 0.08581 0.08929 0.09947 0.10593 0.06981 
#>    23.2    24.2    25.2    26.2    27.2    28.2    29.2    30.2    31.2    32.2 
#> 0.06975 0.08019 0.08294 0.08884 0.09311 0.09721 0.10606 0.11024 0.11982 0.11991 
#> 
#> $exceptions$D1S1656
#>       8       9      10      11      12      13    13.3      14    14.3      15 
#> 0.02644 0.03559 0.04017 0.05389 0.05847 0.06304 0.05389 0.07677 0.03559 0.08134 
#>    15.3      16    16.3      17    17.3      18    18.3    19.3 
#> 0.04932 0.09351 0.05389 0.09964 0.06304 0.10879 0.07219 0.08134 
#> 
#> $exceptions$D2S1338
#>      16      17      18      19      20      21      22      23      24      25 
#> 0.04892 0.05731 0.06614 0.06862 0.07423 0.08291 0.06979 0.07556 0.08495 0.09530 
#>      26 
#> 0.10131 
#> 
#> 
#> $get_expected_stutter_ratio
#> function (locus, allele) 
#> {
#>     exception <- NULL
#>     if (!is.null(stutter$exceptions[[locus]])) {
#>         if (as.character(allele) %in% names(stutter$exceptions[[locus]])) {
#>             exception <- stutter$exceptions[[locus]][[as.character(allele)]]
#>         }
#>     }
#>     if (isTRUE(exception > 0)) {
#>         return(exception)
#>     }
#>     else {
#>         return(stutter$regression(locus, allele))
#>     }
#> }
#> <bytecode: 0x55ce0ef7cbd0>
#> <environment: 0x55ce096f0860>
#> 
#> attr(,"class")
#> [1] "stutter_type"
#> 
#> $stutter_types$ForwardStutter
#> $name
#> [1] "ForwardStutter"
#> 
#> $delta
#> [1] 1
#> 
#> $applies_to_all_loci
#> [1] TRUE
#> 
#> $regression
#> function (locus, allele) 
#> {
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No stutter regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     max(min_stutter_ratio, intercept + slope * as.numeric(allele))
#> }
#> <bytecode: 0x55ce0edc0388>
#> <environment: 0x55ce096dbb78>
#> 
#> $get_expected_stutter_ratio
#> function (locus, allele) 
#> {
#>     exception <- NULL
#>     if (!is.null(stutter$exceptions[[locus]])) {
#>         if (as.character(allele) %in% names(stutter$exceptions[[locus]])) {
#>             exception <- stutter$exceptions[[locus]][[as.character(allele)]]
#>         }
#>     }
#>     if (isTRUE(exception > 0)) {
#>         return(exception)
#>     }
#>     else {
#>         return(stutter$regression(locus, allele))
#>     }
#> }
#> <bytecode: 0x55ce0ef7cbd0>
#> <environment: 0x55ce095c7978>
#> 
#> attr(,"class")
#> [1] "stutter_type"
#> 
#> $stutter_types$`2bpBackStutter`
#> $name
#> [1] "2bpBackStutter"
#> 
#> $delta
#> [1]  0 -2
#> 
#> $applies_to_all_loci
#> [1] FALSE
#> 
#> $applies_to_loci
#> [1] "SE33"    "D1S1656"
#> 
#> $repeat_length_by_marker
#>  D3S1358      vWA  D16S539   CSF1PO     TPOX  D8S1179   D21S11   D18S51 
#>        4        4        4        4        4        4        4        4 
#>   D2S441  D19S433     TH01      FGA D22S1045   D5S818  D13S317   D7S820 
#>        4        4        4        4        3        4        4        4 
#>     SE33 D10S1248  D1S1656  D12S391  D2S1338 
#>        4        4        4        4        4 
#> 
#> $regression
#> function (locus, allele) 
#> {
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No stutter regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     max(min_stutter_ratio, intercept + slope * as.numeric(allele))
#> }
#> <bytecode: 0x55ce0edc0388>
#> <environment: 0x55ce095a7d30>
#> 
#> $get_expected_stutter_ratio
#> function (locus, allele) 
#> {
#>     exception <- NULL
#>     if (!is.null(stutter$exceptions[[locus]])) {
#>         if (as.character(allele) %in% names(stutter$exceptions[[locus]])) {
#>             exception <- stutter$exceptions[[locus]][[as.character(allele)]]
#>         }
#>     }
#>     if (isTRUE(exception > 0)) {
#>         return(exception)
#>     }
#>     else {
#>         return(stutter$regression(locus, allele))
#>     }
#> }
#> <bytecode: 0x55ce0ef7cbd0>
#> <environment: 0x55ce094d00b8>
#> 
#> attr(,"class")
#> [1] "stutter_type"
#> 
#> $stutter_types$DoubleBackStutter
#> $name
#> [1] "DoubleBackStutter"
#> 
#> $delta
#> [1] -2
#> 
#> $applies_to_all_loci
#> [1] TRUE
#> 
#> $regression
#> function (locus, allele) 
#> {
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No stutter regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     max(min_stutter_ratio, intercept + slope * as.numeric(allele))
#> }
#> <bytecode: 0x55ce0edc0388>
#> <environment: 0x55ce094b47f0>
#> 
#> $get_expected_stutter_ratio
#> function (locus, allele) 
#> {
#>     exception <- NULL
#>     if (!is.null(stutter$exceptions[[locus]])) {
#>         if (as.character(allele) %in% names(stutter$exceptions[[locus]])) {
#>             exception <- stutter$exceptions[[locus]][[as.character(allele)]]
#>         }
#>     }
#>     if (isTRUE(exception > 0)) {
#>         return(exception)
#>     }
#>     else {
#>         return(stutter$regression(locus, allele))
#>     }
#> }
#> <bytecode: 0x55ce0ef7cbd0>
#> <environment: 0x55ce0937d738>
#> 
#> attr(,"class")
#> [1] "stutter_type"
#> 
#> 
#> $size_regression
#> function (locus, allele) 
#> {
#>     if (has_exceptions) {
#>         locus_exceptions <- exceptions[[locus]]
#>         if (!is.null(locus_exceptions)) {
#>             if (allele %in% names(locus_exceptions)) {
#>                 size <- locus_exceptions[[allele]]
#>                 return(size)
#>             }
#>         }
#>     }
#>     regression_locus <- regression_df_by_locus[[locus]]
#>     if (is.null(regression_locus)) {
#>         stop("No size regression available for locus ", locus)
#>     }
#>     intercept <- regression_locus$Intercept
#>     slope <- regression_locus$Slope
#>     allele_numeric <- as.numeric(allele)
#>     if (has_repeat_length_by_marker) {
#>         repeat_length <- repeat_length_by_marker[[locus]]
#>         if (!is.null(repeat_length)) {
#>             allele_numeric <- repeats_to_decimals(allele_numeric, 
#>                 repeat_length)
#>         }
#>     }
#>     intercept + slope * allele_numeric
#> }
#> <bytecode: 0x55ce0eb8eb40>
#> <environment: 0x55ce113d0360>
#> 
#> $sex_locus_name
#> [1] "AMEL"
#> 
#> $add_expected_stutter
#> function (...) 
#> allele_specific_stutter_model_add_expected_stutter(stutter_model, 
#>     ...)
#> <bytecode: 0x55ce0f5f01c0>
#> <environment: 0x55ce09375530>
#> 
#> attr(,"class")
#> [1] "stutter_model"
#> 
#> $stutter_variability
#> $stutter_variability$BackStutter
#> $stutter_variability$BackStutter$k2_prior
#> [1] 1.884 7.686
#> 
#> $stutter_variability$BackStutter$inversely_proportional_to_parent
#> [1] TRUE
#> 
#> $stutter_variability$BackStutter$max_stutter_ratio
#> [1] 0.3
#> 
#> 
#> $stutter_variability$ForwardStutter
#> $stutter_variability$ForwardStutter$k2_prior
#> [1] 2.144 4.507
#> 
#> $stutter_variability$ForwardStutter$inversely_proportional_to_parent
#> [1] FALSE
#> 
#> $stutter_variability$ForwardStutter$max_stutter_ratio
#> [1] 0.15
#> 
#> 
#> $stutter_variability$`2bpBackStutter`
#> $stutter_variability$`2bpBackStutter`$k2_prior
#> [1] 2.189 1.431
#> 
#> $stutter_variability$`2bpBackStutter`$inversely_proportional_to_parent
#> [1] FALSE
#> 
#> $stutter_variability$`2bpBackStutter`$max_stutter_ratio
#> [1] 0.1
#> 
#> 
#> $stutter_variability$DoubleBackStutter
#> $stutter_variability$DoubleBackStutter$k2_prior
#> [1] 3.429 2.032
#> 
#> $stutter_variability$DoubleBackStutter$inversely_proportional_to_parent
#> [1] FALSE
#> 
#> $stutter_variability$DoubleBackStutter$max_stutter_ratio
#> [1] 0.05
#> 
#> 
#> 
#> attr(,"class")
#> [1] "pg_model"