Hi Alex,
That's an awful lot of parameters to estimate in your model --- I hope you have tons of data!   It would make it easier to answer your question if you could be clearer about your variables.  I'm guessing that you're fitting a mixture of codon models crossed with HKY85, with stationary nucleotide frequencies to be estimated (else, why constrain your global variables to 1?). I'm also guessing that your 20 category-specific variables correspond to amino-acid stationary distributions that modify the codon substitution matrix, modeling variation over sites.

You can certainly carry out this sort of analysis using HyPhy.  The simplest approach would be to define your five separate codon substitution models and then assign all five to the same likelihood function with your "mixture" parameters to specify how to calculate site likelihood.  I'll post an example in a bit.

- Art.

Thank you, that looks very much like what I want to do.

This batch file seems to estimate a different set of branch lengths for each rate matrix.  How do I make it estimate only a single set of branch lengths?  Also, how would I tell HyPhy to use the branch lengths specified in the nexus file rather than estimating them?

Thanks,
Alex

After making those changes (at least my interpretation of those changes), the branch lengths of the two trees are still different, and they are different from those of the reference tree.  However, they appear to be different from the reference tree by only a scaling factor that is the same for each branch within a tree.

Are these scaling factors being estimated?  If so, how can I fix them to 1.0?  On the other hand if they are deterministic consequences of, for example, normalizing the expected rate of the component rate matrices, then how exactly are they calculated?

I have a separate question about how the likelihood is maximized.  Is HyPhy taking advantage of the 'mixture' structure of my likelihood function as it searches for the MLE?  It is possible to analytically find the maximum likelihood values of some of the parameters in the model conditional on the the data and the values of the other parameters.

Thanks,
Alex



[code]
VERBOSITY_LEVEL = 11;

DataSet                         spectrinData = ReadDataFile ("sim2.nex");

DataSetFilter           filteredData = CreateFilter (spectrinData,1);

HarvestFrequencies (observedFreqs, filteredData, 1, 1, 1);

fprintf (stdout, observedFreqs, "\n");

/* stationary distribution in nucleotide frequencies */
global  eqFreqA = 0.25;
global  eqFreqC = 0.25;
global  eqFreqG = 0.25;
global  eqFreqT := 1.0 - eqFreqA - eqFreqC - eqFreqG;   eqFreqT :> 0;

/* another distribution for nucleotide frequencies */
global  eqFreqA2 = 0.25;
global  eqFreqC2 = 0.25;
global  eqFreqG2 = 0.25;
global  eqFreqT2 := 1.0 - eqFreqA2 - eqFreqC2 - eqFreqG2;       eqFreqT2 :> 0;

/* transition/transversion ratio */
global  kappa = 1.0;    kappa :> 0;

HKY85RateMatrix =
               {{*,a,kappa*a,a}
                {a,*,a,kappa*a}
                {kappa*a,a,*,a}
                {a,kappa*a,a,*}};

/* construct vectors from global variables */
estFreqs1 = {
                       {eqFreqA,
                       eqFreqC,
                       eqFreqG,
                       eqFreqT}
                       };

estFreqs2 = {
                       {eqFreqA2,
                       eqFreqC2,
                       eqFreqG2,
                       eqFreqT2}
                       };
/* assign models to trees */
Model   m1 = (HKY85RateMatrix, estFreqs1);

ACCEPT_BRANCH_LENGTHS = 1;
Tree    givenTree       = DATAFILE_TREE;
branchNames  = BranchName (givenTree, -1);  /* vector populated with names for each branch in tree */
branchLengths = BranchLength (givenTree, -1);  /* vector populated with branch lengths */
for (k = 0; k <  Columns(branchNames)-1; k = k+1)
{
  ExecuteCommands("givenTree."+branchNames[k]+".a:="+branchLengths[k]+";");
}

Model   m2 = (HKY85RateMatrix, estFreqs2);

ACCEPT_BRANCH_LENGTHS = 1;
Tree    otherTree       = DATAFILE_TREE;
branchNames  = BranchName (otherTree, -1);  /* vector populated with names for each branch in tree */
branchLengths = BranchLength (otherTree, -1);  /* vector populated with branch lengths */
for (k = 0; k <  Columns(branchNames)-1; k = k+1)
{
  ExecuteCommands("otherTree."+branchNames[k]+".a:="+branchLengths[k]+";");
}

global  P = 0.5; /* mixing proportion */
P :> 0; P :< 1;

/* construct likelihood function using mixture of two rate distributions */
LikelihoodFunction      lf = (filteredData, givenTree, filteredData, otherTree, "Log(SITE_LIKELIHOOD[0]*P+(1-P)*SITE_LIKELIHOOD[1])");

Optimize(res, lf);

fprintf(stdout, res);
fprintf(stdout, lf);

[/code]

In response to the second part of your message:

HyPhy is computing likelihood in the mixture model for each site of the alignment using the formula provided as the fifth argument in LikelihoodFunction() so, as you say, HyPhy is using the mixture structure of the likelihood function in its model fitting procedure.  Though I have to admit that I find a mixture of stationary frequency distributions rather odd. 

I don't understand what you're asking for, with regards to finding MLEs conditional on the values of other parameters.  I suppose you could either implement meta-parameters in the model to represent conditional dependence between the existing parameters.

Best,
- Art.

Hi Alex,

You can constrain the mixture component parameters by the following:

Code:

P := 0.5;
 

As for your example, HyPhy is reporting only 4 free parameters, namely three nucleotide frequencies and kappa.  I'm not picking up this extra parameter you're talking about    Scaling of the tree to the expected number of substitutions should not require any additional parameters, as this is only a function of the stationary distribution of nucleotide frequencies and parameters in the rate matrix.

Cheers,
- Art.

Hi Alex,

Branch lengths are scaled in units of expected number of substitutions, which is obtained from the substitution rate matrix and is initially equal to 1.  Though you've constrained branch lengths to values specified in the data file, the units are being rescaled as parameter estimates (nucleotide frequencies, kappa) are updated.  Does that make sense?

Cheers,
- Art.