inst/doc/comparison.R

In phylosem: Phylogenetic Structural Equation Model

## ----include = FALSE----------------------------------------------------------
have_packages = all(sapply( c("ggplot2","phylolm","phylopath","Rphylopars","phyr","TreeTools"), FUN=requireNamespace))
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = have_packages  # https://r-pkgs.org/vignettes.html#sec-vignettes-eval-option
)
# devtools::build_rmd("vignettes/comparison.Rmd") # to test build
# setwd( R'(C:\Users\James.Thorson\Desktop\Git\phylosem\vignettes)' ); rmarkdown::render( "comparison.Rmd", rmarkdown::pdf_document()) # to build PDF, and perhaps tinytex::latexmk() to install dependencies

## ----setup, echo=TRUE, warning=FALSE, message=FALSE---------------------------
library(phylosem)

## ----package_warning, include=!have_packages----------------------------------
message("Must install ggplot2, phylopath, phylolm, Rphylopars, TreeTools")

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6----
# Settings
Ntree = 100
sd_x = 0.3
sd_y = 0.3
b0_x = 1
b0_y = 0
b_xy = 1

# Simulate tree
set.seed(1)
tree = ape::rtree(n=Ntree)

# Simulate data
x = b0_x + sd_x * phylolm::rTrait(n = 1, phy=tree)
ybar = b0_y + b_xy*x
y_normal = ybar + sd_y * phylolm::rTrait(n = 1, phy=tree)

# Construct, re-order, and reduce data
Data = data.frame(x=x,y=y_normal)[]

# Compare using BM model
start_time = Sys.time()
plm_bm = phylolm::phylolm(y ~ 1 + x, data=Data, phy=tree, model="BM" )
Sys.time() - start_time
knitr::kable(summary(plm_bm)$coefficients, digits=3)

