Changing mutation model
thoree opened this issue · 3 comments
Below I use attr
to change mutation rates. This leaves
the mutation matrix unchanged so it doesn't make sense
and should perhap not be allowed:
library(pedtools)
locus = list(list(name = "M", alleles = 1:2, mutmod = "equal", rate = 0.05))
x = setMarkers(nuclearPed(1), locus = locus)
attr(mutmod(x, marker = 1)[[1]], "rate" ) = 0.0
attr(mutmod(x, marker = 1)[[2]], "rate" ) = 0.0
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#> 1 2
#> 1 0.95 0.05
#> 2 0.05 0.95
#>
#> Model: equal
#> Rate: 0
#> Frequencies: 0.5, 0.5
#>
#> Stationary: Yes
#> Reversible: Yes
#> Lumpable: Always
# The likelihood performs fine as it only uses the mutation matrix:
amat = matrix(c(1,1,2,2), ncol = 2, byrow = T)
x = setAlleles(x, ids = c(1,3), markers = "M", alleles = amat)
pedprobr::likelihood(x,1)
#> [1] 0.00625
Created on 2019-07-31 by the reprex package (v0.3.0)
First of all, I don't think there is a simple/useful way to stop you from tinkering with the attributes of an R object.
Secondly, as you already noticed, it is the mutation matrix that matters in computations; the rate
parameter is mainly used to create this matrix. Changing the rate means changing the whole model; hence it makes sense (to me, at least) that the way to do this is to replace the whole model. (See code example below.)
Thirdly, in the particular case r=0
you could of course just remove the mutation model entirely.
library(pedtools)
# Attach marker locus with mutation model
locus = list(alleles = 1:2, mutmod = "equal", rate = 0.05)
x = setMarkers(nuclearPed(1), locus = list(locus))
# Setting a new rate
mutmod(x, marker = 1) = pedmut::mutationModel("equal", alleles = alleles(x, 1), rate = 0)
# Inspect it:
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#> 1 2
#> 1 1 0
#> 2 0 1
#>
#> Model: equal
#> Rate: 0
#>
#> Lumpable: Always
# Or simply remove the mutation model
mutmod(x, marker = 1) = NULL
Created on 2019-08-01 by the reprex package (v0.3.0)
Thanks. One could also go "trivial":
library(pedtools)
# Attach marker locus with mutation model
locus = list(alleles = 1:2, mutmod = "equal", rate = 0.05)
x = setMarkers(nuclearPed(1), locus = list(locus))
# Setting a new rate
mutmod(x, marker = 1) = "trivial"
# Inspect it:
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#> 1 2
#> 1 1 0
#> 2 0 1
#>
#> Model: trivial
#> Rate: NA
#>
#> Lumpable: Always
Created on 2019-08-02 by the reprex package (v0.3.0)
Ah, yes! I had completely forgotten about the "trivial" model. :)