Nothing
test_that("likelihood() catches input errors", {
x = singleton(1)
y = nuclearPed(1)
m = marker(x)
expect_error(likelihood(x, y), "Invalid input for argument `markers`")
expect_error(likelihood(x, y, m), "Invalid input for argument `markers`")
expect_error(likelihood(x, y, marker = m), "Invalid input for argument `peelOrder`")
expect_error(likelihood(x, marker = list(1)), "Invalid input for argument `markers`.")
expect_error(likelihood(x, marker = NA), "`markers` is a logical of length 1")
expect_error(likelihood(x, marker = 1), "Marker index out of range: 1")
expect_error(likelihood(x, marker = "a"), "Unknown marker name: a")
peds = list(x,y)
expect_error(likelihood(peds, m), fixed = T,
"`likelihood.list()` requires `markers` to be a vector of marker names or indices.")
expect_error(likelihood(peds, list()), fixed = T,
"`likelihood.list()` requires `markers` to be a vector of marker names or indices.")
expect_error(likelihood(peds, 1), "Marker index out of range: 1")
peds2 = list(setMarkers(x,marker(x)), y)
expect_error(likelihood(peds2), "The pedigree components have different number of markers attached")
expect_error(likelihood(peds2,1), "Marker index out of range: 1")
})
test_that("empty markers give likelihood 1", {
x = nuclearPed(1)
expect_equal(likelihood(x, marker(x)), 1)
expect_equal(likelihood(x, marker(x, alleles=1:2)), 1)
expect_equal(likelihood(x, marker(x, alleles=letters[1:3])), 1)
x2 = cousinPed(1)
expect_equal(likelihood(x2, marker(x2)), 1)
x3 = halfCousinPed(0)
expect_equal(likelihood(x3, marker(x3)), 1)
x4 = cousinPed(0, child=T)
expect_equal(likelihood(x4, marker(x4)), 1)
})
test_that("nuclear likelihoods are correct, HW-like", {
x = nuclearPed(1)
p = 0.1; q = 1-p
m1 = marker(x, '1'=1, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m1), p^2)
m2 = marker(x, '1'=1:2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m2), 2*p*q)
m3 = marker(x, '3'=1, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m3), p^2)
m4 = marker(x, '3'=2:1, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m4), 2*p*q)
m5 = marker(x, '1'=1, '2'=2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m5), p^2*q^2)
m6 = marker(x, '1'=1, '2'=2, '3'=1:2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m6), p^2*q^2)
m7 = marker(x, '1'=1, '2'=1, '3'=2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m7), 0)
m8 = marker(x, '1'=1, '3'=2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m8), 0)
})
test_that("nuclear likelihoods are correct, multiallelic markers", {
x = nuclearPed(1)
p = 0.1; q = 1-p
m1 = marker(x, alleles=1:3, afreq=c(p,q,0))
expect_equal(likelihood(x, m1), 1)
m2 = marker(x, '1'=1, alleles=1:3, afreq=c(p,q,0))
expect_equal(likelihood(x, m2), p^2)
m3 = marker(x, '1'=2, alleles=1:3, afreq=c(p,p,q-p))
expect_equal(likelihood(x, m3), p^2)
#m4 = marker(x, '1'=2, alleles=1:3, afreq=c(p,p,q-p))
#expect_equal(likelihood(x, m4), p^2)
m5 = marker(x, '1'=1, '2'=10, alleles=1:10)
expect_equal(likelihood(x, m5), .1^4)
})
test_that("likelihoods are correct in trio SNP examples", {
x = nuclearPed(1)
p = 0.2; q = 1 - p
# All heterozygous
m1 = marker(x, '1'=1:2, '2'=1:2, '3'=1:2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m1), 2*p^2*q^2)
# Same, without mom
m2 = marker(x, '1'=1:2, '3'=1:2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m2), p*q)
})
test_that("likelihoods are correct in sibling SNP example", {
x = nuclearPed(2)
p = 0.2; q = 1 - p
# Both homozygous
m = marker(x, '3'=2, '4'=2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m), 1/4*q^2*(1+q)^2)
})
test_that("likelihoods are correct in 3 generat. SNP example", {
x = nuclearPed(1)
x = addParents(x, 1)
x = addParents(x, 2)
p = 0.2; q = 1 - p
# Both homozygous
m = marker(x, '4'=1:2, '7'=1:2, '3'=1:2, alleles=1:2, afreq=c(p,q))
expect_equal(likelihood(x, m), (2*p*q)^2 * (3/8 + p*q/2))
})
test_that("likelihoods are correct in looped ped", {
x = fullSibMating(1)
p = 0.1; q = 0.2; r = 1-p-q
m1 = marker(x, '5'=1:2, '6'=1:2, alleles=1:3, afreq=c(p,q,r))
m2 = marker(x, '5'=1:2, '6'=1:2, alleles=1:4, afreq=c(p,q,r/2,r/2))
expect_equal(likelihood(x, m1, verbose=F),
likelihood(x, m2, verbose=F))
m3 = marker(x, '5'=2, '6'=2, alleles=1:3, afreq=c(p,q,r))
m4 = marker(x, '5'=2, '6'=2, alleles=4:1, afreq=c(r/2,r/2,q,p))
expect_equal(likelihood(x, m3, verbose=F),
likelihood(x, m4, verbose=F))
})
# test_that("pedprobr and Familias give identical results", {
# skip("")
# library(forrel)
# library(pedlikCompare)
# testRandom = function() {
# x = randomPed(6)
# if(!is.ped(x)) return(NULL)
# ids = sample(labels(x), size = sample(pedsize(x), 1))
# nals = sample.int(5, size = 1) + 1
# afr = runif(nals, .2, .8)
# afr = afr/sum(afr)
# sims = markerSim(x, N=10, ids = ids, alleles = sample.int(10, size = nals),
# afreq = afr, mutmod = "prop", rate = 0.1, verbose = FALSE)
# for(i in 1:10) {
# tes = compare(sims, i, verbose = FALSE, programs = c("pedprobr", "Familias"))
# if(nrow(tes) < 2) break
# if(!all.equal(tes[1,2], tes[2,2])) {
# message("MISMATCH!")
# print(selectMarkers(x, i))
# }
# }
# }
# for(i in 1:100) testRandom()
# })
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.