start_time = Sys.time()
psem_bm = phylosem( sem = "x -> y, p",
          data = Data,
          tree = tree,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time
knitr::kable(summary(psem_bm)$coefficients, digits=3)

# Compare using OU
start_time = Sys.time()
plm_ou = phylolm::phylolm(y ~ 1 + x, data=Data, phy=tree, model="OUrandomRoot" )
Sys.time() - start_time

start_time = Sys.time()
psem_ou = phylosem( sem = "x -> y, p",
          data = Data,
          tree = tree,
          estimate_ou = TRUE,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time

knitr::kable(summary(psem_ou)$coefficients, digits=3)
knitr::kable(summary(plm_ou)$coefficients, digits=3)
knitr::kable(c( "phylolm_alpha"=plm_ou$optpar, 
                "phylosem_alpha"=exp(psem_ou$parhat$lnalpha) ), digits=3)

# Compare using Pagel's lambda
start_time = Sys.time()
plm_lambda = phylolm::phylolm(y ~ 1 + x, data=Data, phy=tree, model="lambda" )
Sys.time() - start_time

start_time = Sys.time()
psem_lambda = phylosem( sem = "x -> y, p",
          data = Data,
          tree = tree,
          estimate_lambda = TRUE,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time

knitr::kable(summary(psem_lambda)$coefficients, digits=3)
knitr::kable(summary(plm_lambda)$coefficients, digits=3)
knitr::kable(c( "phylolm_lambda"=plm_lambda$optpar, 
                "phylosem_lambda"=plogis(psem_lambda$parhat$logitlambda) ), digits=3)

# Compare using Pagel's kappa
start_time = Sys.time()
plm_kappa = phylolm::phylolm(y ~ 1 + x, data=Data, phy=tree, model="kappa", lower.bound = 0, upper.bound = 3 )
Sys.time() - start_time

start_time = Sys.time()
psem_kappa = phylosem( sem = "x -> y, p",
          data = Data,
          tree = tree,
          estimate_kappa = TRUE,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time

knitr::kable(summary(psem_kappa)$coefficients, digits=3)
knitr::kable(summary(plm_kappa)$coefficients, digits=3)
knitr::kable(c( "phylolm_kappa"=plm_kappa$optpar, 
                "phylosem_kappa"=exp(psem_kappa$parhat$lnkappa) ), digits=3)

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6, out.width = "75%"----
# Settings
Ntree = 100
sd_x = 0.3
sd_y = 0.3
b0_x = 1
b0_y = 0
b_xy = 1

# Simulate tree
set.seed(1)
tree = ape::rtree(n=Ntree)

# Simulate data
x = b0_x + sd_x * phylolm::rTrait(n = 1, phy=tree)
ybar = b0_y + b_xy*x
y_normal = ybar + sd_y * phylolm::rTrait(n = 1, phy=tree)
y_pois = rpois( n=Ntree, lambda=exp(y_normal) )

# Construct, re-order, and reduce data
Data = data.frame(x=x,y=y_pois)

# Compare using phylolm::phyloglm
pglm = phylolm::phyloglm(y ~ 1 + x, data=Data, phy=tree, method="poisson_GEE" )
knitr::kable(summary(pglm)$coefficients, digits=3)

#
pglmm = phyr::pglmm_compare(
  y ~ 1 + x,
  family = "poisson",
  data = Data,
  phy = tree )
knitr::kable(summary(pglmm), digits=3)

#
pgsem = phylosem( sem = "x -> y, p",
          data = Data,
          family = c("fixed","poisson"),
          tree = tree,
          control = phylosem_control(quiet = TRUE) )
knitr::kable(summary(pgsem)$coefficients, digits=3)

## ----echo=TRUE, message=FALSE, fig.width=6, fig.height=6, out.width = "75%", eval=FALSE----
#  # Comare using Bayesian implementation in brms
#  library(brms)
#  Amat <- ape::vcv.phylo(tree)
#  Data$tips <- rownames(Data)
#  mcmc <- brm(
#    y ~ 1 + x + (1 | gr(tips, cov = A)),
#    data = Data, data2 = list(A = Amat),
#    family = 'poisson',
#    cores = 4
#  )
#  knitr::kable(fixef(mcmc), digits = 3)
#  
#  # Plot them together
#  library(ggplot2)
#  pdat <- rbind.data.frame(
#    coef(summary(pglm))[, 1:2],
#    data.frame(Estimate = pglmm$B, StdErr = pglmm$B.se),
#    setNames(as.data.frame(fixef(mcmc))[1:2], c('Estimate', 'StdErr')),
#    setNames(summary(pgsem)$coefficients[2:3, 3:4], c('Estimate', 'StdErr'))
#  )
#  pdat$Param <- rep(c('Intercept', 'Slope'), 4)
#  pdat$Method <- rep( c('phylolm::phyloglm', 'phyr::pglmm_compare',
#                       'brms::brm', 'phylosem::phylosem'), each = 2)
#  figure = ggplot(pdat, aes(
#    x = Estimate, xmin = Estimate - StdErr,
#    xmax = Estimate + StdErr, y = Param, color = Method
#  )) +
#    geom_pointrange(position = position_dodge(width = 0.6)) +
#    theme_classic() +
#    theme(panel.grid.major.x = element_line(), panel.grid.minor.x = element_line())

## ----include = FALSE, eval=FALSE----------------------------------------------
#  saveRDS( figure, file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\phylosem\vignettes)',"brms.RDS") )
#  saveRDS( pdat, file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\phylosem\vignettes)',"pdat.RDS") )

## ----eval=have_packages, message=FALSE, fig.width=6, fig.height=6, out.width = "75%", echo=FALSE----
library(ggplot2)
#figure = readRDS( file.path(system.file("brms",package="phylosem"),"brms.RDS") )
#figure
pdat = readRDS( file.path(system.file("brms",package="phylosem"),"pdat.RDS") )
ggplot(pdat, aes(
  x = Estimate, xmin = Estimate - StdErr,
  xmax = Estimate + StdErr, y = Param, color = Method
)) +
  geom_pointrange(position = position_dodge(width = 0.6)) +
  theme_classic() +
  theme(panel.grid.major.x = element_line(), panel.grid.minor.x = element_line())

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6----
# Settings
Ntree = 100
sd_x = 0.3
sd_y = 0.3
b0_x = 1
b0_y = 0
b_xy = 1

# Simulate tree
set.seed(1)
tree = ape::rtree(n=Ntree)

# Simulate data
x = b0_x + sd_x * phylolm::rTrait(n = 1, phy=tree)
ybar = b0_y + b_xy*x
y_normal = ybar + sd_y * phylolm::rTrait(n = 1, phy=tree)
y_binom = rbinom( n=Ntree, size=1, prob=plogis(y_normal) )

# Construct, re-order, and reduce data
Data = data.frame(x=x,y=y_binom)

#
pglmm = phyr::pglmm_compare(
  y ~ 1 + x,
  family = "binomial",
  data = Data,
  phy = tree )
knitr::kable(summary(pglmm), digits=3)

#
pgsem = phylosem( sem = "x -> y, p",
          data = Data,
          family = c("fixed","binomial"),
          tree = tree,
          control = phylosem_control(quiet = TRUE) )
knitr::kable(summary(pgsem)$coefficients, digits=3)

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6----
library(phylopath)
library(phylosem)

# make copy of data that's rescaled
rhino_scaled = rhino
  rhino_scaled[,c("BM","NL","LS","DD","RS")] = scale(rhino_scaled[,c("BM","NL","LS","DD","RS")])

# Fit and plot using phylopath
dag <- DAG(RS ~ DD, LS ~ BM, NL ~ BM, DD ~ NL)
start_time = Sys.time()
result <- est_DAG( DAG = dag,
                    data = rhino,
                    tree = rhino_tree,
                    model = "BM",
                    measurement_error = FALSE )
Sys.time() - start_time
plot(result)

# Fit and plot using phylosem
model = "
  DD -> RS, p1
  BM -> LS, p2
  BM -> NL, p3
  NL -> DD, p4
"
start_time = Sys.time()
psem = phylosem( sem = model,
          data = rhino_scaled[,c("BM","NL","DD","RS","LS")],
          tree = rhino_tree,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time
plot( as_fitted_DAG(psem) )

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6----
library(sem)
library(TreeTools)

# Simulation parameters
n_obs = 50
# Intercepts
a1 = 1
a2 = 2
a3 = 3
a4 = 4
# Slopes
b12 = 0.3
b23 = 0
b34 = 0.3
# Standard deviations
s1 = 0.1
s2 = 0.2
s3 = 0.3
s4 = 0.4

# Simulate data
E1 = rnorm(n_obs, sd=s1)
E2 = rnorm(n_obs, sd=s2)
E3 = rnorm(n_obs, sd=s3)
E4 = rnorm(n_obs, sd=s4)
Y1 = a1 + E1
Y2 = a2 + b12*Y1 + E2
Y3 = a3 + b23*Y2 + E3
Y4 = a4 + b34*Y3 + E4
Data = data.frame(Y1=Y1, Y2=Y2, Y3=Y3, Y4=Y4)

# Specify path diagram (in this case, using correct structure)
equations = "
  Y2 = b12 * Y1
  Y4 = b34 * Y3
"
model <- specifyEquations(text=equations, exog.variances=TRUE, endog.variances=TRUE)

# Fit using package:sem
start_time = Sys.time()
Sem <- sem(model, data=Data)
Sys.time() - start_time

# Specify star phylogeny
tree_null = TreeTools::StarTree(n_obs)
  tree_null$edge.length = rep(1,nrow(tree_null$edge))
  rownames(Data) = tree_null$tip.label

# Fit using phylosem
start_time = Sys.time()
psem = phylosem( data = Data,
          sem = equations,
          tree = tree_null,
          control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time

## ----eval=have_packages, echo=FALSE, results='asis'---------------------------
knitr::kable(coef(Sem, standardized=TRUE)[1:2], digits=3)
knitr::kable(coef(psem, standardized=TRUE)[1:2,], digits=3)

## ----eval=have_packages, echo=FALSE, results='asis'---------------------------
knitr::kable(coef(Sem, standardized=FALSE), digits=3)
knitr::kable(coef(psem, standardized=FALSE), digits=3)

## ----eval=have_packages, echo=TRUE, message=FALSE, fig.width=6, fig.height=6----
library(Rphylopars)

# Format data, within no values for species t1
Data = rhino[,c("BM","NL","DD","RS","LS")]
  rownames(Data) = tree$tip.label
Data['t1',] = NA

# fit using phylopars
start_time = Sys.time()
pars <- phylopars( trait_data = cbind(species=rownames(Data),Data),
                  tree = tree,
                  pheno_error = FALSE,
                  phylo_correlated = TRUE,
                  pheno_correlated = FALSE)
Sys.time() - start_time

# Display estimates for missing values
knitr::kable(cbind( "Estimate"=pars$anc_recon["t1",], "Var"=pars$anc_var["t1",] ), digits=3)

# fit using phylosem
start_time = Sys.time()
psem = phylosem( data = Data,
                 tree = tree,
                 sem = "",
                 covs = "BM, NL, DD, RS, LS",
                 control = phylosem_control(quiet = TRUE) )
Sys.time() - start_time

# Display estimates for missing values
knitr::kable(cbind(
  "Estimate"=as.list(psem$sdrep,"Estimate")$x_vj[ match("t1",tree$tip.label), ],
  "Var"=as.list(psem$sdrep,"Std. Error")$x_vj[ match("t1",tree$tip.label), ]^2
), digits=3)