You are here

RRWs

 

Setting up a continuous phylogeographic analysis using relaxed random walks

This tutorial describes how to modify an XML file that specifies a homogenous Brownian diffusion phylogeographic model to set up relaxed random walks (RRWs, [1]). For the necessary XML additions, we refer to the tutorial describing the standard continuous phylogeographic analysis [ref]. The example XML file for this analysis can be found here (gamma RRW: [2], cauchy RRW: [3], lognormal RRW: [4]) and considers raccoon rabies diffusion in north-eastern United States [5]. RRWs relax the Brownian diffusion assumption that constraints the process to maintain a constant variance over time. They build on uncorrelated relaxed clock models [6]. By assigning to each branch in the tree a rate scalar to rescale the precision matrix, RRWs allow the underlying process to vary from branch to branch in the tree.

In analogy to relaxed clock models, add a discretizedBranchRates block (before the multivariateTraitLikelihood). For the (one parameter) gamma and Cauchy RRW, this should look like this:

	<discretizedBranchRates id="branchRatesDiffusion">
		<treeModel idref="treeModel"/>
		<distribution>
			<onePGammaDistributionModel id="gammaModel">
         		<shape>
        			<parameter id="halfDF" value="0.5"/>
        		</shape>
			</onePGammaDistributionModel>
		</distribution>
		<rateCategories>
			<parameter id="branchRatesDiffusion.categories" dimension="92"/>
		</rateCategories>
	</discretizedBranchRates>

For the lognormal RRW, we require a logNormalDistributionModel with a mean set to 1.0 in the distribution block instead of the onePGammaDistributionModel:

	<discretizedBranchRates id="branchRatesDiffusion">
		<treeModel idref="treeModel"/>
		<distribution>
			<logNormalDistributionModel meanInRealSpace="true">
				<mean>
					<parameter id="ucldDiffusion.mean" value="1.0" lower="0.0" upper="10.0"/>
				</mean>
				<stdev>
					<parameter id="ucldDiffusion.stdev" value="0.1" lower="0.0" upper="10.0"/>
				</stdev>
			</logNormalDistributionModel>
		</distribution>
		<rateCategories>
			<parameter id="branchRatesDiffusion.categories" dimension="92"/>
		</rateCategories>
	</discretizedBranchRates>

Add the associated statistics:

	<rateStatistic id="meanRateDiffusion" name="meanRate" mode="mean" internal="true" external="true">
		<treeModel idref="treeModel"/>
		<discretizedBranchRates idref="branchRatesDiffusion"/>
	</rateStatistic>
	
	<rateStatistic id="coefficientOfVariationDiffusion" name="coefficientOfVariationDiffusion" mode="coefficientOfVariation" internal="true" external="true">
		<treeModel idref="treeModel"/>
		<discretizedBranchRates idref="branchRatesDiffusion"/>
	</rateStatistic>
	
	<rateCovarianceStatistic id="covarianceDiffusion" name="covarianceDiffusion">
		<treeModel idref="treeModel"/>
		<discretizedBranchRates idref="branchRatesDiffusion"/>
	</rateCovarianceStatistic>

Add the standard operators for the categories of the discretized distributions. For the gamma RRW, we require a scaleOperator on the 'halfDF' parameter (degrees of freedom divided by two). The gamma distribution reduces to Cauchy distribution when the 'halfDF' parameter is fixed to 0.5 (the starting value specified for "halfDF"). Therfore we do not require an operator for the Cauchy distribution.

	<operators id="operators" optimizationSchedule="log">
		.
		.
		.
		<scaleOperator scaleFactor="0.75" weight="6">
			<parameter idref="halfDF"/>
		</scaleOperator>
		<swapOperator size="1" weight="30" autoOptimize="false">
			<parameter idref="branchRatesDiffusion.categories"/>
		</swapOperator>
		<randomWalkIntegerOperator windowSize="2" weight="10.0">
			<parameter idref="branchRatesDiffusion.categories"/>
		</randomWalkIntegerOperator>
		<uniformIntegerOperator weight="10.0">
			<parameter idref="branchRatesDiffusion.categories"/>
		</uniformIntegerOperator>
		.
		.
		.
	</operators>

In case of the lognormal RRW, we need a scaleOperator on the standard deviation of the lognormal distribution (the mean remains fixed to 1.0) instead of the "halfDF":

 
	<operators id="operators" optimizationSchedule="log">
		.
		.
		.
 		<scaleOperator scaleFactor="0.75" weight="3">
			<parameter idref="ucldDiffusion.stdev"/>
		</scaleOperator>
		<swapOperator size="1" weight="30" autoOptimize="false">
			<parameter idref="branchRatesDiffusion.categories"/>
		</swapOperator>
		<randomWalkIntegerOperator windowSize="2" weight="10.0">
			<parameter idref="branchRatesDiffusion.categories"/>
		</randomWalkIntegerOperator>
		<uniformIntegerOperator weight="10.0">
			<parameter idref="branchRatesDiffusion.categories"/>
		</uniformIntegerOperator>
		.
		.
		.
 	</operators>

We also specify diffuse priors on these additional parameters. For the gamma RWW, this is a gamma prior with shape=0.001 and scale=1000.0:

		 <gammaPrior offset="0.0" shape="0.001" scale="1000.0">
                 	 <parameter idref="halfDF"/>
 		 </gammaPrior>

For the lognormal RRW, we can use an exponential prior with mean = 10:

		 <exponentialPrior mean="10" offset="0">
 	 	 	 <parameter idref="ucldDiffusion.stdev"/>
	 	 </exponentialPrior>

We can log these parameters to the screen; gamma RRW:

		 <log id="screenLog" logEvery="50000">
			.
			.
			.
		 	 <column label="halfDF" sf="6" width="12">
	 	 	 	 <parameter idref="halfDF"/>
	 	 	 </column>
			.
			.
			.
	 	 </log>

lognormal RRW:

	 	 <log id="screenLog" logEvery="50000">
			.
			.
			.
 	 	 	 <column label="ucldDiffusion" sf="6" width="12">
	 	 	 	 <parameter idref="ucldDiffusion.stdev"/>
	 	 	 </column>
			.
			.
			.
		</log>

And to the logfile; gamma RRW: 

		<log id="fileLog" logEvery="50000" fileName="RacRABV_gammaRRW.log">
			.
			.
			.
			<parameter idref="halfDF"/>
			.
			.
			.
		</log>

lognormal RRW:

		<log id="fileLog" logEvery="50000" fileName="RacRABV_LogNRRW.log">
			.
			.
			.
			<parameter idref="ucldDiffusion.stdev"/>
			.
			.
			.
		</log>