HowTo/PGLS - Revision history

Hilmar: /* Fitting a Brownian Motion model in PGLS */

2008-02-06T19:21:28Z

Fitting a Brownian Motion model in PGLS

← Older revision		Revision as of 19:21, 6 February 2008
Line 17:		Line 17:
	geospiza.data<-read.table("geospiza.txt",row.names=1)		geospiza.data<-read.table("geospiza.txt",row.names=1)

	where the row names are the taxon names. NOTE: you must associate the taxon names with the trait values so that the values can be correctly tied to the tips of the tree! See also [~~https://www.nescent.org/wg_phyloinformatics/~~R_Hackathon/Ancestral_State_Reconstruction ~~this page~~]. With the dataset in active memory, we can create the dataframe as follows:		where the row names are the taxon names. NOTE: you must associate the taxon names with the trait values so that the values can be correctly tied to the tips of the tree! See also [[R_Hackathon 1/Ancestral_State_Reconstruction\|Ancestral State Reconstruction]]. With the dataset in active memory, we can create the dataframe as follows:

	wingL<-geospiza.data$wingL		wingL<-geospiza.data$wingL

Hilmar at 19:20, 6 February 2008

2008-02-06T19:20:54Z

← Older revision		Revision as of 19:20, 6 February 2008
Line 1:		Line 1:
	~~'''~~Phylogenetic Generalized Least Squares~~'''~~		==Phylogenetic Generalized Least Squares==

	PGLS is a powerful method for analyzing continuous data that has been applied to estimating adaptive optima (Butler and King 2004) and estimating the relationships among traits (e.g., body size and geographic range size in carnivores). PGLS allows the user to specify different ways in which the tree structure is expected to affect the covariance in trait values across taxa. For example, the user might assume that the trait evolves by Brownian motion and thus that the trait covariance between any pair of taxa decreases linearly with the time (in branch length) since their divergence. Alternately, the user might apply a Ornstein-Uhlenbeck model where the expected covariance decreases exponentially, as governed by the parameter alpha (Martins and Hansen 1997). These methods are implemented in the ape package.		PGLS is a powerful method for analyzing continuous data that has been applied to estimating adaptive optima (Butler and King 2004) and estimating the relationships among traits (e.g., body size and geographic range size in carnivores). PGLS allows the user to specify different ways in which the tree structure is expected to affect the covariance in trait values across taxa. For example, the user might assume that the trait evolves by Brownian motion and thus that the trait covariance between any pair of taxa decreases linearly with the time (in branch length) since their divergence. Alternately, the user might apply a Ornstein-Uhlenbeck model where the expected covariance decreases exponentially, as governed by the parameter alpha (Martins and Hansen 1997). These methods are implemented in the ape package.

	~~'''~~Fitting a Brownian Motion model in PGLS~~'''~~		===Fitting a Brownian Motion model in PGLS===

	Let's return to the ''Geospiza'' dataset (within the geiger package) to try PGLS. We assume that you have already loaded the necessary packages (geiger for the data and ape for the function) ~~as described on~~ [~~https://www.nescent.org/wg_phyloinformatics/~~R_Hackathon/TransitionProbability ~~this page~~]. Let's say we want to test whether there is a significant relationship between wing length and tarsus length, accounting for possible dependence among the data points (trait values) due to phylogenetic relatedness.		Let's return to the ''Geospiza'' dataset (within the geiger package) to try PGLS. We assume that you have already loaded the necessary packages (geiger for the data and ape for the function) [[R_Hackathon/TransitionProbability\|as described]]. Let's say we want to test whether there is a significant relationship between wing length and tarsus length, accounting for possible dependence among the data points (trait values) due to phylogenetic relatedness.

	First, let's prepare our tree for the analysis. Our original tree has 14 taxa but we only have data for 13, so we must prune the tree. We can take out the species for which we have no data (olivacea) as follows:		First, let's prepare our tree for the analysis. Our original tree has 14 taxa but we only have data for 13, so we must prune the tree. We can take out the species for which we have no data (olivacea) as follows:
Line 32:		Line 32:
	The summary shows you that the log likelihood of this model is 9.6 and its AIC score is -13.21. We also see that the effect of tarsus length on wing length is 1.37, and their relationship is signficant (p=0.001).		The summary shows you that the log likelihood of this model is 9.6 and its AIC score is -13.21. We also see that the effect of tarsus length on wing length is 1.37, and their relationship is signficant (p=0.001).

	~~'''~~Fitting an Ornstein-Uhlenbeck Motion model in PGLS~~'''~~		===Fitting an Ornstein-Uhlenbeck Motion model in PGLS===

	Again we will first build the correlation structure, this time assuming if the traits evolve as expected under the Ornstein-Uhlenbeck process with a variance-restraining parameter, alpha. Ape automatically estimates the best fitting value of alpha for your data.		Again we will first build the correlation structure, this time assuming if the traits evolve as expected under the Ornstein-Uhlenbeck process with a variance-restraining parameter, alpha. Ape automatically estimates the best fitting value of alpha for your data.
Line 41:		Line 41:

	The summary shows you that the estimated alpha is 8.13, the log likelihood of this model is 10.73, and its AIC score is -13.47. This model estimates that the effect of tarsus length on wing length is 0.73, and their relationship is still significant (p=0.01).		The summary shows you that the estimated alpha is 8.13, the log likelihood of this model is 10.73, and its AIC score is -13.47. This model estimates that the effect of tarsus length on wing length is 0.73, and their relationship is still significant (p=0.01).

			[[Category:HowTo]]
			[[Category:Comparative Methods Help]]

Hilmar: R Hackathon/PGLS moved to HowTo/PGLS

2008-02-06T19:18:55Z

R Hackathon/PGLS moved to HowTo/PGLS

New page

'''Phylogenetic Generalized Least Squares'''

PGLS is a powerful method for analyzing continuous data that has been applied to estimating adaptive optima (Butler and King 2004) and estimating the relationships among traits (e.g., body size and geographic range size in carnivores). PGLS allows the user to specify different ways in which the tree structure is expected to affect the covariance in trait values across taxa. For example, the user might assume that the trait evolves by Brownian motion and thus that the trait covariance between any pair of taxa decreases linearly with the time (in branch length) since their divergence. Alternately, the user might apply a Ornstein-Uhlenbeck model where the expected covariance decreases exponentially, as governed by the parameter alpha (Martins and Hansen 1997). These methods are implemented in the ape package.

'''Fitting a Brownian Motion model in PGLS'''

Let's return to the ''Geospiza'' dataset (within the geiger package) to try PGLS. We assume that you have already loaded the necessary packages (geiger for the data and ape for the function) as described on [https://www.nescent.org/wg_phyloinformatics/R_Hackathon/TransitionProbability this page]. Let's say we want to test whether there is a significant relationship between wing length and tarsus length, accounting for possible dependence among the data points (trait values) due to phylogenetic relatedness.

First, let's prepare our tree for the analysis. Our original tree has 14 taxa but we only have data for 13, so we must prune the tree. We can take out the species for which we have no data (olivacea) as follows:

data(geospiza)
attach(geospiza)
geospiza13.tree<-drop.tip(geospiza.tree,"olivacea")

Now let's create a data frame containing just our traits of interest, with the row.names matching the tip.labels. Most users will bring their data in from a tab-delimited table, for example by typing:

geospiza.data<-read.table("geospiza.txt",row.names=1)

where the row names are the taxon names. NOTE: you must associate the taxon names with the trait values so that the values can be correctly tied to the tips of the tree! See also [https://www.nescent.org/wg_phyloinformatics/R_Hackathon/Ancestral_State_Reconstruction this page]. With the dataset in active memory, we can create the dataframe as follows:

wingL<-geospiza.data$wingL
tarsusL<-geospiza.data$tarsusL
DF.geospiza<-data.frame(wingL,tarsusL,row.names=geospiza13.tree$tip.label)
DF.geospiza

When you type the final command, you will see that you have a data frame where the row names are the tip.labels in the pruned tree. Now we will first build the correlation structure expected if the traits evolve by Brownian motion and fit a generalized least squares model assuming this correlation structure.

bm.geospiza<-corBrownian(phy=geospiza13.tree)
bm.gls<-gls(wingL~tarsusL,correlation=bm.geospiza,data=DF.geospiza)
summary(bm.gls)

The summary shows you that the log likelihood of this model is 9.6 and its AIC score is -13.21. We also see that the effect of tarsus length on wing length is 1.37, and their relationship is signficant (p=0.001).

'''Fitting an Ornstein-Uhlenbeck Motion model in PGLS'''

Again we will first build the correlation structure, this time assuming if the traits evolve as expected under the Ornstein-Uhlenbeck process with a variance-restraining parameter, alpha. Ape automatically estimates the best fitting value of alpha for your data.

ou.geospiza<-corMartins(1,phy=geospiza13.tree)
ou.gls<-gls(wingL~tarsusL,correlation=bm.geospiza,data=DF.geospiza)
summary(ou.gls)

The summary shows you that the estimated alpha is 8.13, the log likelihood of this model is 10.73, and its AIC score is -13.47. This model estimates that the effect of tarsus length on wing length is 0.73, and their relationship is still significant (p=0.01).