magnusdv / pedprobr

Pedigree probabilities in R
GNU General Public License v2.0
4 stars 0 forks source link

`twoMarkerDistribution` gives unexpected result #15

Closed mkruijver closed 1 year ago

mkruijver commented 1 year ago

For the example below an unexpected result is obtained from twoMarkerDistribution

require(pedtools)
#> Loading required package: pedtools
ped <- nuclearPed(nch = 2) |> addMarker() |> addMarker()

pedprobr::twoMarkerDistribution(ped, id = c("1","2"), 
                                partialmarker1 = 1, partialmarker2 = 2, 
                                rho = 0.5)
#> Known genotypes:
#>  id <1> <2>      
#>   1 -/- -/-  <===
#>   2 -/- -/-  <===
#>   3 -/- -/-      
#>   4 -/- -/-      
#> 
#> Allele frequencies, marker 1:
#>    1   2
#>  0.5 0.5
#> 
#> Allele frequencies, marker 2:
#>    1   2
#>  0.5 0.5
#> 
#> Recombination rate : 0.5
#> Chromosome type    : autosomal
#> Target individual  : 1 2 
#> Marginal likelihood: 1 
#> Calculations needed: 81 
#> 
#> Analysis finished in 0.0557 secs
#>            1/1 1/2 2/2
#> 1/1 0.00390625   0   0
#> 1/2 0.00390625   0   0
#> 2/2 0.00390625   0   0

# compute the expected result
p1 <- pedprobr::oneMarkerDistribution(ped, ids = 1, partialmarker = 1, verbose = FALSE)
p2 <- pedprobr::oneMarkerDistribution(ped, ids = 2, partialmarker = 2, verbose = FALSE)

outer(outer(p1, p2),outer(p1, p2))
#> , , 1/1, 1/1
#> 
#>            1/1       1/2        2/2
#> 1/1 0.00390625 0.0078125 0.00390625
#> 1/2 0.00781250 0.0156250 0.00781250
#> 2/2 0.00390625 0.0078125 0.00390625
#> 
#> , , 1/2, 1/1
#> 
#>           1/1      1/2       2/2
#> 1/1 0.0078125 0.015625 0.0078125
#> 1/2 0.0156250 0.031250 0.0156250
#> 2/2 0.0078125 0.015625 0.0078125
#> 
#> , , 2/2, 1/1
#> 
#>            1/1       1/2        2/2
#> 1/1 0.00390625 0.0078125 0.00390625
#> 1/2 0.00781250 0.0156250 0.00781250
#> 2/2 0.00390625 0.0078125 0.00390625
#> 
#> , , 1/1, 1/2
#> 
#>           1/1      1/2       2/2
#> 1/1 0.0078125 0.015625 0.0078125
#> 1/2 0.0156250 0.031250 0.0156250
#> 2/2 0.0078125 0.015625 0.0078125
#> 
#> , , 1/2, 1/2
#> 
#>          1/1     1/2      2/2
#> 1/1 0.015625 0.03125 0.015625
#> 1/2 0.031250 0.06250 0.031250
#> 2/2 0.015625 0.03125 0.015625
#> 
#> , , 2/2, 1/2
#> 
#>           1/1      1/2       2/2
#> 1/1 0.0078125 0.015625 0.0078125
#> 1/2 0.0156250 0.031250 0.0156250
#> 2/2 0.0078125 0.015625 0.0078125
#> 
#> , , 1/1, 2/2
#> 
#>            1/1       1/2        2/2
#> 1/1 0.00390625 0.0078125 0.00390625
#> 1/2 0.00781250 0.0156250 0.00781250
#> 2/2 0.00390625 0.0078125 0.00390625
#> 
#> , , 1/2, 2/2
#> 
#>           1/1      1/2       2/2
#> 1/1 0.0078125 0.015625 0.0078125
#> 1/2 0.0156250 0.031250 0.0156250
#> 2/2 0.0078125 0.015625 0.0078125
#> 
#> , , 2/2, 2/2
#> 
#>            1/1       1/2        2/2
#> 1/1 0.00390625 0.0078125 0.00390625
#> 1/2 0.00781250 0.0156250 0.00781250
#> 2/2 0.00390625 0.0078125 0.00390625

Created on 2023-09-11 with reprex v2.0.2

mkruijver commented 1 year ago

In hindsight this is probably caused by id not being of length one. twoMarkerDistribution does not support this case but does not throw an error.

magnusdv commented 1 year ago

You are absolutely right; twoMarkerDistribution() currently only works for a single individual, and should definitely throw an error otherwise.

Thanks for exposing it.