So... I realise now that this rev model forces equilibrium frequencies of 0.5 for both 1 and 0 characters - (for two characters with the rate of change the same from 0 to 1 as 1 to 0, the equilibrium frequency has to be 0.5 for both of them...)
I found that including "t" into the substitution matrix, and not defining "t=1" initially, gives me a tree whose branch lengths don't change as I change the rate parameter (which is what I'm looking for - I want to optimise the rate parameter without changing branch lenghts or topology).
To get a tree whose branch lengths don't change after "Optimize" (i.e. branches that stay the same for different values of the rate parameter), I have to place the constraints on branch lengths AFTER optimising the likelihood function
Somehow I'm surprised that the constraints need to be placed afterwards like this - I'd have thought they need to be constrained BEFORE doing the optimisation - something I'm missing here...
I realise I just don't get what "t" is doing here, and how it makes a difference whether or not it is the description of Q.
Also, in the script below, I find that changing the rate has no effect on the likelihood of the model - not sure whether this is because I'm somehow asking a nonsensical question here, am coding things wrong, or that I'm just dealing with a likelihood that doesn't depend on rate....??
Any hints/suggestions on these issues much appreciated.
Sorry for all the Newbie questions - I hope if the conceptual problems I seem to have can be ironed out, I'll be able to sort stuff out much better on my own...?!
Aidan
Code:dssRev = "$BASESET:01\n#a\n101110101111111\n#b\n101011000101111\n#c\n011001111001000\n#d\n111001110100001";
DataSet myDataRev = ReadFromString(dssRev);
DataSetFilter myFilterRev = CreateFilter (myDataRev,1);
ACCEPT_BRANCH_LENGTHS = 1;
ACCEPT_ROOTED_TREES = 1;
LIKELIHOOD_FUNCTION_OUTPUT = 1;
global rateRev = 1;
/*tRev=1; t=1?*/
/*t=1;*/
fprintf (stdout, "\n\nStart New Analysis\n");
binaryRevQ = {{*,rateRev*t}
{rateRev*t,*}};
eqFrequenciesRev = {{0.5}{0.5}};
/*For a fixed topology and fixed set of branch lengths, get the ML estimate of the rate of a reversible model*/
/* With rev model, need to have final parameter set to 0...? */
Model binaryRev = (binaryRevQ, eqFrequenciesRev, 0);
Tree myTreeRev = (((a:0.1,b:0.3):0.1,c:0.1):0.1,d:1);
/*get all branch names and branch lengths*/
branchNamesRev = BranchName (myTreeRev,-1);
branchLengthsInitialRev = BranchLength (myTreeRev,-1);
LikelihoodFunction lfRev = (myFilterRev,myTreeRev);
fprintf (stdout, "\nLikelihood function before optimisation\n",lfRev);
treeString = Format (myTreeRev, 1, 1);
fprintf (stdout, "tree string before optimisation: ", treeString, "\n");
Optimize (res,lfRev);
for (k = 0; k < Columns (branchNamesRev)-1; k=k+1)
{
outString = "myTreeRev."+branchNamesRev[k]+".t:="+branchLengthsInitialRev[k]+"/rateRev;";
ExecuteCommands (outString);
fprintf (stdout, "executing: ", outString, "\n");
}
fprintf (stdout, "\nLikelihood function after optimisation\n",lfRev);
treeString = Format (myTreeRev, 1, 1);
fprintf (stdout, "\ntree string with single timeRev rate:",rateRev," ", treeString, "\n");
rateRev = 0.01;
fprintf (stdout, "\nLikelihood function with rateRev set to: ",rateRev, "\n",lfRev);
treeString = Format (myTreeRev, 1, 1);
fprintf (stdout, "\ntree string with single timeRev rate:",rateRev," ", treeString, "\n");