phylosem (version 1.1.4)

phylosem: Fit phylogenetic structural equation model

Description

Fits a phylogenetic structural equation model

Usage

phylosem(
  sem,
  tree,
  data,
  family = rep("fixed", ncol(data)),
  covs = colnames(data),
  estimate_ou = FALSE,
  estimate_lambda = FALSE,
  estimate_kappa = FALSE,
  data_labels = rownames(data),
  tmb_inputs = NULL,
  control = phylosem_control()
)

Value

An object (list) of class `phylosem`. Elements include:

data

Copy of argument data

SEM_model

SEM model parsed from sem using specifyModel or specifyEquations

obj

TMB object from MakeADFun

tree

Copy of argument tree

tmb_inputs

The list of inputs passed to MakeADFun

opt

The output from nlminb

sdrep

The output from sdreport

report

The output from obj$report()

parhat

The output from obj$env$parList() containing maximum likelihood estimates and empirical Bayes predictions

Arguments

sem

structural equation model structure, passed to either specifyModel or specifyEquations and then parsed to control the set of path coefficients and variance-covariance parameters

tree

phylogenetic structure, using class as.phylo

data

data-frame providing variables being modeled. Missing values are inputted as NA. If an SEM includes a latent variable (i.e., variable with no available measurements) then it still must be inputted as a column of data with entirely NA values.

family

Character-vector listing the distribution used for each column of data, where each element must be fixed, normal, binomial, or poisson. family="fixed" is default behavior and assumes that a given variable is measured exactly. Other options correspond to different specifications of measurement error.

covs

optional: a character vector of one or more elements, with each element giving a string of variable names, separated by commas. Variances and covariances among all variables in each such string are added to the model. For confirmatory factor analysis models specified via cfa, covs defaults to all of the factors in the model, thus specifying all variances and covariances among these factors. Warning: covs="x1, x2" and covs=c("x1", "x2") are not equivalent: covs="x1, x2" specifies the variance of x1, the variance of x2, and their covariance, while covs=c("x1", "x2") specifies the variance of x1 and the variance of x2 but not their covariance.

estimate_ou

Boolean indicating whether to estimate an autoregressive (Ornstein-Uhlenbeck) process using additional parameter lnalpha, corresponding to the model="OUrandomRoot" parameterization from phylolm as listed in tools:::Rd_expr_doi("10.1093/sysbio/syu005")

estimate_lambda

Boolean indicating whether to estimate additional branch lengths for phylogenetic tips (a.k.a. the Pagel-lambda term) using additional parameter logitlambda

estimate_kappa

Boolean indicating whether to estimate a nonlinear scaling of branch lengths (a.k.a. the Pagel-kappa term) using additional parameter lnkappa

data_labels

For each row of data, listing the corresponding name from tree$tip.label. Default pulls data_labels from rownames(data)

tmb_inputs

optional tagged list that overrides the default constructor for TMB inputs (use at your own risk)

control

Output from phylosem_control, used to define user settings, and see documentation for that function for details.

Details

Note that parameters logitlambda, lnkappa, and lnalpha if estimated are each estimated as having a single value that applies to all modeled variables. This differs from default behavior in phylolm, where these parameters only apply to the "response" and not "predictor" variables. This also differs from default behavior in phylopath, where a different value is estimated in each call to phylolm during the d-separation estimate of path coefficients. However, it is consistent with default behavior in Rphylopars, and estimates should be comparable in that case. These additional parameters are estimated with unbounded support, which differs somewhat from default bounded estimates in phylolm, although parameters should match if overriding phylolm defaults to use unbounded support. Finally, phylosem allows these three parameters to be estimated in any combination, which is expanded functionality relative to the single-option functionality in phylolm.

Also note that phylopath by default uses standardized coefficients. To achieve matching parameter estimates between phylosem and phylopath, standardize each variable to have a standard deviation of 1.0 prior to fitting with phylosem.

References

**Introducing the package, its features, and comparison with other software (to cite when using phylosem):**

