Dear Jennifer,
Quote: I have a question about the output from the HYK85_CI.bf. I have modified this batch file so that the substitution rate matrix used is the one chosen by data monkey for my dataset (011010).
In which B replaces kappa_inv everywhere in the batch file.
The problem is the output values for B here are the inverse of my estimates of B determined by the sample.bf in which the rate matrix has been replaced by 011010.
Is it correct to assume that the output from HKY_CI.bf is the inverse of the paramater being estimated? And so I should take the inverse of the estimates of "B" (kappa_inv) to have the true value of B.
I wrote HKY_CI.bf so that the model is parameterized by the inverse of kappa (i.e. the transversion/transition ratio in place of the more common transition/transversion ratio). This is done for technical (speed of convergence) reasons. However the file 'undoes' the reciprocation when it prints out the output. In your case, you probably don't want to inverse 'B' when it's output. You will notice the last two lines in HKY_CI.bf do the inversion:
Code:fprintf (stdout, "\n\nMLE Transition/Transversion estimate: ", 1/kappa_inv, "\nProfile likelihood 95% CI :", Format (1/cmx[2], 10, 5), " - ", Format(1/cmx[0], 10,5), "\n");
fprintf (stdout, "\n\nBS Transition/Transversion mean: ", bs_iterates/meanBS[0], "\nBootstrap 95% CI :", Format (1/upperBS, 10, 5), " - ", Format(1/lowerBS, 10,5), "\n");
.
In your case you should probably write:
Code:fprintf (stdout, "\n\nMLE B estimate: ", B, "\nProfile likelihood 95% CI :", Format (cmx[0], 10, 5), " - ", Format(cmx[2], 10,5), "\n");
fprintf (stdout, "\n\nBS B mean: ", meanBS[0]/bs_iterates, "\nBootstrap 95% CI :", Format (lowerBS, 10, 5), " - ", Format(upperBS, 10,5), "\n");
HTH,
Sergei