Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rearranging profileSim output #12

Closed
thoree opened this issue Nov 7, 2018 · 2 comments
Closed

Rearranging profileSim output #12

thoree opened this issue Nov 7, 2018 · 2 comments

Comments

@thoree
Copy link
Contributor

thoree commented Nov 7, 2018

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
Copy link
Owner

magnusdv commented Nov 8, 2018

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:

  • The seed can be set in the function call
  • By default all indivis are included (so the ids parameter is not necessary in this case)
  • Since the only difference between ch and fa is the label, you can make ch in one line using relabel().
  • The transferMarkers() part is a one-liner!
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

@magnusdv magnusdv closed this as completed Nov 8, 2018
@thoree
Copy link
Contributor Author

thoree commented Nov 8, 2018

Great - thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants