Example 16 with Regression
use of dr.stats.Regression in project beast-mcmc by beast-dev.
the class TempestFrame method writeTimeTreeFile.
protected void writeTimeTreeFile(PrintStream ps) throws IOException {
FlexibleTree tree = new FlexibleTree(tempestPanel.getTreeAsViewed());
Regression r = tempestPanel.getTemporalRooting().getRootToTipRegression(tempestPanel.getTreeAsViewed());
for (int i = 0; i < tree.getInternalNodeCount(); i++) {
NodeRef node = tree.getInternalNode(i);
double height = tree.getNodeHeight(node);
tree.setNodeHeight(node, height / r.getGradient());
}
TreeUtils.setHeightsFromDates(tree);
NexusExporter nexusExporter = new NexusExporter(new PrintStream(ps));
nexusExporter.exportTree(tree);
}
Also used :
NodeRef(dr.evolution.tree.NodeRef)
NexusExporter(dr.app.tools.NexusExporter)
FlexibleTree(dr.evolution.tree.FlexibleTree)
Regression(dr.stats.Regression)
Example 17 with Regression
use of dr.stats.Regression in project beast-mcmc by beast-dev.
the class ContinuousDiffusionStatistic method getStatisticValue.
public double getStatisticValue(int dim) {
double treeLength = 0;
double treeDistance = 0;
double totalMaxDistanceFromRoot = 0;
// can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
double maxDistanceFromRootCumulative = 0;
double maxBranchDistanceFromRoot = 0;
// can only be used when cumulative and not associated with discrete state (not based on the distances on the branches from the root up that point)
double maxDistanceOverTimeFromRootWA = 0;
double maxBranchDistanceOverTimeFromRootWA = 0;
List<Double> rates = new ArrayList<Double>();
List<Double> distances = new ArrayList<Double>();
List<Double> times = new ArrayList<Double>();
List<Double> traits = new ArrayList<Double>();
List<double[]> traits2D = new ArrayList<double[]>();
//double[] diffusionCoefficients = null;
List<Double> diffusionCoefficients = new ArrayList<Double>();
double waDiffusionCoefficient = 0;
double lowerHeight = heightLowers[dim];
double upperHeight = Double.MAX_VALUE;
if (heightLowers.length == 1) {
upperHeight = heightUpper;
} else {
if (dim > 0) {
if (!cumulative) {
upperHeight = heightLowers[dim - 1];
}
}
}
for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
MultivariateTraitTree tree = traitLikelihood.getTreeModel();
BranchRateModel branchRates = traitLikelihood.getBranchRateModel();
String traitName = traitLikelihood.getTraitName();
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef node = tree.getNode(i);
if (node != tree.getRoot()) {
NodeRef parentNode = tree.getParent(node);
boolean testNode = true;
if (branchset.equals(BranchSet.CLADE)) {
try {
testNode = inClade(tree, node, taxonList);
} catch (TreeUtils.MissingTaxonException mte) {
throw new RuntimeException(mte.toString());
}
} else if (branchset.equals(BranchSet.BACKBONE)) {
if (backboneTime > 0) {
testNode = onAncestralPathTime(tree, node, backboneTime);
} else {
try {
testNode = onAncestralPathTaxa(tree, node, taxonList);
} catch (TreeUtils.MissingTaxonException mte) {
throw new RuntimeException(mte.toString());
}
}
}
if (testNode) {
if ((tree.getNodeHeight(parentNode) > lowerHeight) && (tree.getNodeHeight(node) < upperHeight)) {
double[] trait = traitLikelihood.getTraitForNode(tree, node, traitName);
double[] parentTrait = traitLikelihood.getTraitForNode(tree, parentNode, traitName);
double[] traitUp = parentTrait;
double[] traitLow = trait;
double timeUp = tree.getNodeHeight(parentNode);
double timeLow = tree.getNodeHeight(node);
double rate = (branchRates != null ? branchRates.getBranchRate(tree, node) : 1.0);
// System.out.println(rate);
MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
double[] precision = diffModel.getPrecisionParameter().getParameterValues();
History history = null;
if (stateString != null) {
history = setUpHistory(markovJumpLikelihood.getHistoryForNode(tree, node, SITE), markovJumpLikelihood.getStatesForNode(tree, node)[SITE], markovJumpLikelihood.getStatesForNode(tree, parentNode)[SITE], timeLow, timeUp);
}
if (tree.getNodeHeight(parentNode) > upperHeight) {
timeUp = upperHeight;
traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
if (stateString != null) {
history.truncateUpper(timeUp);
}
}
if (tree.getNodeHeight(node) < lowerHeight) {
timeLow = lowerHeight;
traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, trueNoise);
if (stateString != null) {
history.truncateLower(timeLow);
}
}
if (dimension > traitLow.length) {
System.err.println("specified trait dimension for continuous trait summary, " + dimension + ", is > dimensionality of trait, " + traitLow.length + ". No trait summarized.");
} else {
traits.add(traitLow[(dimension - 1)]);
}
if (traitLow.length == 2) {
traits2D.add(traitLow);
}
double time;
if (stateString != null) {
time = history.getStateTime(stateString);
// System.out.println("tine before = "+(timeUp - timeLow)+", time after= "+time);
} else {
time = timeUp - timeLow;
}
treeLength += time;
times.add(time);
//setting up continuous trait values for heights in discrete trait history
if (stateString != null) {
history.setTraitsforHeights(traitUp, traitLow, precision, rate, trueNoise);
}
double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
double timeFromRoot = (tree.getNodeHeight(tree.getRoot()) - timeLow);
if (useGreatCircleDistances && (trait.length == 2)) {
// Great Circle distance
double distance;
if (stateString != null) {
distance = history.getStateGreatCircleDistance(stateString);
} else {
distance = getGreatCircleDistance(traitLow, traitUp);
}
distances.add(distance);
if (time > 0) {
treeDistance += distance;
double dc = Math.pow(distance, 2) / (4 * time);
diffusionCoefficients.add(dc);
waDiffusionCoefficient += (dc * time);
rates.add(distance / time);
}
SphericalPolarCoordinates rootCoord = new SphericalPolarCoordinates(rootTrait[0], rootTrait[1]);
double tempDistanceFromRootLow = rootCoord.distance(new SphericalPolarCoordinates(traitUp[0], traitUp[1]));
if (tempDistanceFromRootLow > totalMaxDistanceFromRoot) {
totalMaxDistanceFromRoot = tempDistanceFromRootLow;
if (stateString != null) {
double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, true);
if (stateTimeDistance[0] > 0) {
maxDistanceFromRootCumulative = tempDistanceFromRootLow * (stateTimeDistance[0] / timeFromRoot);
maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
maxBranchDistanceFromRoot = stateTimeDistance[1];
maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
}
} else {
maxDistanceFromRootCumulative = tempDistanceFromRootLow;
maxDistanceOverTimeFromRootWA = tempDistanceFromRootLow / timeFromRoot;
double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, true);
maxBranchDistanceFromRoot = timeDistance[1];
maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
}
//distance between traitLow and traitUp for maxDistanceFromRootCumulative
if (timeUp == upperHeight) {
if (time > 0) {
maxDistanceFromRootCumulative = distance;
maxDistanceOverTimeFromRootWA = distance / time;
maxBranchDistanceFromRoot = distance;
maxBranchDistanceOverTimeFromRootWA = distance / time;
}
}
}
} else {
double distance;
if (stateString != null) {
distance = history.getStateNativeDistance(stateString);
} else {
distance = getNativeDistance(traitLow, traitUp);
}
distances.add(distance);
if (time > 0) {
treeDistance += distance;
double dc = Math.pow(distance, 2) / (4 * time);
diffusionCoefficients.add(dc);
waDiffusionCoefficient += dc * time;
rates.add(distance / time);
}
double tempDistanceFromRoot = getNativeDistance(traitLow, rootTrait);
if (tempDistanceFromRoot > totalMaxDistanceFromRoot) {
totalMaxDistanceFromRoot = tempDistanceFromRoot;
if (stateString != null) {
double[] stateTimeDistance = getStateTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, precision, branchRates, false);
if (stateTimeDistance[0] > 0) {
maxDistanceFromRootCumulative = tempDistanceFromRoot * (stateTimeDistance[0] / timeFromRoot);
maxDistanceOverTimeFromRootWA = maxDistanceFromRootCumulative / stateTimeDistance[0];
maxBranchDistanceFromRoot = stateTimeDistance[1];
maxBranchDistanceOverTimeFromRootWA = stateTimeDistance[1] / stateTimeDistance[0];
}
} else {
maxDistanceFromRootCumulative = tempDistanceFromRoot;
maxDistanceOverTimeFromRootWA = tempDistanceFromRoot / timeFromRoot;
double[] timeDistance = getTimeAndDistanceFromRoot(tree, node, timeLow, traitLikelihood, traitName, traitLow, false);
maxBranchDistanceFromRoot = timeDistance[1];
maxBranchDistanceOverTimeFromRootWA = timeDistance[1] / timeDistance[0];
}
//distance between traitLow and traitUp for maxDistanceFromRootCumulative
if (timeUp == upperHeight) {
if (time > 0) {
maxDistanceFromRootCumulative = distance;
maxDistanceOverTimeFromRootWA = distance / time;
maxBranchDistanceFromRoot = distance;
maxBranchDistanceOverTimeFromRootWA = distance / time;
}
}
}
}
}
}
}
}
}
if (summaryStat == summaryStatistic.DIFFUSION_RATE) {
if (summaryMode == Mode.AVERAGE) {
return DiscreteStatistics.mean(toArray(rates));
} else if (summaryMode == Mode.MEDIAN) {
return DiscreteStatistics.median(toArray(rates));
} else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
final double mean = DiscreteStatistics.mean(toArray(rates));
return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
//weighted average
} else {
return treeDistance / treeLength;
}
} else if (summaryStat == summaryStatistic.TRAIT) {
if (summaryMode == Mode.MEDIAN) {
return DiscreteStatistics.median(toArray(traits));
} else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
// don't compute mean twice
final double mean = DiscreteStatistics.mean(toArray(traits));
return Math.sqrt(DiscreteStatistics.variance(toArray(traits), mean)) / mean;
// default is average. A warning is thrown by the parser when trying to use WEIGHTED_AVERAGE
} else {
return DiscreteStatistics.mean(toArray(traits));
}
} else if (summaryStat == summaryStatistic.TRAIT2DAREA) {
double area = getAreaFrom2Dtraits(traits2D, 0.99);
return area;
} else if (summaryStat == summaryStatistic.DIFFUSION_COEFFICIENT) {
if (summaryMode == Mode.AVERAGE) {
return DiscreteStatistics.mean(toArray(diffusionCoefficients));
} else if (summaryMode == Mode.MEDIAN) {
return DiscreteStatistics.median(toArray(diffusionCoefficients));
} else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
// don't compute mean twice
final double mean = DiscreteStatistics.mean(toArray(diffusionCoefficients));
return Math.sqrt(DiscreteStatistics.variance(toArray(diffusionCoefficients), mean)) / mean;
} else {
return waDiffusionCoefficient / treeLength;
}
//wavefront distance
//TODO: restrict to non state-specific wavefrontDistance/rate
} else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
return maxDistanceFromRootCumulative;
// return maxBranchDistanceFromRoot;
} else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE_PHYLO) {
return maxBranchDistanceFromRoot;
//wavefront rate, only weighted average TODO: extend for average, median, COEFFICIENT_OF_VARIATION?
} else if (summaryStat == summaryStatistic.WAVEFRONT_RATE) {
return maxDistanceOverTimeFromRootWA;
// return maxBranchDistanceOverTimeFromRootWA;
} else if (summaryStat == summaryStatistic.DIFFUSION_DISTANCE) {
return treeDistance;
//DIFFUSION_TIME
} else if (summaryStat == summaryStatistic.DISTANCE_TIME_CORRELATION) {
if (summaryMode == Mode.SPEARMAN) {
return getSpearmanRho(convertDoubles(times), convertDoubles(distances));
} else if (summaryMode == Mode.R_SQUARED) {
Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
return r.getRSquared();
} else {
Regression r = new Regression(convertDoubles(times), convertDoubles(distances));
return r.getCorrelationCoefficient();
}
} else {
return treeLength;
}
}
Also used :
SphericalPolarCoordinates(dr.geo.math.SphericalPolarCoordinates)
Regression(dr.stats.Regression)
MultivariateTraitTree(dr.evolution.tree.MultivariateTraitTree)
NodeRef(dr.evolution.tree.NodeRef)
BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel)
TreeUtils(dr.evolution.tree.TreeUtils)
Aggregations
Regression (dr.stats.Regression)17
NodeRef (dr.evolution.tree.NodeRef)3
NexusExporter (dr.app.tools.NexusExporter)2
FlexibleTree (dr.evolution.tree.FlexibleTree)2
Variate (dr.stats.Variate)2
NumberFormatter (dr.util.NumberFormatter)2
PrintWriter (java.io.PrintWriter)2
Node (jebl.evolution.graphs.Node)2
Taxon (jebl.evolution.taxa.Taxon)2
RootedTree (jebl.evolution.trees.RootedTree)2
MultivariateTraitTree (dr.evolution.tree.MultivariateTraitTree)1
TreeUtils (dr.evolution.tree.TreeUtils)1
BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)1
SphericalPolarCoordinates (dr.geo.math.SphericalPolarCoordinates)1