Thorson, J. T., & van der Bijl, W. (In press). phylosem: A fast and simple R package for phylogenetic inference and trait imputation using phylogenetic structural equation models. Journal of Evolutionary Biology. tools:::Rd_expr_doi("10.1111/jeb.14234")

*Statistical methods for phylogenetic structural equation models*

Thorson, J. T., Maureaud, A. A., Frelat, R., Merigot, B., Bigman, J. S., Friedman, S. T., Palomares, M. L. D., Pinsky, M. L., Price, S. A., & Wainwright, P. (2023). Identifying direct and indirect associations among traits by merging phylogenetic comparative methods and structural equation models. Methods in Ecology and Evolution, 14(5), 1259-1275. tools:::Rd_expr_doi("10.1111/2041-210X.14076")

*Earlier development of computational methods, originally used for phlogenetic factor analysis:*

Thorson, J. T. (2020). Predicting recruitment density dependence and intrinsic growth rate for all fishes worldwide using a data-integrated life-history model. Fish and Fisheries, 21(2), 237-251. tools:::Rd_expr_doi("10.1111/faf.12427")

Thorson, J. T., Munch, S. B., Cope, J. M., & Gao, J. (2017). Predicting life history parameters for all fishes worldwide. Ecological Applications, 27(8), 2262-2276. tools:::Rd_expr_doi("10.1002/eap.1606")

*Earlier development of phylogenetic path analysis:*

van der Bijl, W. (2018). phylopath: Easy phylogenetic path analysis in R. PeerJ, 6, e4718. tools:::Rd_expr_doi("10.7717/peerj.4718")

von Hardenberg, A., & Gonzalez-Voyer, A. (2013). Disentangling evolutionary cause-effect relationships with phylogenetic confirmatory path analysis. Evolution; International Journal of Organic Evolution, 67(2), 378-387. tools:::Rd_expr_doi("10.1111/j.1558-5646.2012.01790.x")

*Interface involving SEM `arrow notation` is repurposed from:*

Fox, J., Nie, Z., & Byrnes, J. (2020). Sem: Structural equation models. R package version 3.1-11. https://CRAN.R-project.org/package=sem

*Coercing output to phylo4d depends upon:*

Bolker, B., Butler, M., Cowan, P., de Vienne, D., Eddelbuettel, D., Holder, M., Jombart, T., Kembel, S., Michonneau, F., & Orme, B. (2015). phylobase: Base package for phylogenetic structures and comparative data. R Package Version 0.8.0. https://CRAN.R-project.org/package=phylobase

*Laplace approximation for parameter estimation depends upon:*

Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., & Bell, B. M. (2016). TMB: Automatic differentiation and Laplace approximation. Journal of Statistical Software, 70(5), 1-21. tools:::Rd_expr_doi("10.18637/jss.v070.i05")

Examples

# Load data set
data(rhino, rhino_tree, package="phylopath")

# Run phylosem
model = "
  DD -> RS, p1
  BM -> LS, p2
  BM -> NL, p3
  NL -> DD, p4
"
psem = phylosem( sem = model,
          data = rhino[,c("BM","NL","DD","RS","LS")],
          tree = rhino_tree )

# Convert and plot using phylopath
library(phylopath)
my_fitted_DAG = as_fitted_DAG(psem)
coef_plot( my_fitted_DAG )
plot( my_fitted_DAG )

# Convert to phylo4d to extract estimated traits and Standard errors
# for all ancestors and tips in the tree.
# In this rhino example, note that species are labeled s1-s100
# and ancestral nodes are not named.
(traits_est = as_phylo4d(psem))
(traits_SE = as_phylo4d(psem, what="Std. Error"))

# Convert to sem and plot
library(sem)
my_sem = as_sem(psem)
pathDiagram( model = my_sem,
                  style = "traditional",
                  edge.labels = "values" )
effects( my_sem )

# Plot using semPlot
if( require(semPlot) ){
  myplot = semPlotModel( my_sem )
  semPaths( my_sem,
                   nodeLabels = myplot@Vars$name )
}