From 8bf5e2d0c414150143204b64f3b84104dfdb15a7 Mon Sep 17 00:00:00 2001 From: Denise Date: Wed, 20 Apr 2016 11:38:16 +0200 Subject: [PATCH] When having sampled ancestors, the removal probability r doesn't affect the sampling proportion s. (that is, s=psi/(psi+mu)). --- .../speciation/BirthDeathMigrationModel.java | 28 +++++++++---------- .../BirthDeathMigrationModelUncoloured.java | 13 ++++----- ...ecewiseBirthDeathSamplingDistribution.java | 3 -- 3 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/beast/evolution/speciation/BirthDeathMigrationModel.java b/src/beast/evolution/speciation/BirthDeathMigrationModel.java index ee1f58a..a38765e 100755 --- a/src/beast/evolution/speciation/BirthDeathMigrationModel.java +++ b/src/beast/evolution/speciation/BirthDeathMigrationModel.java @@ -79,7 +79,7 @@ public void initAndValidate() { if (originBranch==null) throw new RuntimeException("Error: Origin specified but originBranch missing!"); - updateOrigin(coltree.getRoot()); + updateOrigin(coltree.getRoot()); if (!Boolean.valueOf(System.getProperty("beast.resume")) && treeInput.get().getRoot().getHeight() >= origin.get().getValue()) throw new RuntimeException("Error: origin("+T+") must be larger than tree height("+coltree.getRoot().getHeight()+")!"); @@ -197,8 +197,8 @@ void setupIntegrators(){ // set up ODE's and integrators state = i/totalIntervals; rho[i]= rhoChanges>0? - rhoSamplingChangeTimes.contains(times[i]) ? rhos[rhos.length > n ? (rhoChanges+1)*state+index(times[i%totalIntervals], rhoSamplingChangeTimes) : state] : 0. - : rhos[0]; + rhoSamplingChangeTimes.contains(times[i]) ? rhos[rhos.length > n ? (rhoChanges+1)*state+index(times[i%totalIntervals], rhoSamplingChangeTimes) : state] : 0. + : rhos[0]; } } @@ -347,10 +347,10 @@ public double calculateTreeLogLikelihood(TreeInterface tree) { collectTimes(T); setRho(); - if (updateRates() < 0 || (times[totalIntervals-1] > T)) { + if (updateRates() < 0 || (times[totalIntervals-1] > T)) { logP = Double.NEGATIVE_INFINITY; return logP; - } + } double[] noSampleExistsProp = new double[n]; @@ -563,20 +563,18 @@ public void transformParameters(){ // // death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i]; // - if (!SAModel || removalAffectsSamplingProportion.get()) - psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] - * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ; - - if (!SAModel) - death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i]; + if (!SAModel) { + psi[i] = p[p.length > n ? (samplingChanges + 1) * state + index(times[i % totalIntervals], samplingRateChangeTimes) : state] + * ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state]; + death[i] = ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state] - psi[i]; + } else { r[i] = removalProbabilities[removalProbabilities.length > n ? (rChanges+1)*state+index(times[i%totalIntervals], rChangeTimes) : state]; - if (!removalAffectsSamplingProportion.get()) - psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] - * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - / (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]); + psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] + * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] + / (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]); death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i]*r[i]; diff --git a/src/beast/evolution/speciation/BirthDeathMigrationModelUncoloured.java b/src/beast/evolution/speciation/BirthDeathMigrationModelUncoloured.java index 119d0eb..32eb52e 100755 --- a/src/beast/evolution/speciation/BirthDeathMigrationModelUncoloured.java +++ b/src/beast/evolution/speciation/BirthDeathMigrationModelUncoloured.java @@ -695,18 +695,17 @@ public void transformParameters(){ birth[i] = R[R.length > n ? (birthChanges+1)*state+index(times[i%totalIntervals], birthRateChangeTimes) : state] * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ; - if (!SAModel || removalAffectsSamplingProportion.get()) - psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] - * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ; + if (!SAModel) { + psi[i] = p[p.length > n ? (samplingChanges + 1) * state + index(times[i % totalIntervals], samplingRateChangeTimes) : state] + * ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state]; - if (!SAModel) - death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i]; + death[i] = ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state] - psi[i]; + } else { r[i] = removalProbabilities[removalProbabilities.length > n ? (rChanges+1)*state+index(times[i%totalIntervals], rChangeTimes) : state]; - if (!removalAffectsSamplingProportion.get()) - psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] + psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state] * ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] / (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]); diff --git a/src/beast/evolution/speciation/PiecewiseBirthDeathSamplingDistribution.java b/src/beast/evolution/speciation/PiecewiseBirthDeathSamplingDistribution.java index 617a924..4fc2cce 100644 --- a/src/beast/evolution/speciation/PiecewiseBirthDeathSamplingDistribution.java +++ b/src/beast/evolution/speciation/PiecewiseBirthDeathSamplingDistribution.java @@ -119,9 +119,6 @@ public abstract class PiecewiseBirthDeathSamplingDistribution extends SpeciesTre public Input stateNumber = new Input<>("stateNumber", "The number of states or locations", Input.Validate.REQUIRED); - public Input removalAffectsSamplingProportion = - new Input<>("removalAffectsSamplingProportion", "In R0 param, is samplingProportion = samplingRate/(r*samplingRate+deathRate)? Default=true. (Alternative: samplingProportion = samplingRate/(samplingRate+deathRate)) ", true); - public Input adjustTimesInput = new Input<>("adjustTimes", "Origin of MASTER sims which has to be deducted from the change time arrays"); //