code p-value calculation

  pvalue/p-value calculation

HYPERGEOMETRIC TEST

http://mengnote.blogspot.com/2012/12/calculate-correct-hypergeometric-p.html

# x, q     vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
# m     the number of white balls in the urn.
# n     the number of black balls in the urn.
# k     the number of balls drawn from the urn.
# phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)

sucessDrawn=1
sucessTotal=2

sampleSize=24
populatonSizeWithoutSample=1751-sampleSize

phyper(sucessDrawn-1, sucessTotal, populatonSizeWithoutSample, sampleSize, lower.tail = FALSE, log.p = FALSE)

FISHER EXACT TEST

Bonferroni correction: Pvalue* (#of test cases)

http://www.une.edu.au/WebStat/unit_materials/c7_anova/oneway_bonferroni_adjust.htm

http://www.langsrud.com/stat/fisher.htm

Code done by  Benoit Marchand

/*
 * pvalcalc.cpp
 *      Author:  Benoit Marchand
 */
#include
#include


    /*
     *                 CLASSA    CLASSB
     *                        -----------------
     *  PRESENT    |    k    |        |    TOTALpresent= n
     *  ABSENT      |          |        |    TOTALabsent
     *                         ------------------------------------------------
     *     K= TOTALclassA|        |  N = TOTALclass A + B
     * 

/*
 *  k - #PRESENT in CLASSA
 *  n - #PRESENT in CLASSA + CLASSB
 *
 *  K - #TOTAL sample in CLASSA
 *  N - #TOTAL sample in CLASSA + CLASSB
 */
double calcOneMotifPValue(int k,int n,int K,int N)
{
        double pv,t1,x;

/*
        k      n-k         |  n
        K-k    (N-K)-(n-k) |  N-n
        -------------------------
        K      N-K         |  N


        This calculates the Right Hand Side PValue.
        (LHS would be k=0,n.)
*/
        t1=lgamma(((float) K+1))+lgamma(((float) N-K+1))-lgamma(((float) N+1))+lgamma(((float) n+1))+lgamma(((float) N-n+1));
        for (pv=-1e100;k<=n && k<=K && (n-k)<=(N-K);k++) {
                x=t1-lgamma(((float) k+1))-lgamma(((float) K-k+1))-lgamma(((float) n-k+1))-lgamma(((float) N-K-n+k+1));
                pv=(pv>=x?pv+log(1.0+exp(x-pv)):x+log(1.0+exp(pv-x)));
        }

        // Bonferroni correction (pval * nSEQF)
        //rankDB[id].pvalue=((float) rankDB[id].nSEQF)*exp(pv);
        //rankDB[id].pvalue=exp(pv);
        return(exp(pv));
}


int main()
{
    printf("%.8f", calcOneMotifPValue( 75,100, 155, 200) );
}






     *                 CLASSA    CLASSB
     *                        -----------------
     *  PRESENT    |    k    |        |    TOTALpresent=n
     *  ABSENT      |          |        |    TOTALabsent
     *                         -------------------------------
     *         TOTALclassA  |        |      TOTALclass A + B
     *                 K                                    N
     */

    /*
     *  k - #PRESENT in CLASSA
     *  n - #PRESENT in CLASSA + CLASSB
     *
     *  K - #TOTAL sample in CLASSA
     *  N - #TOTAL sample in CLASSA + CLASSB
     */
    public static double calcRightsidePValue(int hitA,int hitAB,int totA,int totAB)
    {
            double pv,t1,x;

    /*
            k      n-k         |  n
            K-k    (N-K)-(n-k) |  N-n
            -------------------------
            K      N-K         |  N


            This calculates the Right Hand Side PValue.
            (LHS would be k=0,n.)
    */
            t1= Gamma.logGamma(((float) totA+1))+Gamma.logGamma(((float) totAB-totA+1))-Gamma.logGamma(((float) totAB+1))+Gamma.logGamma(((float) hitAB+1))+Gamma.logGamma(((float) totAB-hitAB+1));
            for (pv=-1e100;hitA<=hitAB && hitA<=totA && (hitAB-hitA)<=(totAB-totA);hitA++) {
                    x=t1-Gamma.logGamma(((float) hitA+1))-Gamma.logGamma(((float) totA-hitA+1))-Gamma.logGamma(((float) hitAB-hitA+1))-Gamma.logGamma(((float) totAB-totA-hitAB+hitA+1));
                   
                    pv=( pv>=x  ? pv+  Math.log(1.0+Math.exp(x-pv))  :  x+ Math.log(1.0+ Math.exp(pv-x)));
            }

            // Bonferroni correction (pval * nSEQF)
            //rankDB[id].pvalue=((float) rankDB[id].nSEQF)*exp(pv);
            //rankDB[id].pvalue=exp(pv);
            return( Math.exp(pv));
    }

Fisher 's exact test in R package: stats
===========================
myContin = matrix(c(1,4,7,4), nrow = 2);
fisher.test(myContin ,alternative = "two.sided" ) ; // two sided

fisher.test(myContin ,alternative = "greater" ) ; // right-hand sided
fisher.test(myContin ,alternative = "less" ) ; // left-hand sided