founder_inbreeding and likelihood

[thoree]

The below example shows a well known fact: inbred founders is not yet accomodated in likelihood ?! I would just like that too mention that this feature would be important also for LR calculations in forensics and Hilde KB's project

library(pedtools)
library(pedprobr)
H = nuclearPed(father ="FA", mother = "MO", child  ="CH")
m = marker(H, FA = 1:2)
l1 = likelihood(H, m)
founder_inbreeding(H, "FA") = 0.25
l2 = likelihood(H, m)
l1 == l2
#> [1] TRUE
single = singleton("FA")
m = marker(single, FA = 1:2)
founder_inbreeding(single, "FA") = 0.25
l3 = likelihood(single, m)
l1 == l3
#> [1] TRUE

Created on 2018-08-15 by the reprex package (v0.2.0).

[magnusdv]

Yes, I plan to implement this shortly. Thanks for simple examples - I may reuse them in unit tests.

[magnusdv]

A topic for discussion:
Situations where the concept of inbred founders leads to problematic behaviour of the likelihood.

Terminology: If P is a pedigree with inbred founders, let us say that a realisation of P is any pedigree which i) extends P (i.e. contains P as a subgraph), ii) has outbred founders, and iii) preserves ("realises") all inbreeding coefficients in P.

For the likelihood to be well-defined, the following requirement seems natural:

• For any realisation P2 of P, the likelihood (given the same marker data) should be the same.

This is clearly satisfied in many typical situations, but here is one where it's not: If the likelihood is computed under an unbalanced mutation model. Different choices of P2 will then in general have different likelihoods.

These cases should probably issue a warning message, at least.
Any thoughts, @thoree?

[thoree]

Good point. I agree, consistency in the sense "For any realisation P2 of P, the likelihood (given the same marker data) should be the same" cannot generally be expected for unbalanced mutation models. A warning to this effect seems in order.

[magnusdv]

I'm happy to report that a first version of this is now implemented! See 4639bcf for details.
Only autosomal markers for now, and not linked markers.

Here is an example validating the procedure in a very simple case:

library(pedtools)
library(pedprobr)

# Create inbred individual
x = addChildren(nuclearPed(2, sex=1:2), 3, 4, 1)

# Add a marker
p = 0.1; q = 1-p
m = marker(x, '5'=1, alleles=1:2, afreq=c(p,q))
x = setMarkers(x, m)

# Extract singleton and set founder inbreeding
y = subset(x, 5)
founder_inbreeding(y, 5) = 0.25

# These should now all be equal:
likelihood(x, 1, verbose = F)
#> [1] 0.0325

likelihood(y, 1)
#> [1] 0.0325

0.25 * p + 0.75 * p^2 # theory
#> [1] 0.0325

[magnusdv]

Back to the issue with muations in pedigrees with inbred founder. Here is another problem with this:
Consider an inbred singleton, for which a heterozygous genotype is observed. What should the likelihood be? I don't this has a well-behaved answer.

[thoree]

I see your interesting point Magnus @magnusdv. Perhaps, at least for the time being, mutations should not be allowed whenever inbred founders are specified. Below is an example to check my understanding (it's long but I enjoy learning to code), where the Jacquard coefficients of parents coincide but the likelihood
of the their heterozygous child differ with a balanced mutation model:

library(pedtools)
library(pedprobr)
library(forrel)
library(ribd)
x1 = halfCousinsPed(0, child =TRUE)
p = c(0.4,0.6); a = 1:2; R = 0.1 ; phi = 0.25
M = mutationMatrix(a, "proportional", p, rate = R)
m1 = marker(x1,"6" = c(1,2), afreq = p, alleles = a, mutmat = M)
x2 = cousinsPed(0, removal = 1)
x2 = addChildren(x2, father =3, mother = 6)
x2 = relabel(x2, c(1:3,5,4,7,6))
m2 = marker(x2,"6" = c(1,2), afreq = p, alleles = a, mutmat = M)
# plot(x1, m1, skip.empty.genotypes = TRUE)
# plot(x2, m2, skip.empty.genotypes = TRUE)
founder_inbreeding(x1,1) = phi
founder_inbreeding(x1,2) = phi
founder_inbreeding(x1,4) = phi
founder_inbreeding(x2,1) = phi
founder_inbreeding(x2,2) = phi
founder_inbreeding(x2,4) = phi
l1 = likelihood(x1,m1)
#> Tip: To optimize speed, consider breaking loops before calling 'likelihood'. See ?breakLoops.
#> Loop breakers: 3
l2 = likelihood(x2,m2)
#> Tip: To optimize speed, consider breaking loops before calling 'likelihood'. See ?breakLoops.
#> Loop breakers: 3
l1 -l2
#> [1] -0.006137471
all(ibd_identity(x1, parents(x1,6)) == 
      ibd_identity(x2, parents(x2,6)))
#> Function calls:
#>   phi2  = 6
#>   phi3  = 8 (recurse: 3)
#>   phi4  = 2
#>   phi22 = 5 (recurse: 2)
#> Time used: 0.0 seconds
#> Comparison with `identity` package:
#> skipped. (Pedigree has inbred founders.)
#> Function calls:
#>   phi2  = 7
#>   phi3  = 11 (recurse: 5)
#>   phi4  = 2
#>   phi22 = 5 (recurse: 2)
#> Time used: 0.0 seconds
#> Comparison with `identity` package:
#> skipped. (Pedigree has inbred founders.)
#> [1] TRUE

Created on 2018-08-17 by the reprex package (v0.2.0).

[magnusdv]

Closing this, as the code examples here are getting a bit dated.

The quirks of likelihoods in pedigrees with founder inbreeding remain interesting, of course, so new issues on this are welcome.