Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## Code for installation of packages inside vignette taken from PCMBase vignette by Venelin Mitov
if(!requireNamespace("ggplot2")) {
message("Building the vignette requires ggplot2 R-package. Trying to install.")
status.ggplot2 <- try({
install.packages("ggplot2")
}, silent = TRUE)
if(class(status.ggplot2 == "try-error")) {
stop(
"The ggplot2 installation did not succeed. The vignette cannot be built.")
}
}
## -----------------------------------------------------------------------------
library(ggplot2)
library(ape)
library(mvSLOUCH)
## -----------------------------------------------------------------------------
load("./Carnivora_mvSLOUCH_objects.RData")
## -----------------------------------------------------------------------------
RNGversion(min(as.character(getRversion()),"3.6.1"))
set.seed(5, kind = "Mersenne-Twister", normal.kind = "Inversion")
## ----eval=TRUE, echo=TRUE, results="hide", message=FALSE----------------------
b_correct_dryad_download<-FALSE
temp <- tempfile()
tryCatch({
download.file("datadryad.org/api/v2/datasets/doi%253A10.5061%252Fdryad.77tm4/download",temp)
b_correct_dryad_download<-TRUE
},error=function(e){cat("Could not download data file from Dryad! No analysis can be done! Vignette not built!")},
warning=function(w){cat("Problem with downloading data file from Dryad! No analysis can be done! Vignette not built!")})
if (b_correct_dryad_download){
b_correct_dryad_download<-FALSE
tryCatch({
dfcarnivores_postcranial <- read.table(unz(temp, "Carnivore Postcranial Morphology Data Samuels et al. 2013.txt"),header=TRUE,sep="\t",stringsAsFactors =FALSE)
b_correct_dryad_download<-TRUE
},error=function(e){cat("Corrupt data file from Dryad! No analysis can be done! Vignette not built!")},
warning=function(w){cat("Problem with accessing data file from Dryad! No analysis can be done! Vignette not built!")})
}
## -----------------------------------------------------------------------------
unlink(temp)
Urocyon_cinereoargenteus_duplicated<-which(dfcarnivores_postcranial[,2]=="Urocyon cinereoargenteus")[2]
dfcarnivores_postcranial<-dfcarnivores_postcranial[-Urocyon_cinereoargenteus_duplicated,]
v_species_to_remove<-c("Bdeogale crassicauda","Lycalopex sp.","Daphoenus vetus","Barbourofelis loveorum","Archaeocyon leptodus","Canis armbrusteri","Canis dirus","Canis latrans orcutti","Canis lupus (Pleistocene)","Desmocyon thomsoni","Hesperocyon gregarius","Mesocyon coryphaeus","Paraenhydrocyon josephi","Phlaocyon leucosteus","Vulpes macrotis (Pleistocene)","Vulpes vulpes (Pleistocene)","Homotherium ischyrus","Homotherium serum","Leopardus wiedii (Pleistocene)","Lynx rufus (Pleistocene)","Miracinonyx inexpectatus","Panthera atrox","Puma concolor (Pleistocene)","Puma lacustris","Smilodon fatalis","Miacis gracilis","Mephitis mephitis (Pleistocene)","Spilogale gracilis (Pleistocene)","Spilogale putorius (Pleistocene)","Gulo gulo (Pleistocene)","Martes americana (Pleistocene)","Martes nobilis (Pleistocene)","Mustela nigripes (Pleistocene)","Neovison frenata (Pleistocene)","Neovison vison (Pleistocene)","Satherium piscinarium","Taxidea taxus (Pleistocene)","Dinictis felina","Dinictis major","Hoplophoneus primaevus","Nimravus brachyops","Agriotherium africanum","Arctodus simus","Ursus arctos (Pleistocene)")
v_indices_to_remove<-which(sapply(dfcarnivores_postcranial[,2],function(x,v_species_to_remove){is.element(x,v_species_to_remove)},v_species_to_remove=v_species_to_remove,simplify=TRUE))
dfcarnivores_postcranial<-dfcarnivores_postcranial[-v_indices_to_remove,]
m_names_change<-rbind(c("Alopex lagopus","Vulpes lagopus"),c("Lycalopex gymnocerus","Lycalopex gymnocercus"),c("Caracal serval","Leptailurus serval"),c("Felis silvestris libyca","Felis silvestris"),c("Atilax palundinosus","Atilax paludinosus"),c("Gallerella pulverulenta","Galerella pulverulenta"),c("Gallerella sanguinea","Galerella sanguinea"),c("Conepatus mesoleucus","Conepatus leuconotus"),c("Mephitis macroura vittata","Mephitis macroura"),c("Aonyx cinereus","Aonyx cinerea"),c("Paradoxurus hemaphroditus","Paradoxurus hermaphroditus"))
for (i in 1:nrow(m_names_change)){
dfcarnivores_postcranial[which(dfcarnivores_postcranial[,2]==m_names_change[i,1]),2]<-m_names_change[i,2]
}
dfcarnivores_postcranial[,2]<-gsub(" ", "_", dfcarnivores_postcranial[,2])
row.names(dfcarnivores_postcranial)<-dfcarnivores_postcranial[,2]
dat<-dfcarnivores_postcranial[,c("Ecology","HuL","HuPCL","RaL")]
## -----------------------------------------------------------------------------
head(dat)
## -----------------------------------------------------------------------------
b_correct_tree_download<-FALSE
tryCatch({
phyltree_mammals<-ape::read.tree("http://www.biodiversitycenter.org/files/ttol/8.TTOL_mammals_unsmoothed.nwk")
b_correct_tree_download<-TRUE
},error=function(e){cat("Could not download tree file! No analysis can be done! Vignette not built!")},
warning=function(w){cat("Problem with downloading tree file! No analysis can be done! Vignette not built!")}
)
## -----------------------------------------------------------------------------
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Uncia_uncia")]<-"Panthera_uncia"
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Parahyaena_brunnea")]<-"Hyaena_brunnea"
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Bdeogale_crassicauda")]<-"Bdeogale_jacksoni"
tips_todrop<-setdiff(phyltree_mammals$tip.label,rownames(dat))
PrunedTree<-ape::drop.tip(phyltree_mammals,tips_todrop)
Tree<-ape::di2multi(PrunedTree)
dat<-dat[Tree$tip.label,]
## -----------------------------------------------------------------------------
mvSLOUCH::phyltree_paths(Tree)$tree_height
## -----------------------------------------------------------------------------
tree_height<-mvSLOUCH::phyltree_paths(Tree)$tree_height
ScaledTree<-Tree
ScaledTree$edge.length<-ScaledTree$edge.length/tree_height
mvSLOUCH::phyltree_paths(ScaledTree)$tree_height
## -----------------------------------------------------------------------------
isTRUE(all.equal(ScaledTree$tip.label,rownames(dat)))
## -----------------------------------------------------------------------------
row.names(dat)==ScaledTree$tip.label
## -----------------------------------------------------------------------------
regimes<-dat$Ecology[order(match(row.names(dat), ScaledTree$tip.label))]
## -----------------------------------------------------------------------------
regimesFitch<-mvSLOUCH::fitch.mvsl(ScaledTree,regimes)
## -----------------------------------------------------------------------------
regimesFitch<-mvSLOUCH::fitch.mvsl(ScaledTree,regimes,deltran=TRUE)
## -----------------------------------------------------------------------------
reg.col<-regimesFitch$branch_regimes
reg.col[reg.col=="generalist"]<-"purple"
reg.col[reg.col=="arboreal"]<-"green"
reg.col[reg.col=="cursorial"]<-"red"
reg.col[reg.col=="scansorial"]<-"orange"
reg.col[reg.col=="semiaquatic"]<-"blue"
reg.col[reg.col=="semifossorial"]<-"brown"
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# plot(ScaledTree, cex = 1, edge.color = reg.col, edge.width=3.5, type="fan", font=4)
## ----eval=TRUE, echo=FALSE, out.width = "100%", fig.pos="h"-------------------
knitr::include_graphics("./ScaledTree_fan.png", auto_pdf=TRUE)
## -----------------------------------------------------------------------------
mvData<-data.matrix(dat[,c("HuPCL","RaL","HuL")])
mvData<-log(mvData)
## -----------------------------------------------------------------------------
mvStree<-mvSLOUCH::phyltree_paths(ScaledTree)
## -----------------------------------------------------------------------------
BMestim<-mvSLOUCH::BrownianMotionModel(mvStree,mvData)
## -----------------------------------------------------------------------------
BMestim$ParamsInModel$Sxx
BMestim$ParamsInModel$vX0
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1BM<-mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3))
## -----------------------------------------------------------------------------
OU1BM$FinalFound$ParamsInModel$A
## -----------------------------------------------------------------------------
OU1BM$FinalFound$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OU1BM$FinalFound$ParamSummary$optimal.regression
OU1BM$FinalFound$ParamSummary$evolutionary.regression
## -----------------------------------------------------------------------------
OU1BM$FinalFound$ParamsInModel$mPsi
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUBMestim <- mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3), regimes = regimesFitch$branch_regimes, root.regime = "generalist")
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUBMestim$FinalFound
# OUBMestim$MaxLikFound
## -----------------------------------------------------------------------------
OUBMestim$FinalFound$ParamsInModel$mPsi
## -----------------------------------------------------------------------------
OUBMestim$FinalFound$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OUBMestim$FinalFound$ParamSummary$optimal.regression
OUBMestim$FinalFound$ParamSummary$evolutionary.regression
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1OU <- mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3))
## -----------------------------------------------------------------------------
OU1OU$FinalFound$ParamsInModel$A
## -----------------------------------------------------------------------------
OU1OU$FinalFound$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OU1OU$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## -----------------------------------------------------------------------------
OU1BM$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
OUBMestim$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## -----------------------------------------------------------------------------
OUBMestim$FinalFound$ParamsInModel$A
## -----------------------------------------------------------------------------
OUBMestim$FinalFound$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OU1OU$FinalFound$ParamSummary$evolutionary.regression
## -----------------------------------------------------------------------------
OU1OU$FinalFound$ParamsInModel$mPsi
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUOUestim <- mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3))
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamsInModel$mPsi
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamSummary$evolutionary.regression
## -----------------------------------------------------------------------------
cbind(BMestim$ParamSummary$dof,
OU1BM$FinalFound$ParamSummary$dof,
OU1OU$FinalFound$ParamSummary$dof,
OUBMestim$FinalFound$ParamSummary$dof,
OUOUestim$FinalFound$ParamSummary$dof)
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamsInModel$Syy
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUBMestim.mod <- mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3), regimes = regimesFitch$branch_regimes, root.regime = "generalist", Syytype = "LowerTri", Atype = "Diagonal")
## -----------------------------------------------------------------------------
OUBMestim.mod$MaxLikFound
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1 <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)
## -----------------------------------------------------------------------------
capture.output(OU1,file = "OU1.txt")
## -----------------------------------------------------------------------------
OU1$BestModel$model
## -----------------------------------------------------------------------------
OU1$BestModel$BestModel$ParamSummary$phyl.halflife$halflives
## -----------------------------------------------------------------------------
OU1$BestModel$key.properties$R2_non_phylogenetic_conditional_on_predictors
## -----------------------------------------------------------------------------
OU1$BestModel$BestModel$ParamsInModel$A
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1$BestModel
## -----------------------------------------------------------------------------
OU1$BestModel$i
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1$testedModels[[11]]
## -----------------------------------------------------------------------------
OU1$testedModels[[21]]
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OU1$model.setups
## -----------------------------------------------------------------------------
OU1$BestModel$BestModel$ParamSummary$aic.c
## -----------------------------------------------------------------------------
OU1$BestModel$BestModel$ParamSummary$aic
OU1$BestModel$BestModel$ParamSummary$sic
OU1$BestModel$BestModel$ParamSummary$bic
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUf <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)
## -----------------------------------------------------------------------------
capture.output(OUf, file = "OUf.txt")
## -----------------------------------------------------------------------------
OUf$BestModel$model
## -----------------------------------------------------------------------------
OUf$BestModel$BestModel$ParamSummary$phyl.halflife$halflives
OUf$BestModel$key.properties$R2_non_phylogenetic_conditional_on_predictors
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamSummary$aic.c
OUf$BestModel$BestModel$ParamSummary$aic.c
## -----------------------------------------------------------------------------
OUOUestim$FinalFound$ParamSummary$aic.c-OUf$BestModel$BestModel$ParamSummary$aic.c
OU1$BestModel$BestModel$ParamSummary$aic.c-OUf$BestModel$BestModel$ParamSummary$aic.c
## -----------------------------------------------------------------------------
climb.reg <- regimesFitch$branch_regimes
climb.reg[climb.reg=="arboreal"] <- "climber"
climb.reg[climb.reg=="scansorial"] <- "climber"
## -----------------------------------------------------------------------------
climb.col <- climb.reg
climb.col[climb.col=="generalist"] <- "purple"
climb.col[climb.col=="climber"] <- "green"
climb.col[climb.col=="cursorial"] <- "red"
climb.col[climb.col=="semiaquatic"] <- "blue"
climb.col[climb.col=="semifossorial"] <- "brown"
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# plot(ScaledTree, cex = 1, edge.color = climb.col, edge.width=3.5, type="fan", font=4)
## ----eval=TRUE, echo=FALSE, out.width = "100%", fig.pos="h"-------------------
knitr::include_graphics("./ScaledTree2_fan.png", auto_pdf=TRUE)
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUc <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = climb.reg, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)
## -----------------------------------------------------------------------------
OUc$BestModel$model
## -----------------------------------------------------------------------------
OUc$BestModel$BestModel$ParamSummary$aic.c
## -----------------------------------------------------------------------------
strok.reg<-regimesFitch$branch_regimes
strok.reg[strok.reg=="semiaquatic"]<-"stroker"
strok.reg[strok.reg=="semifossorial"]<-"stroker"
## -----------------------------------------------------------------------------
strok.col<-strok.reg
strok.col[strok.col=="generalist"]<-"purple"
strok.col[strok.col=="stroker"]<-"blue"
strok.col[strok.col=="cursorial"]<-"red"
strok.col[strok.col=="arboreal"]<-"green"
strok.col[strok.col=="scansorial"]<-"orange"
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# plot(ScaledTree, cex = 1, edge.color = strok.col, edge.width=3.5, type="fan", font=4)
## ----eval=TRUE, echo=FALSE, out.width = "100%", fig.pos="h"-------------------
knitr::include_graphics("./ScaledTree3_fan.png", auto_pdf=TRUE)
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUs <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)
## -----------------------------------------------------------------------------
OUs$BestModel$model
## -----------------------------------------------------------------------------
OUs$BestModel$BestModel$ParamSummary$aic.c
## -----------------------------------------------------------------------------
OUf$BestModel$BestModel$ParamSummary$dof
OUs$BestModel$BestModel$ParamSummary$dof
## -----------------------------------------------------------------------------
red.reg<-regimesFitch$branch_regimes
red.reg[red.reg=="arboreal"]<-"climber"
red.reg[red.reg=="scansorial"]<-"climber"
red.reg[red.reg=="semiaquatic"]<-"stroker"
red.reg[red.reg=="semifossorial"]<-"stroker"
## -----------------------------------------------------------------------------
red.col<-red.reg
red.col[red.col=="generalist"]<-"purple"
red.col[red.col=="climber"]<-"green"
red.col[red.col=="cursorial"]<-"red"
red.col[red.col=="stroker"]<-"blue"
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# plot(ScaledTree, cex = 1, edge.color = red.col, edge.width=3.5, type="fan", font=4)
## ----eval=TRUE, echo=FALSE, out.width = "100%", fig.pos="h"-------------------
knitr::include_graphics("./ScaledTree4_fan.png", auto_pdf=TRUE)
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUr <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = red.reg, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)
## -----------------------------------------------------------------------------
OUr$BestModel$model
## -----------------------------------------------------------------------------
OUr$BestModel$BestModel$ParamSummary$aic.c
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUOUstart<-mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3), Atype = "Diagonal", diagA = NULL)
## -----------------------------------------------------------------------------
OUOUstart$FinalFound$ParamsInModel$A
OUOUstart$FinalFound$ParamsInModel$Syy
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OptOUs1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OUOUstart$FinalFound$ParamsInModel$A, Syy = OUOUstart$FinalFound$ParamsInModel$Syy))
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OptOUf1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OUOUstart$FinalFound$ParamsInModel$A, Syy = OUOUstart$FinalFound$ParamsInModel$Syy))
## -----------------------------------------------------------------------------
OptOUs1$FinalFound$LogLik
OptOUf1$FinalFound$LogLik
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OptOUs2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OptOUs1$FinalFound$ParamsInModel$A, Syy = OptOUs1$FinalFound$ParamsInModel$Syy))
# OptOUf2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OptOUf1$FinalFound$ParamsInModel$A, Syy = OptOUf1$FinalFound$ParamsInModel$Syy))
## -----------------------------------------------------------------------------
OptOUs2$FinalFound$LogLik
OptOUf2$FinalFound$LogLik
## -----------------------------------------------------------------------------
OptOUs2$FinalFound$LogLik
OptOUf2$FinalFound$LogLik
## -----------------------------------------------------------------------------
OUs$BestModel$BestModel$LogLik
OUf$BestModel$BestModel$LogLik
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# OUOUreStart<-mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3), Atype = "Diagonal", diagA = NULL)
## -----------------------------------------------------------------------------
OUOUreStart$FinalFound$ParamsInModel$A
OUOUreStart$FinalFound$ParamsInModel$Syy
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# FinalOUs1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OUOUreStart$FinalFound$ParamsInModel$A, Syy = OUOUreStart$FinalFound$ParamsInModel$Syy))
# FinalOUf1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OUOUreStart$FinalFound$ParamsInModel$A, Syy = OUOUreStart$FinalFound$ParamsInModel$Syy))
## -----------------------------------------------------------------------------
FinalOUs1$FinalFound$LogLik
FinalOUf1$FinalFound$LogLik
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# FinalOUs2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype =OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = FinalOUs1$FinalFound$ParamsInModel$A, Syy = FinalOUs1$FinalFound$ParamsInModel$Syy))
# FinalOUf2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = FinalOUf1$FinalFound$ParamsInModel$A, Syy = FinalOUf1$FinalFound$ParamsInModel$Syy))
## -----------------------------------------------------------------------------
FinalOUs2$FinalFound$LogLik
FinalOUf2$FinalFound$LogLik
## -----------------------------------------------------------------------------
FinalOUf2$FinalFound$ParamSummary$aic.c - FinalOUs2$FinalFound$ParamSummary$aic.c
## -----------------------------------------------------------------------------
FinalOUs2$FinalFound$ParamSummary$RSS$R2
## -----------------------------------------------------------------------------
FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval
## -----------------------------------------------------------------------------
DFs<-data.frame(
Ecology=factor(colnames(FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Estimated.Point)),
RaL=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Estimated.Point["RaL",],
upper=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Upper.end["RaL",],
lower=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Lower.end["RaL",]
)
## -----------------------------------------------------------------------------
ggplot2::ggplot(DFs, ggplot2::aes(Ecology, RaL))+
ggplot2::geom_point(size=4, colour=c("green","red","purple","orange","blue")) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=lower,ymax=upper),width=0.1,lwd=1.5, colour=c("green","red","purple","orange","blue"))+
ggplot2::xlab("Locomotor habits")+
ggplot2::ylab("RaL(log)") + ggplot2::coord_flip()
## -----------------------------------------------------------------------------
FinalOUs2$FinalFound$ParamSummary$evolutionary.regression
FinalOUs2$FinalFound$ParamSummary$corr.matrix
## -----------------------------------------------------------------------------
FinalOUs2$FinalFound$ParamSummary$trait.regression
FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix
## ----eval=FALSE, echo=TRUE----------------------------------------------------
# BT<-mvSLOUCH::parametric.bootstrap(estimated.model = FinalOUs2, phyltree = mvStree,
# values.to.bootstrap = c("evolutionary.regression", "corr.matrix", "trait.regression", "conditional.corr.matrix"),
# regimes = strok.reg, root.regime = "generalist", predictors = c(3), numboot = 1000,
# Atype = OUs$BestModel$model$Atype,
# Syytype = OUs$BestModel$model$Syytype,
# diagA = OUs$BestModel$model$diagA)
## -----------------------------------------------------------------------------
BTU.EvoReg <- FinalOUs2$FinalFound$ParamSummary$evolutionary.regression
BTU.EvoReg[] <- 0L
BTL.EvoReg <- BTU.EvoReg
for(i in 1:nrow(FinalOUs2$FinalFound$ParamSummary$evolutionary.regression)){
BT.EvoReg<-quantile(sapply(BT$bootstrapped.parameters$evolutionary.regression,function(x)
x[i]), c(0.025, 0.975))
BTL.EvoReg[i,] <- BT.EvoReg[1]
BTU.EvoReg[i,] <- BT.EvoReg[2]
}
## -----------------------------------------------------------------------------
BTL.EvoReg
BTU.EvoReg
## -----------------------------------------------------------------------------
BTU.CorrMat <- rep(NA,length(as.vector(FinalOUs2$FinalFound$ParamSummary$corr.matrix)))
BTL.CorrMat<-BTU.CorrMat
for(i in 1:length(as.vector(FinalOUs2$FinalFound$ParamSummary$corr.matrix))){
BT.CorrMat<-quantile(sapply(BT$bootstrapped.parameters$corr.matrix,function(x) x[i]),c(0.025,0.975))
BTL.CorrMat[i] <- BT.CorrMat[1]
BTU.CorrMat[i] <- BT.CorrMat[2]
}
BTL.CorrMat <- matrix(BTL.CorrMat, nrow =
nrow(FinalOUs2$FinalFound$ParamSummary$corr.matrix))
BTU.CorrMat <- matrix(BTU.CorrMat, nrow =
nrow(FinalOUs2$FinalFound$ParamSummary$corr.matrix))
dimnames(BTL.CorrMat) <- dimnames(BTU.CorrMat)<-
list(row.names(FinalOUs2$FinalFound$ParamSummary$corr.matrix),
colnames(FinalOUs2$FinalFound$ParamSummary$corr.matrix))
## -----------------------------------------------------------------------------
BTL.CorrMat
BTU.CorrMat
## -----------------------------------------------------------------------------
NA.TrtReg<-lapply(1:length(BT$bootstrapped.parameters$trait.regression), function(x)
rep(NA,length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression))))
BTU.TrtReg <- rep(NA, length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression)))
BTL.TrtReg <- BTU.TrtReg
for(i in 1:length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression))){
BT.TrtReg<-quantile(sapply(relist(unlist(BT$bootstrapped.parameters$trait.regression),
NA.TrtReg), function(x) x[i]), c(0.025, 0.975))
BTL.TrtReg[i] <- BT.TrtReg[1]
BTU.TrtReg[i] <- BT.TrtReg[2]
}
BTL.TrtReg <- relist(BTL.TrtReg, FinalOUs2$FinalFound$ParamSummary$trait.regression)
BTU.TrtReg <- relist(BTU.TrtReg, FinalOUs2$FinalFound$ParamSummary$trait.regression)
## -----------------------------------------------------------------------------
BTL.TrtReg
BTU.TrtReg
## -----------------------------------------------------------------------------
BTU.CondCorr <-
rep(NA,length(as.vector(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix)))
BTL.CondCorr<-BTU.CondCorr
for(i in 1:length(as.vector(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))){
BT.CondCorr<-quantile(sapply(BT$bootstrapped.parameters$conditional.corr.matrix,function(x) x[i]),c(0.025,0.975))
BTL.CondCorr[i] <- BT.CondCorr[1]
BTU.CondCorr[i] <- BT.CondCorr[2]
}
BTL.CondCorr <- matrix(BTL.CondCorr, nrow =
nrow(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))
BTU.CondCorr<-matrix(BTU.CondCorr,nrow =
nrow(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))
dimnames(BTL.CondCorr)<-dimnames(BTU.CondCorr)<-list(row.names(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix),
colnames(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))
## -----------------------------------------------------------------------------
BTL.CondCorr
BTU.CondCorr
## ----conditional_print_treeerror, eval=!b_correct_tree_download, echo=FALSE----
# cat("Error: Could not download tree file! No analysis can be done! Vignette not built!")
## ----conditional_print_Dryaderror, eval=!b_correct_dryad_download, echo=FALSE----
# cat("Error: Could not download data file from Dryad! No analysis can be done! Vignette not built!")
# unlink(temp)
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.