Closed thoree closed 5 years ago
ibd_kappa suprises me below, not ibd_identity:
ibd_kappa
ibd_identity
library(ribd, quietly = TRUE) x = nuclearPed(2, father = "FA", mother = "MO", children = c("B1", "B2")) f1 = 0.1 f2 = 0.2 f = c(f1, f2) founder_inbreeding(x, c("FA", "MO")) = f ibd_kappa(x, ids = c("B1", "B2")) #> Inbred individuals: B1, B2 #> The kappa coefficients are undefined for inbred individuals. #> [1] NA NA NA f00 = (1-f1)*(1-f2) f10 = f1*(1-f2) f01 = (1-f1)*f2 f11 = f1*f2 k0 = f00*0.25 k1 = f00*0.5+f10*0.5+f01*0.5 k2 = f00*0.25+f10*0.5+f01*0.5+f11 k = c(k0, k1, k2) j = ibd_identity(x, ids= c("B1", "B2"))[9:7] k-j #> [1] -1.387779e-16 2.775558e-16 0.000000e+00
Created on 2019-01-22 by the reprex package (v0.2.1)
Thanks. This was caused by rounding errors; the current version should be more robust.
ibd_kappa
suprises me below, notibd_identity
:Created on 2019-01-22 by the reprex package (v0.2.1)