Hi everyone,

For rate category variables, I'm looking for a way to combine a discretised gamma with equal-probability rate categories, (see code snippet), with a fixed-location discrete component.

This would be much like the "Gamma with a class of invariant sites" model, except I want to place the discrete component at a particular non-zero rate r and have probability p (so that the Gamma distribution is multiplied by 1-p).

I don't see any code in "Examples/Categories" that implements a Gamma+Invariant model, and my attempts to achieve this effect have bombed out with errors like "can't have a category variable depend on a category variable".

Does anyone have an idea on how to go about doing this? 


global alpha = 0.5;
alpha:>0.01;
alpha:<100;
beta:=alpha;

category c = (
    5,      /* number of rates */
    EQUAL,     /* probs. of rates */
    MEAN,     /* sampling method */
    GammaDist(_x_,alpha,beta),     /* density */
    CGammaDist(_x_,alpha,beta),    /* CDF */
    0,         /* left bound */
    1e5,       /* right bound */
    CGammaDist(_x_,alpha+1,beta)*alpha/beta /* "antiderivative" of x f(x) */
    );

Hi Sergei,

Thank you for the tips!  I wouldn't have thought to use the extra parameter to the LikelihoodFunction to carry out the mixture as a mixture of two complete models.

The resulting likelihoods behave well when comparing two mixtures (constraining some parameters reduces the likelihood, as expected). 

However, a simple (non-mixture) model like plain gamma or even plain point mass (!) comes out more likely than the mixture (by a few log-likelihood points), although the mixture is more general.  (point mass significantly less likely than gamma of course)

I have a few hypothesis, but after a lot of testing I suspect the mixture model is converging on some particular local optimum, or is converging too slowly to reach the true optimum in a reasonable step limit.

I think the right approach is to optimise params for the gamma model, then for the point-mass model (keeping the TV_TS from gamma?), then optimise P such that the sum of "Log(SITE_LIKELIHOOD[0]*P+(1-P)*SITE_LIKELIHOOD[1])" over sites is maximised.  However, if HyPhy is optimising everything simultaneously it could explain the local or slow convergence (?)

Does this sound right to you?  Would you like to see the test code?

To check that the model mixture isn't converging, I compared it to something conceptually equivalent but implemented without  model mixture.

Since I don't know how to make gamma+point mass using just a category variable (would like to know how, please!), I used gamma+invariant.  The first uses  5-cat category var (first being invariant sites), and the second a 4-cat gamma mixed with an invariant model.  The results are for P51.nex.   

For this code, the model-mixture gets a worse likelihood, and reject having invariant sites by setting IP=0 (!).   I also noticed that with gamma+point-mass the model mixture rejects putting a point-mass anywhere (P=0), and produces worse-than-gamma likelihoods.

What do you think is causing the above?

[code]
Gamma+Invariant using category variable
Log Likelihood = -3285.62637458753;
Shared Parameters:
IP=0.205486
alpha=0.369268
TV_TS=0.0989337
beta=alpha/(1-IP)=0.464773

Gamma+Invariant using model mixture
Log Likelihood = -3286.09788196153;
Shared Parameters:
IP=0
alpha=0.225829
TV_TS=0.0992244
beta=alpha/(1-IP)=0.225829
[/code]

With code being:

[code]
/********************************************************************
Read input data
********************************************************************/

/* Read in P51 HyPhy example alignment */
DataSet nucleotideSequences = ReadDataFile ( "p51.nex" );

/* Filter the data (all data, as nucleotides) */
DataSetFilter filteredData = CreateFilter ( nucleotideSequences, 1 );

/* Get observed nucleotide frequencies */
HarvestFrequencies ( observedFreqs, filteredData, 1, 1, 1 );

/********************************************************************
Gamma+Invariant LF using category variable
********************************************************************/

fprintf ( stdout, "Gamma+Invariant using category variable\n\n" );

