Nothing
#### Examples-code.R original: 30-June-2000;
#### a few later (August, October, November, March-2001) modifications.
### This file contains all the commands used to generate the examples
### shown in Díaz-Uriarte & Garland, in review.
### All the files that are read are in the 'Examples' directory
### of the package (which is also where this file you are reading
### was originally located).
### For reasons of space, the simulated data files (*.sim) are small
### files with only 100 simulated data sets (in contrast to the 5,000
### simulations used in the paper). To obtain more realistic
### results, you might want to run PDSIMUL and generate a larger
### number of simulations. (For completeness: the seeds we used
### were 1000 for 49lbr.inp, 100 for 49ms.pdi and 20 for 49hmt.pdi).
### To reproduce the results of the paper, the simplest would be
### to either copy all these files to a new directory and launch
### R from there, or launch R from the 'Examples' directory (this
### might not work under Linux if you don't have writing permissions
### to that directory and want to generate postscript files without
### specifying another directory).
### You can run this file either one command at a time (cutting and
### pasting or, better yet, using ESS), or you can directly source the
### file into R.
### To use "source" do:
### > source("Examples-code.R",echo=TRUE,print.eval=TRUE,max.deparse.length=500)
### that will run the file, echo the commands, and return output.
options(digits=4) # number of output digits
library(PHYLOGR) # load our package
###########################
###########################
## This reads in the simulated data file, and adds a column with the
## herbivore vs. carnivore indicator
ANCOVA.sim <- read.sim.data("49lbr.sim", "49lbr.inp",
variable.names = c("body.mass","home.range"),
other.variables =
data.frame(clade=c(rep("Carnivore",19),
rep("Herbivore",30))))
ancova.fit <- lm.phylog(home.range ~ body.mass + clade, data = ANCOVA.sim) # fit the model
summary(ancova.fit)
#############################
##############################
garland.janis.ic <- cbind(read.table("49ms.fic")[,c(3,4)],
read.table("49hmt.fic")[,c(3,4)]) # read fic files
branch.lengths <- read.table("49ms.fic")[,5] # the branch lengths
garland.janis.ic <- garland.janis.ic/branch.lengths # standardize contrasts
names(garland.janis.ic) <- c("body.mass", "running.speed",
"hind.l.length", "mtf.ratio") # name the variables
garland.janis.ic[garland.janis.ic$body.mass<0,] <- -1 * garland.janis.ic[ garland.janis.ic$body.mass<0, ] # to positivize contrasts
garland.janis.ic$branch.lengths <- branch.lengths # add branch lengths to the file, for completennes
garland.janis.ic$names.contr <- as.factor(read.table("49ms.fic")[,1])
garland.janis.ic$clade.contr <- as.factor( c("root",rep("Carnivore",18), rep("Herbivore",29))) # create indicator for clade
## FIXME: avoid using attach. Use with, etc.
attach(garland.janis.ic)
# correlation of residuals from regressions on body mass
resid.hll <- resid(lm(hind.l.length ~ body.mass - 1)) # this model formula uses a regression through the origin (the -1).
resid.speed <- resid(lm(running.speed ~ body.mass - 1))
rho1 <- cor.origin(resid.hll,resid.speed)
rho1 # the correlation through the origin
## the p-value, using the standard formula but providing the correct df's.
2*(1-pt(sqrt(46)*abs(rho1)/sqrt(1-rho1^2), 46))
# the same, but using multiple regression
plot(body.mass,running.speed)
plot(body.mass,hind.l.length)
plot(hind.l.length,running.speed)
fit1 <- lm(running.speed ~ body.mass + hind.l.length - 1) # note also regression through the origin
plot(fit1)
summary(fit1)
# Excluding the polar bear-grizzly contrast
# using residuals:
resid.hll.nbear <- resid(lm(hind.l.length ~ body.mass - 1, subset = names.contr!="Tm-Ur"))
resid.speed.nbear <- resid(lm(running.speed ~ body.mass - 1, subset=names.contr!="Tm-Ur"))
rho2 <- cor.origin(resid.hll.nbear,resid.speed.nbear)
rho2 # the correlation through the origin
2*(1-pt(sqrt(45)*abs(rho2)/sqrt(1-rho2^2),45)) # the p-value, using the standard formula but providing the correct df's.
# with multiple regression
fit2 <- lm(running.speed ~ body.mass + hind.l.length - 1, subset=names.contr!="Tm-Ur")
plot(fit2)
summary(fit2)
# is MT/F relevant after including hindl.lenght? remember: need to look only at coefficient for mtf.
summary(lm(running.speed ~ body.mass + hind.l.length + mtf.ratio - 1, subset=names.contr!="Tm-Ur"))
anova(lm(running.speed ~ body.mass + hind.l.length + mtf.ratio - 1, subset=names.contr!="Tm-Ur")) # anova summary
GarlandJanis.sim <- read.sim.data(c("49ms.sim","49hmt.sim"), pdi.files = c("49ms.pdi","49hmt.pdi"), variable.names = c("body.mass","running.speed","hind.l.length","mtf.ratio"), other.variables = data.frame(clade = c(rep("Carnivore",19),rep("Herbivore",30))))
plot(lm(running.speed ~ body.mass + hind.l.length, data = GarlandJanis.sim, subset= sim.counter==0)) # regression residual plots
speed.hll <- lm.phylog(running.speed ~ body.mass + hind.l.length, data = GarlandJanis.sim)
summary(speed.hll)
#postscript(file="Figure1.ps") # to generate a postscript file; you won't see the figure in the screen.
par(mfrow=c(2,3)) # 2x3 figures in the page, so that all fit in one indow
par(cex.lab=1.5) # to make labels larger
plot(speed.hll) # plot the analyses from lm.phylog
#dev.off()
speed.hll.nbear <- lm.phylog(running.speed ~ body.mass + hind.l.length, data=GarlandJanis.sim, exclude.tips=c("Tm","Ur"))
summary(speed.hll.nbear)
speed.hll.nbear.int <- lm.phylog(running.speed ~ body.mass + hind.l.length*clade, data=GarlandJanis.sim, exclude.tips=c("Tm","Ur"))
summary(speed.hll.nbear.int)
attach(GarlandJanis.sim) # to allow subseting of Tips for excluding either Carnivores or Herbivores
speed.mass.length.carn <- lm.phylog(running.speed ~ body.mass + hind.l.length, data=GarlandJanis.sim, exclude.tips= c(as.character(Tips[clade=="Herbivore" & sim.counter==0]), "Tm","Ur")) # only carnivores, except polar bear and grizzly
summary(speed.mass.length.carn)
# Testing whether MT/F ratio adds anything
speed.mtf <- lm.phylog(running.speed ~ body.mass + hind.l.length + mtf.ratio, data=GarlandJanis.sim, exclude.tips=c("Tm","Ur"))
summary(speed.mtf) # we are only interested here in the F value of mtf.ratio
############### PCA ###############
library(PHYLOGR) # load the package, if you haven't already done so.
LacertidSim <- read.sim.data(c("ifsmi.sim", "ihshw.sim", "iclag.sim", "icfxx.sim"), pdi.files=c("ifsmi.pdi", "ihshw.pdi", "iclag.pdi", "icfxx.pdi"), variable.names= c("svl","svl.matur", "hatsvl", "hatweight", "clutch.size", "age.mat","cl.freq", "xx"))
LacertidSim <- LacertidSim[,-10] # exclude last column
LacertidPCA <- prcomp.phylog(LacertidSim)
summary(LacertidPCA)
#par(mfrow=c(2,4))
#par(cex.lab=1.2)
#plot(LacertidPCA) # shown in Figure 2
## or we can use:
#postscript(file="Figure2.ps")
par(mfrow=c(2,4))
par(cex.lab=1.2)
plot(LacertidPCA)
#dev.off()
LacertidIC <- cbind(read.table("ifsmi.fic")[,c(3,4)], read.table("ihshw.fic")[,c(3,4)], read.table("iclag.fic")[,c(3,4)], read.table("icfxx.fic")[,3]) # the (unstandardized) contrasts
stand <- read.table("ifsmi.fic")[,5] # the square root of the sum of branch lengths
LacertidIC <- LacertidIC/stand # compute the standardized contrasts
LacertidIC$contr <- read.table("ifsmi.fic")[,1] # name of the contrast
names(LacertidIC) <- c("svl","svl.matur", "hatsvl", "hatweight", "clutch.size", "age.mat","cl.freq","ICcontr")
cor.for.ic.pca <- matrix(nrow=7,ncol=7)
for (i in 1:7) for (j in 1:7) cor.for.ic.pca[i,j] <- cor.origin(LacertidIC[[i]],LacertidIC[[j]]) # to obtain the correlation matrix based on regressions through the origin
ic.pca <- svd(cor.for.ic.pca) #so that the eigenvectors are normalized
cor.with.factors <- t(sqrt(ic.pca$d) * t(ic.pca$u))
cor.with.factors <- as.data.frame(cor.with.factors)
ic.pca$d # these are the eigenvalues
100*ic.pca$d/length(ic.pca$d) # percentage of variance
names(cor.with.factors) <- paste("PC",seq(1:7),sep="") # nicer names for output
row.names(cor.with.factors)<-names(LacertidIC)[-8]
cor.with.factors # the correlation between the variables and the components
######################################
############# ####################
############# GLS ####################
############# ####################
######################################
#### Lacertids
LacertidOriginal <- LacertidSim[LacertidSim$sim.counter==0,-1] # we could read the pdi data again, but this is another simple way, since the original data are in the file with the simulated data
Lacertid.varcov <- read.phylog.matrix("ifsmi.dsc") # from pddist, option 5, in matrix form, with header
# attach the data set, so variables can be accessed directly
# without reference to the data set
detach(GarlandJanis.sim)
attach(LacertidOriginal)
t(apply(LacertidOriginal[,-c(1,2)],2,function(y){summary(phylog.gls.fit(svl,y,Lacertid.varcov))[[4]][2,c(1,2)]}))
#### Garland & Janis
Garland.Janis.Orig <- read.pdi.data(c("49ms.pdi","49hmt.pdi"), variable.names = c("body.mass", "running.speed", "hind.l.length", "mtf.ratio"))
Garland.Janis.Orig$clade <- as.factor(c(rep("Carnivore",19), rep("Herbivore",30)))
Garland.Janis.Cov <- read.phylog.matrix("49ms.dsc")
detach(LacertidOriginal)
detach(garland.janis.ic)
attach(Garland.Janis.Orig)
# all data
fit.gls.gj <- phylog.gls.fit(cbind(body.mass,hind.l.length), running.speed, Garland.Janis.Cov) # the model fitting call
summary(fit.gls.gj) # summary of the gls model; same as with IC
# without the grizzly and polar bear
summary(phylog.gls.fit(cbind(body.mass,hind.l.length), running.speed, Garland.Janis.Cov, exclude.tips=c("Tm","Ur"))) # as expected, not identically equal to IC
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.