Java Code Examples for org.apache.commons.math3.stat.regression.SimpleRegression#getR()

public static void main(String[] args) {
    SimpleRegression sr = getData("data/Data1.dat");
    double m = sr.getSlope();
    double b = sr.getIntercept();
    double r = sr.getR();  // correlation coefficient
    double r2 = sr.getRSquare();
    double sse = sr.getSumSquaredErrors();
    double tss = sr.getTotalSumSquares();

    System.out.printf("y = %.6fx + %.4f%n", m, b);
    System.out.printf("r = %.6f%n", r);
    System.out.printf("r2 = %.6f%n", r2);
    System.out.printf("EV = %.5f%n", tss - sse);
    System.out.printf("UV = %.4f%n", sse);
    System.out.printf("TV = %.3f%n", tss);
}
 

public static double correlation(IDoubleVector x, IDoubleVector y) {
    final int len = x.getLength();
    if (len != y.getLength()) {
        throw new AdeCoreIllegalArgumentException("Mismatching lengths");
    }
    if (len < 2) {
        throw new AdeCoreIllegalArgumentException("Vectors must have length >=2");
    }
    final SimpleRegression regression = new SimpleRegression();
    for (int i = 0; i < len; i++) {
        regression.addData(x.get(i), y.get(i));
    }
    return regression.getR();
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between two arrays.
 *
 * <p>Throws MathIllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2.  Returns {@code NaN} if either of the arrays
 * has zero variance (i.e., if one of the arrays does not contain at least two distinct
 * values).</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between the two arrays.
 *
 * </p>Throws IllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

/**
 * Computes the Pearson's product-moment correlation coefficient between two arrays.
 *
 * <p>Throws MathIllegalArgumentException if the arrays do not have the same length
 * or their common length is less than 2.  Returns {@code NaN} if either of the arrays
 * has zero variance (i.e., if one of the arrays does not contain at least two distinct
 * values).</p>
 *
 * @param xArray first data array
 * @param yArray second data array
 * @return Returns Pearson's correlation coefficient for the two arrays
 * @throws DimensionMismatchException if the arrays lengths do not match
 * @throws MathIllegalArgumentException if there is insufficient data
 */
public double correlation(final double[] xArray, final double[] yArray) {
    SimpleRegression regression = new SimpleRegression();
    if (xArray.length != yArray.length) {
        throw new DimensionMismatchException(xArray.length, yArray.length);
    } else if (xArray.length < 2) {
        throw new MathIllegalArgumentException(LocalizedFormats.INSUFFICIENT_DIMENSION,
                                               xArray.length, 2);
    } else {
        for(int i=0; i<xArray.length; i++) {
            regression.addData(xArray[i], yArray[i]);
        }
        return regression.getR();
    }
}
 

@Override
public void calculateStatistics() {

	double zscoreSum = 0;

	SimpleRegression regression = new SimpleRegression();

	for (int i = 0 ; i < a1Counts.size() ; ++i){

		regression.addData(a1Counts.getQuick(i), a2Counts.getQuick(i));

		final double pvalue = pValues.getQuick(i);

		// we used 2 sided test so divide by 2
		//double zscore = normalDist.inverseCumulativeProbability(pvalue/2);
		final double pvalueDiv2 = pvalue / 2;
		final double zscore;
		if (pvalueDiv2 < Double.MIN_NORMAL){
			zscore = LARGEST_ZSCORE;	
		} else {
			zscore = Probability.normalInverse(pvalueDiv2);
		}
		// Min / plus might look counter intuative but i omit 1 - p/2 above so here I have to swap
		if(a1Counts.getQuick(i) < a2Counts.getQuick(i)){
			zscoreSum -= zscore;
		} else {
			zscoreSum += zscore;
		}
	}

	countPearsonR = regression.getR();
	metaZscore = zscoreSum / Math.sqrt(a1Counts.size());
	metaPvalue = 2 * Probability.normal(-Math.abs(metaZscore));
	mle = new AseMleBeta(a1Counts, a2Counts);

}
 

private static DoubleMatrixDataset<String, String> createLocalGeneCorrelation(final DoubleMatrixDataset<String, String> geneZscoresNullGwasCorrelationSubset, final ArrayList<Gene> genes, final int correlationWindow) {

		if (genes.size() != geneZscoresNullGwasCorrelationSubset.rows()) {
			throw new RuntimeException("Genes should match geneZscoresNullGwasCorrelationSubset");
		}

		final DoubleMatrixDataset<String, String> correlations = new DoubleMatrixDataset<>(geneZscoresNullGwasCorrelationSubset.getHashRows(), geneZscoresNullGwasCorrelationSubset.getHashRows());
		final DoubleMatrix2D correlationMatrix = correlations.getMatrix();
		final int geneCount = geneZscoresNullGwasCorrelationSubset.rows();
		final int nullGwasCount = geneZscoresNullGwasCorrelationSubset.columns();
		DoubleMatrix2D geneZscoresNullGwasCorrelationSubsetMatrix = geneZscoresNullGwasCorrelationSubset.getMatrix();

		final SimpleRegression regression = new SimpleRegression();

		for (int i = geneCount; --i >= 0;) {
			for (int j = i + 1; --j >= 0;) {
				regression.clear();

				if (i == j) {
					correlationMatrix.setQuick(i, j, 1);
				} else {

					//Genes should be in the same order as the matrix
					Gene geneI = genes.get(i);
					Gene geneJ = genes.get(j);

					//Only look at position because this is done per chromosome arm
					int geneIStart = geneI.getStart();
					int geneIStop = geneI.getStop();

					int geneJStart = geneJ.getStart();
					int geneJStop = geneJ.getStop();

					if (Math.abs(geneIStart - geneJStart) <= correlationWindow
							|| Math.abs(geneIStart - geneJStop) <= correlationWindow
							|| Math.abs(geneIStop - geneJStart) <= correlationWindow
							|| Math.abs(geneIStop - geneJStop) <= correlationWindow) {
						for (int n = 0; n < nullGwasCount; ++n) {
							regression.addData(geneZscoresNullGwasCorrelationSubsetMatrix.getQuick(i, n), geneZscoresNullGwasCorrelationSubsetMatrix.getQuick(j, n));
						}

						double x = regression.getR();

						correlationMatrix.setQuick(i, j, x);
						correlationMatrix.setQuick(j, i, x); // symmetric
					}

				}
			}
		}

		return correlations;

	}
 

private static DoubleMatrixDataset<String, String> createLocalGeneCorrelation(final DoubleMatrixDataset<String, String> geneZscoresNullGwasCorrelationSubset, final ArrayList<Gene> genes, final int correlationWindow) {

		final DoubleMatrixDataset<String, String> correlations = new DoubleMatrixDataset<>(geneZscoresNullGwasCorrelationSubset.getHashRows(), geneZscoresNullGwasCorrelationSubset.getHashRows());
		final DoubleMatrix2D correlationMatrix = correlations.getMatrix();
		final int geneCount = geneZscoresNullGwasCorrelationSubset.rows();
		final int nullGwasCount = geneZscoresNullGwasCorrelationSubset.columns();
		DoubleMatrix2D geneZscoresNullGwasCorrelationSubsetMatrix = geneZscoresNullGwasCorrelationSubset.getMatrix();

		final SimpleRegression regression = new SimpleRegression();

		for (int i = geneCount; --i >= 0;) {
			for (int j = i + 1; --j >= 0;) {
				regression.clear();

				if (i == j) {
					correlationMatrix.setQuick(i, j, 1);
				} else {

					//Genes should be in the same order as the matrix
					Gene geneI = genes.get(i);
					Gene geneJ = genes.get(j);

					//Only look at position because this is done per chromosome arm
					int geneIStart = geneI.getStart();
					int geneIStop = geneI.getStop();

					int geneJStart = geneJ.getStart();
					int geneJStop = geneJ.getStop();

					if (Math.abs(geneIStart - geneJStart) <= correlationWindow
							|| Math.abs(geneIStart - geneJStop) <= correlationWindow
							|| Math.abs(geneIStop - geneJStart) <= correlationWindow
							|| Math.abs(geneIStop - geneJStop) <= correlationWindow) {
						for (int n = 0; n < nullGwasCount; ++n) {
							regression.addData(geneZscoresNullGwasCorrelationSubsetMatrix.getQuick(i, n), geneZscoresNullGwasCorrelationSubsetMatrix.getQuick(j, n));
						}

						double x = regression.getR();

						correlationMatrix.setQuick(i, j, x);
						correlationMatrix.setQuick(j, i, x); // symmetric
					}

				}
			}
		}

		return correlations;

	}
 

/**
 * @param args the command line arguments
 */
public static void main(String[] args) throws IOException, Exception {

	final File pathwayMatrixFile = new File(args[0]);
	final File significantTermsFile = new File(args[1]);
	final File sumChi2MatrixFile = new File(args[2]);
	final File transQtlEnrichmentsMatrixFile = new File(args[3]);
	
	System.out.println("Pathway file: " + pathwayMatrixFile.getPath());
	System.out.println("Pathway significant terms file: " + significantTermsFile.getPath());
	System.out.println("SumChi2 file: " + sumChi2MatrixFile.getPath());
	System.out.println("Output file: " + transQtlEnrichmentsMatrixFile.getPath());
	
	LinkedHashSet<String> significantTerms = loadSignificantTerms(significantTermsFile);
	
	DoubleMatrixDataset<String, String> pathwayMatrix = DoubleMatrixDataset.loadDoubleData(pathwayMatrixFile.getPath());
	DoubleMatrixDataset<String, String> sumChi2Matrix = DoubleMatrixDataset.loadDoubleData(sumChi2MatrixFile.getPath());

	LinkedHashSet<String> genesInBoth = new LinkedHashSet<String>();

	for (String gene : pathwayMatrix.getHashRows().keySet()) {
		if (sumChi2Matrix.containsRow(gene)) {
			genesInBoth.add(gene);
		}
	}

	pathwayMatrix = pathwayMatrix.viewColSelection(significantTerms);
	
	pathwayMatrix = pathwayMatrix.viewRowSelection(genesInBoth);
	DoubleMatrixDataset<String, String> transQtlEnrichmentsMatrix = new DoubleMatrixDataset<String, String>(pathwayMatrix.getHashCols(), sumChi2Matrix.getHashCols());
	sumChi2Matrix = sumChi2Matrix.viewRowSelection(genesInBoth);

	System.out.println("Genes in both datasets: " + genesInBoth.size());

	System.out.println("Pathways to test: " + pathwayMatrix.columns());
	
	final SimpleRegression regression = new SimpleRegression();
	final DoubleRandomEngine randomEngine = new DRand();
	StudentT tDistColt = new StudentT(sumChi2Matrix.rows() / 2 - 2, randomEngine);
	
	for (String trait : sumChi2Matrix.getColObjects()) {
		
		System.out.println("Trait: " + trait);

		DoubleMatrix1D traitSumChi2 = sumChi2Matrix.getCol(trait);

		for (String pathway : pathwayMatrix.getColObjects()) {

			DoubleMatrix1D pathwayScores = pathwayMatrix.getCol(pathway);

			regression.clear();

			for (int i = 0; i < traitSumChi2.size(); ++i) {

				//System.out.println(traitSumChi2.get(i) + " & " + pathwayScores.get(i));
				
				regression.addData(traitSumChi2.get(i), pathwayScores.get(i));

			}

			double r = regression.getR();

			//System.out.println(trait + " " + pathway + " " + r);
			
			double t = r / (Math.sqrt((1 - r * r) / (double) (traitSumChi2.size() / 2 - 2)));
			double pValue;
			double zScore;
			if (t < 0) {
				pValue = tDistColt.cdf(t);
				if (pValue < 2.0E-323) {
					pValue = 2.0E-323;
				}
				zScore = Probability.normalInverse(pValue);
			} else {
				pValue = tDistColt.cdf(-t);
				if (pValue < 2.0E-323) {
					pValue = 2.0E-323;
				}
				zScore = -Probability.normalInverse(pValue);
			}
			pValue *= 2;

			transQtlEnrichmentsMatrix.setElement(pathway, trait, zScore);

		}

	}
	
	transQtlEnrichmentsMatrix.save(transQtlEnrichmentsMatrixFile);

}