********************************
Spencer, M. and Howe, C.J. (2004).  
Authenticity of ancient DNA results: a statistical approach.  
American Journal of Human Genetics 75: 240-250.

Matlab code for parameter estimation and confidence intervals

********************************

EXAMPLES OF USE

********************************

1.  One kind of control, we want the probability gamma that a positive result is at 
least partly from sample DNA:

[Pmax,mest,mnum,gconf]=aDNA_twokinds_llplot(nb,bb,ns,bs);

where the inputs are

nb is number of blank trials
bb is number of blank successses
ns is number of sample trials
bs is number of sample successes

and the outputs are

Pmax is max likelihood of data from exact ML estimate
mest is vector [theta2 gamma] exact ML estimates
mnum is vector [theta2 gamma] numerical ML estimates (rough, by grid not optimization)
gconf is vector [lo hi] of profile likelihood 95% CI for gamma

**********************************

2.  One kind of control, we want the probability gamma_2 that a positive result is 
entirely from sample DNA:

[Pmax,mest,mnum,gconf]=aDNA_twokinds_noc_llplot(nb,bb,ns,bs);

where the inputs are

nb is number of blank trials
b is number of blank successses
ns is number of sample trials
bs is number of sample successes

and the outputs are

Pmax is max likelihood of data from exact ML estimate
mest is vector [theta2 gamma2] exact ML estimates
mnum is vector [theta2 gamma2] numerical ML estimates (rough, by grid not optimization)
gconf is vector [lo hi] of profile likelihood 95% CI for gamma2

**********************************

3.  Two kinds of control, we want the probability gamma that a positive result is at
least partly from sample DNA:

[Pmax,mest,mnum,gconf]=aDNA_threekinds_llplot(ne,be,np,bp,ns,bs);

where the inputs are

ne is number of extraction blank trials
be is number of extraction blank successses
np is number of PCR blank trials
bp is number of PCR blank successes
ns is number of sample trials
bs is number of sample successes

and the outputs are

Pmax is max likelihood of data, using exact ML estimates
mest is vector [theta2 theta3 gamma] exact ML estimates
mnum is vector [theta2 theta3 gamma] numerical ML estimates (rough, by grid not optimization)
gconf is 95% profile confidence interval for gamma

**********************************

4.  Two kinds of control, we want the probability gamma_2 that a positive result is
entirely from sample DNA:

[Pmax,mest,mnum,gconf]=aDNA_threekinds_noc_llplot(ne,be,np,bp,ns,bs);

where the inputs are

ne is number of extraction blank trials
be is number of extraction blank successses
np is number of PCR blank trials
bp is number of PCR blank successes
ns is number of sample trials
bs is number of sample successes

and the outputs are

Pmax is max likelihood of data, using exact ML estimates
mest is vector [theta2 theta3 gamma2] exact ML estimates
mnum is vector [theta2 theta3 gamma2] numerical ML estimates (rough, by grid not optimization)
gconf is 95% profile confidence interval for gamma2

**********************************

The other .m files in this archive are used by the functions above.


**********************************

OTHER NOTES

1.  The hypothesis test described in Spencer and Howe (2004), equation 10 for gamma=0.  
Gamma=0 is a boundary value, so the distribution of 2*log likelihood under the null hypothesis
will be the mixture 1/2 chi^2_0 + 1/2 chi^2_1 (where chi^2_i is the chi-square distribution with
i degrees of freedom).  The P-value is half the size of the P-value calculated under chi^2_1.  See 
Self and Liang (1987).  If one or other of the confidence bounds on gamma hits zero or one, the confidence
interval may end up being too wide.  Similar but more complicated effects will occur for the case of two kinds of control.

2.  The confidence intervals are calculated assuming that the chi-square distribution is a good approximation to the distribution of
twice the log likelihood ratio.  This may not be true for small samples, so a parametric bootstrap test might be good in these situations.

3.  When reporting results, use mest, not mnum as the parameter estimates.  The mest vector is found analytically, 
while mnum is found by evaluating a grid of parameter values and is just a rough check that
the analytical result is about right.

REFERENCES

Self, S.G. and Liang, K.-Y. (1987).
Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions.
Journal of the American Statistical Association 82:605-610.

Spencer, M. and Howe, C.J. (2004).  
Authenticity of ancient DNA results: a statistical approach.  
American Journal of Human Genetics 75: 240-250.


********************************

Tested under Matlab release 14 service pack 1
Please report problems to Matthew Spencer (matts@mathstat.dal.ca)
26 Aug 2005