magnusdv / forrel

Forensic pedigree analysis and relatedness inference
GNU General Public License v2.0
10 stars 0 forks source link

Rearranging `profileSim` output #12

Closed thoree closed 5 years ago

thoree commented 5 years ago

Below I explore the profileSim function and rearrange the output allowing for transferMarkers and likelihood calculations. Any comments on improved code is appreciated.

library(forrel,quietly = TRUE)
fa = singleton("F")
m = marker(fa, name = "L1")
fa = setMarkers(fa, m)
ch = singleton("C")
m = marker(ch,  name = "L1")
ch = setMarkers(ch, m)
nsim = 5
set.seed(123)
sim = profileSim(list(fa, ch), N = nsim, ids = c("F", "C"), cond = "L1")
simN = vector("list", nsim)
for (i in 1:nsim)
  simN[[i]] = lapply(sim, function(z) z[[i]])
LR(simN, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3]   [,4] [,5]
#> L1 0.25 0.125 0.25 0.0625 0.25
x = nuclearPed(1, father = "F", children = "C")
m = marker(x,  name = "L1")
x = setMarkers(x,m)
for (i in 1:nsim)
  simN[[i]] = transferMarkers(simN[[i]], x, id = c("F","C"))
LR(simN, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3] [,4] [,5]
#> L1 0.25 0.125 0.25    0 0.25

Created on 2018-11-07 by the reprex package (v0.2.1)

magnusdv commented 5 years ago

Ok, nice example. I see that the output of profileSim() was not optimal. I have modified this now, so that you don't have to reorganise manually.

Below is my version of the code, with a couple of simplifications:

library(forrel, quietly = T)

### F and C unrelated singletons
fa = singleton("F")
fa = setMarkers(fa, marker(fa, name = "L1"))
ch = relabel(fa, "C")

# Simulate 5 profiles
simUN = profileSim(list(fa, ch), N = 5, cond = "L1", seed = 123)

# LRs
LR(simUN, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3]   [,4] [,5]
#> L1 0.25 0.125 0.25 0.0625 0.25

### F is the true father
x = nuclearPed(1, father = "F", children = "C")

# Transfer genotypes from each sim
simPO = lapply(simUN, function(s) transferMarkers(s, x))

# LRs
LR(simPO, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3] [,4] [,5]
#> L1 0.25 0.125 0.25    0 0.25
thoree commented 5 years ago

Great - thanks!