Define two likelihood functions We are now ready to set up LikelihoodFunctions and Optimize them. The analysis of the combined data set is routine: LikelihoodFunction theLikFun = (myFilter,myTree); Optimize (paramValues, theLikFun); We also store some results for later use: lnLik0 = paramValues[1][0]; npar0 = paramValues[1][1]+3; fprintf (stdout, theLikFun, "\n\n"); The statement npar0 = paramValues[1][1]+3; requires some explanation. The Optimize function always returns the number of parameters that were optimized as the [1][1] element of its returned matrix of results. Typically, we do not optimize over base frequency values, electing instead to simply use observed frequencies, which are usually very close to the maximum likelihood estimates. Since the frequencies are, in fact, estimated from the data, they need to be included in the parameter count. The value of npar0, therefore, includes the count of independent substitution parameters in the model (the number of which is returned by Optimize) along with the three independent base frequencies estimated from the data. The LikelihoodFunction for the ``partitioned'' analysis simply uses the extended form of the LikelihoodFunction command: LikelihoodFunction theSplitLikFun = (myFilter1,myTree1, myFilter2,myTree2, myFilter3,myTree3); Optimize (paramValues, theSplitLikFun); lnLik1 = paramValues[1][0]; npar1 = paramValues[1][1]+9; Note the addition of the 9 estimated frequencies to the model's parameter count. Finally, we compute the P-value for the test of the combined analysis (null hypothesis) against the split model (alternative hypothesis). Two approaches are used. First, the normal chi-squared approximation to the LRT statistic: LRT = 2*(lnLik1-lnLik0); pValueChi2 = 1-CChi2 (LRT, npar1-npar0).
Sergei L. Kosakovsky Pond and Spencer V. Muse, 1997-2002