global alpha = 1.0; alpha:>0.01; alpha:<100; /* shape parameter */
global beta = 1.0; beta:>0.01; beta:<100; /* beta paramter */
global TV_TS; TV_TS:>0; /* transversion/transition rate ratio */
global ncats = 5; ncats:=5; /* number of rate categories */
global IP = 1/ncats; IP :< 1-(1e-10); /* weight of invariant category */
beta:=alpha/(1-IP); /* mean of rate dist:=1 => (1-IP)*(beta/alpha)+IP*0 := 1 */

catFreqMatrix = { ncats, 1 };
catFreqMatrix[0] := IP;

for (h=1; h<ncats; h=h+1) {
    catFreqMatrix[h] := (1-IP)/(ncats-1);
}

category c = (
    ncats,                                                 /* number of rates */
    catFreqMatrix,                                         /* weights */
    MEAN,                                                  /* sampling method */
    (1-IP)*GammaDist(_x_,alpha,beta)*(_x_>0),              /* density */
    (1-IP)*CGammaDist(_x_,alpha,beta)*(_x_>0)+IP,          /* CDF */
    0,                                                     /* left bound */
    1e5,                                                   /* right bound */
    (1-IP)*CGammaDist(_x_,alpha+1,beta)*alpha/beta*(_x_>0) /* mean cumulative function */
    );

HKY85RateMatrixGamma = 
{{*        ,c*t*TV_TS,c*t      ,c*t*TV_TS}
{c*t*TV_TS,*        ,c*t*TV_TS,c*t      }
{c*t      ,c*t*TV_TS,*        ,c*t*TV_TS}
{c*t*TV_TS,c*t      ,c*t*TV_TS,*        }};

Model HKY85Gamma = ( HKY85RateMatrixGamma, observedFreqs );
Tree givenTree = DATAFILE_TREE;

LikelihoodFunction likeFunc = ( filteredData, givenTree );

Optimize ( paramValues, likeFunc );
fprintf ( stdout, likeFunc, "\n\n" );

/********************************************************************
Gamma+Invariant LF using model mixture
********************************************************************/

fprintf ( stdout, "Gamma+Invariant using model mixture\n\n" );

global alpha = 1.0; alpha:>0.01; alpha:<100; /* shape parameter */
global beta = 1.0; beta:>0.01; beta:<100; /* beta paramter */
global TV_TS; TV_TS:>0; /* transversion/transition rate ratio */
global IP = 0.5;  IP :< 1-(1e-10); /* Weight of invariant */
beta:=alpha/(1-IP); /* mean of rate dist:=1 => (1-IP)*(beta/alpha)+IP*0 := 1 */

category c = (
    ncats-1,                                      /* number of rates */
    EQUAL,                                  /* weights */
    MEAN,                                   /* sampling method */
    GammaDist(_x_,alpha,beta),              /* density */
    CGammaDist(_x_,alpha,beta),             /* CDF */
    0,                                      /* left bound */
    1e5,                                    /* right bound */
    CGammaDist(_x_,alpha+1,beta)*alpha/beta /* mean cumulative function */
    );

HKY85RateMatrixGamma = 
{{*        ,c*t*TV_TS,c*t      ,c*t*TV_TS}
{c*t*TV_TS,*        ,c*t*TV_TS,c*t      }
{c*t      ,c*t*TV_TS,*        ,c*t*TV_TS}
{c*t*TV_TS,c*t      ,c*t*TV_TS,*        }};

Model HKY85Gamma = ( HKY85RateMatrixGamma, observedFreqs );
Tree givenTree1 = DATAFILE_TREE;

RateMatrixInvariant = 
{{1,0,0,0}
{0,1,0,0}
{0,0,1,0}
{0,0,0,1}};

Model Invariant = ( RateMatrixInvariant, observedFreqs );
Tree givenTree2 = DATAFILE_TREE;

LikelihoodFunction likeFunc = ( filteredData, givenTree1, filteredData, givenTree2,
  "Log(SITE_LIKELIHOOD[0]*(1-IP)+SITE_LIKELIHOOD[1]*IP)" );

Optimize ( paramValues, likeFunc );
fprintf ( stdout, likeFunc, "\n\n" );

[/code]