Transition matrix with diagonals equal

49 views
Skip to first unread message

Anne-Marie L

unread,
Apr 12, 2017, 11:51:29 AM4/12/17
to revbayes-users
Hi everyone,

I would like to model microsatellite data, i.e. my data consists of lengths (integers) of N microsatellite loci for each of M species.

I can define a transition matrix for the simplest model, in which microsatellites can only change length by + or - 1, by using a rate matrix with zeros everywhere other than the + or - 1 diagonals. My question is: how can I define a more complex model in which, for example, jumps of a larger size can occur with some probability dependent on the size of the jump? This would correspond to a rate matrix where each diagonal has a different value, something like the following:

0  a  b  c
a  0  a  b
b  a  0  b
c  b  a  0

Thanks in advance,
Anne-Marie

Michael Landis

unread,
Apr 12, 2017, 12:18:35 PM4/12/17
to Anne-Marie L, revbayes-users
Hi Anne-Marie,

Here's one way to accomplish that:

```
# size of the nxn rate matrix
n <- 5

# n-1 rate classes
for (i in 1:(n-1)) {
    k[i] ~ dnExp(1)
    k[i].setValue(n-i) # pretty starting values
}

# populate the transition rates:
for (i in 1:n) {
    for (j in i:n) {
        if (i != j) {
            r[i][j] := k[ j-i ]
            r[j][i] := k[ j-i ]
        } else {
            r[i][i] <- 0.
        }
    }
}

# create the rate matrix
Q := fnFreeK(r, rescaled=false)
```

You'll likely notice that performance degrades as the size of Q gets large, so you have to be a little judicious in choosing n with respect to the variation you observe in your microsats.

Also, you may try alternative parameterizations to improve mixing if you notice problems in practice. Two things that might help are 1) setting rescaled=true in fnFreeK and using an external clock parameter, 2) defining k as a Dirichlet rather than a vector of exponentials (although don't believe this will result in a "flat" rate matrix without special considerations).

Hope this helps!
Michael
--
You received this message because you are subscribed to the Google Groups "revbayes-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to revbayes-user...@googlegroups.com.
To post to this group, send email to revbaye...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/revbayes-users/49fead67-fd20-4dc5-b185-15c68a4e17e3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Anne-Marie L

unread,
Apr 13, 2017, 4:02:49 AM4/13/17
to revbayes-users
Yes, that helps a lot! Thanks so much. It was the option 'rescaled=FALSE' that I had missed. I'll try out the various alternatives and let you know which works best.

Anne-Marie
Reply all
Reply to author
Forward
0 new messages