View Javadoc

    1   package baseCode.math.metaanalysis;
    2
    3   import baseCode.math.Constants;
    4   import cern.colt.list.DoubleArrayList;
    5   import cern.jet.stat.Descriptive;
    6   import cern.jet.stat.Probability;
    7
    8   /***
    9    * Statistics for meta-analysis. Methods from Cooper and Hedges (CH); Hunter and Schmidt (HS).
    10   * <p>
    11   * In this class "conditional variance" means the variance for one data set. Unconditional means "between data set", or
    12   * across data set.
    13   * <p>
    14   * Copyright (c) 2004
    15   * </p>
    16   * <p>
    17   * Institution: Columbia University
    18   * </p>
    19   *
    20   * @author Paul Pavlidis
    21   * @version $Id: MetaAnalysis.java,v 1.5 2005/03/21 18:01:03 pavlidis Exp $
    22   */
    23
    24  public abstract class MetaAnalysis {
    25
    26     /***
    27      * Fisher's method for combining p values. (Cooper and Hedges 15-8)
    28      *
    29      * @param pvals DoubleArrayList
    30      * @return double
    31      */
    32     protected double fisherCombinePvalues( DoubleArrayList pvals ) {
    33        double r = 0.0;
    34        for ( int i = 0, n = pvals.size(); i < n; i++ ) {
    35           r += Math.log( pvals.getQuick( i ) );
    36        }
    37        r *= -2.0;
    38        return Probability.chiSquare( r, 2.0 * pvals.size() );
    39     }
    40
    41     /***
    42      * Fisher's method for combining p values (Cooper and Hedges 15-8)
    43      * <p>
    44      * Use for p values that have already been log transformed.
    45      *
    46      * @param pvals DoubleArrayList
    47      * @return double
    48      */
    49     protected double fisherCombineLogPvalues( DoubleArrayList pvals ) {
    50        double r = 0.0;
    51        for ( int i = 0, n = pvals.size(); i < n; i++ ) {
    52           r += pvals.getQuick( i );
    53        }
    54        r *= -2.0;
    55        return Probability.chiSquare( r, 2.0 * pvals.size() );
    56     }
    57
    58     /***
    59      * The "Q" statistic used to test homogeneity of effect sizes. (Cooper and Hedges 18-6)
    60      *
    61      * @param effectSizes DoubleArrayList
    62      * @param variances DoubleArrayList
    63      * @param globalMean double
    64      * @return double
    65      */
    66     protected double qStatistic( DoubleArrayList effectSizes, DoubleArrayList variances, double globalMean ) {
    67
    68        if ( !( effectSizes.size() == variances.size() ) ) throw new IllegalArgumentException( "Unequal sizes" );
    69
    70        double r = 0.0;
    71        for ( int i = 0, n = effectSizes.size(); i < n; i++ ) {
    72           r += Math.pow( effectSizes.getQuick( i ) - globalMean, 2.0 ) / variances.getQuick( i );
    73        }
    74        return r;
    75     }
    76
    77     /***
    78      * Test for statistical significance of Q.
    79      *
    80      * @param Q - computed using qStatistic
    81      * @param N - number of studies.
    82      * @see qStatistic
    83      * @return The upper tail chi-square probability for Q with N - degrees of freedom.
    84      */
    85     public double qTest( double Q, double N ) {
    86        return Probability.chiSquareComplemented( N - 1, Q );
    87     }
    88
    89     /***
    90      * General formula for weighted mean of effect sizes. Cooper and Hedges 18-1, or HS pg. 100.
    91      * <p>
    92      * In HS, the weights are simply the sample sizes. For CH, the weights are 1/v for a fixed effect model. Under a
    93      * random effects model, we would use 1/(v + v_bs) where v_bs is the between-studies variance.
    94      *
    95      * @param effectSizes
    96      * @param sampleSizes
    97      * @return
    98      */
    99     protected double weightedMean( DoubleArrayList effectSizes, DoubleArrayList weights ) {
    100
    101       if ( !( effectSizes.size() == weights.size() ) ) throw new IllegalArgumentException( "Unequal sizes" );
    102
    103       double wm = 0.0;
    104       for ( int i = 0; i < effectSizes.size(); i++ ) {
    105          wm += weights.getQuick( i ) * effectSizes.getQuick( i );
    106       }
    107
    108       double s = Descriptive.sum( weights );
    109
    110       if ( s == 0.0 ) return 0.0;
    111       return wm /= s;
    112    }
    113
    114    /***
    115     * General formula for weighted mean of effect sizes including quality index scores for each value. Cooper and Hedges
    116     * 18-1, or HS pg. 100.
    117     *
    118     * @param effectSizes
    119     * @param sampleSizes
    120     * @param qualityIndices
    121     * @return
    122     */
    123    protected double weightedMean( DoubleArrayList effectSizes, DoubleArrayList weights, DoubleArrayList qualityIndices ) {
    124
    125       if ( !( effectSizes.size() == weights.size() && weights.size() == qualityIndices.size() ) )
    126             throw new IllegalArgumentException( "Unequal sizes" );
    127
    128       double wm = 0.0;
    129       for ( int i = 0; i < effectSizes.size(); i++ ) {
    130          wm += weights.getQuick( i ) * effectSizes.getQuick( i ) * qualityIndices.getQuick( i );
    131       }
    132       return wm /= ( Descriptive.sum( weights ) * Descriptive.sum( qualityIndices ) );
    133    }
    134
    135    /***
    136     * CH 18-3. Can be used for fixed or random effects model, the variances just have to computed differently.
    137     *
    138     * <pre>
    139     *          v_dot = 1/sum_i=1&circ;k ( 1/v_i)
    140     * </pre>
    141     *
    142     * @param variances
    143     * @return
    144     */
    145    protected double metaVariance( DoubleArrayList variances ) {
    146       double var = 0.0;
    147       for ( int i = 0; i < variances.size(); i++ ) {
    148          var += 1.0 / variances.getQuick( i );
    149       }
    150       if ( var == 0.0 ) {
    151          var = Double.MIN_VALUE;
    152          //   throw new IllegalStateException( "Variance of zero" );
    153       }
    154       return 1.0 / var;
    155    }
    156
    157    /***
    158     * CH 18-3 version 2 for quality weighted. ( page 266 ) in Fixed effects model.
    159     *
    160     * <pre>
    161     *
    162     *
    163     *
    164     *
    165     *
    166     *
    167     *
    168     *          v_dot = [ sum_i=1&circ;k ( q_i &circ; 2 * w_i) ]/[ sum_i=1&circ;k  q_i * w_i ]&circ;2
    169     *
    170     *
    171     *
    172     *
    173     *
    174     *
    175     *
    176     * </pre>
    177     *
    178     * @param variances
    179     * @return
    180     */
    181    protected double metaVariance( DoubleArrayList weights, DoubleArrayList qualityIndices ) {
    182       double num = 0.0;
    183       double denom = 0.0;
    184       for ( int i = 0; i < weights.size(); i++ ) {
    185          num += Math.pow( weights.getQuick( i ), 2 ) * qualityIndices.getQuick( i );
    186          denom += Math.pow( weights.getQuick( i ) * qualityIndices.getQuick( i ), 2 );
    187       }
    188       if ( denom == 0.0 ) {
    189          throw new IllegalStateException( "Attempt to divide by zero." );
    190       }
    191       return num / denom;
    192    }
    193
    194    /***
    195     * Test statistic for H0: effectSize == 0. CH 18-5. For fixed effects model.
    196     *
    197     * @param metaEffectSize
    198     * @param metaVariance
    199     * @return
    200     */
    201    protected double metaZscore( double metaEffectSize, double metaVariance ) {
    202       return Math.abs( metaEffectSize ) / Math.sqrt( metaVariance );
    203    }
    204
    205    /***
    206     * CH sample variance under random effects model, equation 18-20
    207     *
    208     * @param
    209     * @return
    210     */
    211    protected double metaRESampleVariance( DoubleArrayList effectSizes ) {
    212       return Descriptive.sampleVariance( effectSizes, Descriptive.mean( effectSizes ) );
    213    }
    214
    215    /***
    216     * CH estimate of between-studies variance, equation 18-22, for random effects model.
    217     *
    218     * @param effectSizes
    219     * @param variances
    220     * @return
    221     */
    222    //   protected double metaREVariance( DoubleArrayList effectSizes,
    223    //         DoubleArrayList variances ) {
    224    //      return Math.max( metaRESampleVariance( effectSizes )
    225    //            - Descriptive.mean( variances ), 0.0 );
    226    //   }
    227    /***
    228     * CH equation 18-23. Another estimator of the between-studies variance s<sup>2 </sup> for random effects model.
    229     * This is non-zero only if Q is larger than expected under the null hypothesis that the variance is zero.
    230     *
    231     * <pre>
    232     *
    233     *
    234     *
    235     *    s&circ;2 = [Q - ( k - 1 ) ] / c
    236     *
    237     *
    238     *
    239     * </pre>
    240     *
    241     * where
    242     *
    243     * <pre>
    244     *
    245     *
    246     *
    247     *    c = Max(sum_i=1&circ;k w_i - [ sum_i&circ;k w_i&circ;2 / sum_i&circ;k w_i ], 0)
    248     *
    249     *
    250     *
    251     * </pre>
    252     *
    253     * @param effectSizes
    254     * @param variances
    255     * @param weights
    256     * @return
    257     */
    258    protected double metaREVariance( DoubleArrayList effectSizes, DoubleArrayList variances, DoubleArrayList weights ) {
    259
    260       if ( !( effectSizes.size() == weights.size() && weights.size() == variances.size() ) )
    261             throw new IllegalArgumentException( "Unequal sizes" );
    262
    263       // the weighted unconditional variance.
    264       double q = qStatistic( effectSizes, variances, weightedMean( effectSizes, weights ) );
    265
    266       double c = Descriptive.sum( weights ) - Descriptive.sumOfSquares( weights ) / Descriptive.sum( weights );
    267       return Math.max( ( q - ( effectSizes.size() - 1 ) ) / c, 0.0 );
    268    }
    269
    270    /***
    271     * Under a random effects model, CH eqn. 18-24, we replace the conditional variance with the sum of the
    272     * between-sample variance and the conditional variance.
    273     *
    274     * <pre>
    275     *
    276     *
    277     *
    278     *    v_i&circ;* = sigma-hat_theta&circ;2 + v_i.
    279     *
    280     *
    281     *
    282     * </pre>
    283     *
    284     * @param variances Conditional variances
    285     * @param sampleVariance estimated...somehow.
    286     * @return
    287     */
    288    protected DoubleArrayList metaREWeights( DoubleArrayList variances, double sampleVariance ) {
    289       DoubleArrayList w = new DoubleArrayList( variances.size() );
    290
    291       for ( int i = 0; i < variances.size(); i++ ) {
    292          if ( variances.getQuick( i ) <= 0 ) {
    293             throw new IllegalStateException( "Negative or zero variance" );
    294          }
    295          w.add( 1 / ( variances.getQuick( i ) + sampleVariance ) );
    296       }
    297
    298       return w;
    299    }
    300
    301    /***
    302     * Weights under a fixed effects model. Simply w_i = 1/v_i. CH eqn 18-2.
    303     * @param variances
    304     * @return
    305     */
    306    protected DoubleArrayList metaFEWeights( DoubleArrayList variances ) {
    307       DoubleArrayList w = new DoubleArrayList( variances.size() );
    308
    309       for ( int i = 0; i < variances.size(); i++ ) {
    310          double v = variances.getQuick( i );
    311          if ( v <=  Constants.SMALL) {
    312             v = Constants.SMALL;
    313             System.err.println( "Tiny variance    " + v );
    314          }
    315          w.add( 1 / v );
    316       }
    317
    318       return w;
    319    }
    320 }