Sample cases under Hp for sample-dependent error probabilities, $w_t$ and $w_r$

Same latent genotype, Z, with independent errors for true donor (D) and suspect (S).

Usage

sample_data_Hp_wTwR(n, wT, wR, p, ...)

Arguments

n

number of samples

wT

error probability for donor sample

wR

error probability for PoI sample

p

list of genotype probabilities (length is number of loci) or vector of length 3 for single locus

...

Passed on to add_errors_to_genotypes()

Value

list of two matrices, each of size n x loci with genotype in 0/1/2 format resembling the situation in real life.

Examples

sample_data_Hp_wTwR(n = 10, wT = 0.3, wR = 1e-6, p = c(0.25, 0.25, 0.5))
#> $xT
#>       [,1]
#>  [1,]    0
#>  [2,]    1
#>  [3,]    1
#>  [4,]    1
#>  [5,]    1
#>  [6,]    0
#>  [7,]    1
#>  [8,]    2
#>  [9,]    1
#> [10,]    1
#> 
#> $xR
#>       [,1]
#>  [1,]    2
#>  [2,]    1
#>  [3,]    2
#>  [4,]    1
#>  [5,]    2
#>  [6,]    0
#>  [7,]    1
#>  [8,]    2
#>  [9,]    0
#> [10,]    2
#> 
sample_data_Hp_wTwR(n = 10, wT = 0.3, wR = 1e-6, p = list(
  c(0.25, 0.25, 0.5), c(0.1, 0.8, 0.1)))
#> $xT
#>       [,1] [,2]
#>  [1,]    1    0
#>  [2,]    2    1
#>  [3,]    1    1
#>  [4,]    0    2
#>  [5,]    1    1
#>  [6,]    0    1
#>  [7,]    2    2
#>  [8,]    2    2
#>  [9,]    0    0
#> [10,]    2    1
#> 
#> $xR
#>       [,1] [,2]
#>  [1,]    1    0
#>  [2,]    2    1
#>  [3,]    1    1
#>  [4,]    0    1
#>  [5,]    2    2
#>  [6,]    0    1
#>  [7,]    2    1
#>  [8,]    2    1
#>  [9,]    0    1
#> [10,]    2    1
#> 
  
cases <- sample_data_Hp_wTwR(n = 1000, wT = 0, wR = 0, p = c(0.25, 0.25, 0.5))
table(xT = cases$xT, xR = cases$xR)
#>    xR
#> xT    0   1   2
#>   0 253   0   0
#>   1   0 239   0
#>   2   0   0 508
cases <- sample_data_Hp_wTwR(n = 1000, wT = 0.1, wR = 0, p = c(0.25, 0.25, 0.5))
table(xT = cases$xT, xR = cases$xR)
#>    xR
#> xT    0   1   2
#>   0 178  31   2
#>   1  41 218  93
#>   2   3  20 414
cases <- sample_data_Hp_wTwR(n = 1000, wT = 0, wR = 0.1, p = c(0.25, 0.25, 0.5))
table(xT = cases$xT, xR = cases$xR)
#>    xR
#> xT    0   1   2
#>   0 197  42   2
#>   1  18 229  24
#>   2   4  95 389
  
cases <- sample_data_Hp_wTwR(n = 1000, wT = 1e-1, wR = 1e-8, p = c(0.25, 0.25, 0.5))
tab <- table(xT = cases$xT, xR = cases$xR)
tab
#>    xR
#> xT    0   1   2
#>   0 203  19   3
#>   1  54 202 100
#>   2   3  25 391
estimate_w(tab)
#> [1] 0.05907991

cases <- sample_data_Hp_wTwR(n = 1000, wT = 0, wR = 0, p = c(0.1, 0.7, 0.2))
tab <- table(xT = cases$xT, xR = cases$xR)
diag(tab/sum(tab))
#>     0     1     2 
#> 0.100 0.686 0.214