mtc.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * SHORE is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


#ifndef SHORE_STATISTICS_MTC_HPP__
#define SHORE_STATISTICS_MTC_HPP__

#include <algorithm>
#include <cmath>
#include <iterator>
#include <vector>

#include "shore/base/stringops.hpp"

namespace shore {

class mtc
{

 public:
        
        enum mtcmethod
        {
                FDR_BH=0,
                FWER_BONFERRONI=1,
                FWER_HOLM=2,
                FWER_HOCHBERG=3,
                FWER_SIDAK_SS=4,
                FWER_SIDAK_SD=5,
                FDR_BY=6
        };

 private:

        std::vector<double> m_rawp;
        std::vector<size_t> m_rawp_ord;

        template<class Iter>
        void runBhCorrection(Iter res)
        {
                ensureOrdering();

                res[*m_rawp_ord.rbegin()]=m_rawp[*m_rawp_ord.rbegin()];
                for(int i=m_rawp_ord.size()-2;i>=0;--i)
                {
                        const size_t j=m_rawp_ord[i+1];
                        const size_t k=m_rawp_ord[i];

                        const double c=static_cast<double>(m_rawp_ord.size())/(i+1);
                        const double d=std::min(c*m_rawp[k],1.0);

                        res[k]=std::min(res[j],d);
                }
        }

        template<class Iter>
        void runBonferroniCorrection(Iter res)
        {
                for(size_t i=0;i<m_rawp.size();++i)
                {
                        res[i]=std::min(m_rawp[i]*m_rawp.size(),1.0);
                }
        }

        template<class Iter>
        void runHolmCorrection(Iter res)
        {
                ensureOrdering();

                const size_t m=m_rawp_ord.size();
                res[m_rawp_ord[0]]=std::min(m_rawp[m_rawp_ord[0]]*m,1.0);
                for(size_t i=1;i<m;++i)
                {
                        const double scaled=
                                std::min(m_rawp[m_rawp_ord[i]]*(m-i),1.0);

                        res[m_rawp_ord[i]]=
                                std::max(scaled,res[m_rawp_ord[i-1]]);
                }
        }

        template<class Iter>
        void runHochbergCorrection(Iter res)
        {
                ensureOrdering();

                const size_t m=m_rawp_ord.size();
                res[*m_rawp_ord.rbegin()]=m_rawp[*m_rawp_ord.rbegin()];

                for(int i=m-2;i>=0;--i)
                {
                        const double sc=std::min(m_rawp[m_rawp_ord[i]]*(m-i),1.0);
                        res[m_rawp_ord[i]]=
                                std::min(sc,res[m_rawp_ord[i+1]]);
                }
        }

        template<class Iter>
        void runSidakSsCorrection(Iter res)
        {
                const int m=m_rawp.size();

                for(size_t i=0;i<m_rawp.size();++i)
                {
                        res[i]=1.0-std::pow(1.0-m_rawp[i],m);
                }
        }

        template<class Iter>
        void runSidakSdCorrection(Iter res)
        {
                ensureOrdering();

                const int m=m_rawp_ord.size();
                res[m_rawp_ord[0]]=
                        1.0-std::pow(1.0-m_rawp[m_rawp_ord[0]],m);

                for(int i=1;i<m;++i)
                {
                        res[m_rawp_ord[i]]=
                                std::max(res[m_rawp_ord[i-1]],
                                         1.0-std::pow(1.0-m_rawp[m_rawp_ord[i]],m-i));
                }
        }

        template<class Iter>
        void runByCorrection(Iter res)
        {
                ensureOrdering();

                const size_t m=m_rawp_ord.size();

                double a=1.5;
                for(size_t i=3;i<=m;++i) a+=1.0/i;

                res[*m_rawp_ord.rbegin()]=
                        std::min(a*m_rawp[*m_rawp_ord.rbegin()],1.0);

                for(int i=m-2;i>=0;--i)
                {
                        const double sc=std::min(m*a*m_rawp[m_rawp_ord[i]]/(i+1),1.0);
                        res[m_rawp_ord[i]]=
                                std::min(res[m_rawp_ord[i+1]],sc);
                }
        }

        void ensureOrdering();

 public:

        mtc();

        template<class Iter>
        mtc(Iter rawp_start,Iter rawp_end)
        {
                std::copy(rawp_start,rawp_end,std::back_inserter(m_rawp));
        }

        void addRawValue(const double val);

        template<class Iter>
        void correctPvalues(const mtcmethod m,Iter result)
        {
                if(m_rawp.size()<=1)
                {
                        if(!m_rawp.empty()) *result=m_rawp[0];
                }
                else switch(m)
                {
                 case FDR_BH:
                        runBhCorrection(result);
                        break;

                 case FWER_BONFERRONI:
                        runBonferroniCorrection(result);
                        break;

                 case FWER_HOLM:
                        runHolmCorrection(result);
                        break;

                 case FWER_HOCHBERG:
                        runHochbergCorrection(result);
                        break;

                 case FWER_SIDAK_SS:
                        runSidakSsCorrection(result);
                        break;

                 case FWER_SIDAK_SD:
                        runSidakSdCorrection(result);
                        break;

                 case FDR_BY:
                        runByCorrection(result);
                        break;

                 default:
                        throw "unhandled p-value correction method";
                }
        }

        size_t size() const;
};

S_MAKE_ENUM_TRAITS(mtc::mtcmethod)

} // namespace

#endif // SHORE_STATISTICS_MTC_HPP__