# Phylogenetic comparative methods in fisheries science

#### James T. Thorson

Phylogenetic comparative methods (PCM) are rarely in fisheries science, perhaps due to a lack of familiarity with methods and software among fisheries scientists. This underuse is surprising, given that fisheries science and management strongly depends on foundational results from comparative methods regarding fish productivity and life-history parameters, e.g., for defining proxies for biological reference points and informing hard-to-estimate demographic rates such as stock-recruit and natural mortality parameters.

## Preparing the Then et al. mortality database

To demonstrate the potential role for PCM in fisheries science, we re-analyze a foundational dataset compiled by Then et al.. We specifically download the file “Mlifehist_ver1.0.csv” and then include a copy as a data object in package phylosem to simplify the following demonstration. This demonstration then provides basic syntax and output from PCM, showing the relationship between natural mortality rate, longevity, and growth parameters.

To do so, we first load a phylogeny for fishes using package fishtree, which includes several versions of a phylogeny for fishes develoepd by Rabosky et al.. We then associate all trait data with a tip label from that phylogeny, and provide convenient names for the modeled variables.

    # Load packages
    library(phylosem)
    library(fishtree)

    # Download tree
    out = fishtree_complete_phylogeny()
    tree = out[[1]]

    # Load data object
    data( Mlifehist_ver1_0 )
    Data = Mlifehist_ver1_0

    # Reformat to match tree$tip.label
    Data$Genus_species = factor( paste0(Data$Genus, "_", Data$Species) )

    # Drop duplicates ... not dealing with variation among stocks within species
    Data = Data[match(unique(Data$Genus_species),Data$Genus_species), ]

    # log-transform to simplify later syuntax
    Data = cbind( Data, "logM" = log(Data[,'M']),
                        "logK" = log(Data[,'K']),
                        "logtmax" = log(Data[,'tmax']),
                        "logLinf" = log(Data[,'Linf']) )

    # Identify species in both datasets
    species_to_use = intersect( tree$tip.label, Data$Genus_species )
    species_to_drop = setdiff( Data$Genus_species, tree$tip.label )

    # Drop tips not present in trait-data
    # Not strictly necessary, but helpful to simplify later plots
    tree = ape::keep.tip( tree, tip=species_to_use )

    # Drop trait-data not in phylogeny
    # Necessary to define correlation among data
    rows_to_use = which( Data$Genus_species %in% species_to_use )
    Data = Data[rows_to_use,]

    # Only include modeled variables in trait-data passed to phylosem
    rownames(Data) = Data$Genus_species
    Data = Data[,c('logM','logK','logtmax','logLinf')]

## Fitting and selecting among phylogenetic structural equation models

We then define a path diagram specifying a set of linkages among variables. In the following, we use a path diagram that ensures that mortality rate is statistically independent of growth, conditional upon a measurement for longevity. This specification ensures that, if longevity is available, then it is the sole information used to predict mortality rate. However, if longevity is not available, the model reverts to predicting mortality from growth parameters.

We then fit this model using phylogenetic structural equation models. We specifically apply a grid-search across the eight models formed by any combination of modeled transformations of the phylogenetic tree. We then use marginal AIC to select a model, and list estimated path coefficients.

    # Specify SEM structure
    sem_structure = "
      logK -> logtmax, b1
      logLinf -> logtmax, b2
      logtmax -> logM, a
    "

    # Grid-search model selection using AIC for transformations
    Grid = expand.grid( "OU" = c(FALSE,TRUE),
                        "lambda" = c(FALSE,TRUE),
                        "kappa" = c(FALSE,TRUE) )
    psem_grid = NULL
    for( i in 1:nrow(Grid)){
      psem_grid[[i]] = phylosem( data=Data,
                       tree = tree,
                       sem = sem_structure,
                       estimate_ou = Grid[i,'OU'],
                       estimate_lambda = Grid[i,'lambda'],
                       estimate_kappa = Grid[i,'kappa'],
                       control = phylosem_control(quiet = TRUE) )
    }

    # Extract AIC for each model and rank-order by parsimony
    Grid$AIC = sapply( psem_grid, \(m) AIC(m) )
    Grid = Grid[order(Grid$AIC,decreasing=FALSE),]

    # Select model with lowest AIC
    psem_best = psem_grid[[as.numeric(rownames(Grid[1,]))]]

OU | lambda | kappa | AIC
--- | --- | --- | ---
FALSE | TRUE | TRUE | 1552.940
FALSE | TRUE | FALSE | 1553.164
TRUE | TRUE | TRUE | 1559.502
TRUE | FALSE | TRUE | 1561.424
TRUE | TRUE | FALSE | 1563.372
FALSE | FALSE | TRUE | 1643.863
TRUE | FALSE | FALSE | 1717.509
FALSE | FALSE | FALSE | 2358.137

Path | VarName | Estimate | StdErr | t.value | p.value
--- | --- | --- | --- | --- | ---
NA | Intercept_logM | 1.668 | 0.318 | 5.250 | 0.000
NA | Intercept_logK | -1.759 | 0.566 | 3.107 | 0.002
NA | Intercept_logtmax | -0.575 | 0.655 | 0.878 | 0.380
NA | Intercept_logLinf | 6.559 | 0.445 | 14.742 | 0.000
logK -> logtmax | b1 | -0.448 | 0.067 | 6.716 | 0.000
logLinf -> logtmax | b2 | 0.407 | 0.086 | 4.731 | 0.000
logtmax -> logM | a | -1.022 | 0.037 | 27.732 | 0.000
logM <-> logM | V[logM] | 0.056 | 0.014 | 3.966 | 0.000
logK <-> logK | V[logK] | 0.105 | 0.026 | 4.047 | 0.000
logtmax <-> logtmax | V[logtmax] | 0.078 | 0.019 | 4.029 | 0.000
logLinf <-> logLinf | V[logLinf] | 0.082 | 0.020 | 4.122 | 0.000

## Visualizing output

Finally, we can convert output to formats from other packages, and use existing and third-party PCM packages to plot, query, and post-process output.

    # Plot path diagram
    my_fitted_DAG = as_fitted_DAG(psem_best)
    plot(my_fitted_DAG, type="color")

    # Total, direct, and indirect effects
    my_sem = as_sem(psem_best)
    effects(my_sem)
    #>
    #>  Total Effects (column on row)
    #>               logK   logtmax    logLinf
    #> logtmax -0.4480233  0.000000  0.4073488
    #> logM     0.4579994 -1.022267 -0.4164192
    #>
    #>  Direct Effects
    #>               logK   logtmax   logLinf
    #> logtmax -0.4480233  0.000000 0.4073488
    #> logM     0.0000000 -1.022267 0.0000000
    #>
    #>  Indirect Effects
    #>              logK logtmax    logLinf
    #> logtmax 0.0000000       0  0.0000000
    #> logM    0.4579994       0 -0.4164192

We also show how to plot output using phylosignal, although it is currently unavailable on CRAN and therefore commented out:

    # Load for plotting,
    # https://r-pkgs.org/vignettes.html#sec-vignettes-eval-option
    library(phylosignal)

    # Plot using phylobase
    my_phylo4d = as_phylo4d( psem_best )
    barplot(my_phylo4d)

## Sensitivity to using taxonomic tree

After analysis, we can also conduct sensivity analyses. Here, we show how to construct a taxonomic tree and use this in place of phylogenetic information. As before, this requires some code to reformat the data and then the statistical analysis is simple to specify. The estimated path coefficients are very similar to estimates when using phyogeny.

    library(ape)
    Data = Mlifehist_ver1_0

    # Make taxonomic factors
    Data$Genus_species = factor( paste0(Data$Genus, "_", Data$Species) )
    Data$Genus = factor( Data$Genus )
    Data$Family = factor( Data$Family )
    Data$Order = factor( Data$Order )

    # Make taxonomic tree
    tree = ape::as.phylo( ~Order/Family/Genus/Genus_species, data=Data, collapse=FALSE)
    tree$edge.length = rep(1,nrow(tree$edge))
    tree = collapse.singles(tree)
    tmp = root(tree, node=ape::Ntip(tree)+1 )

    # Drop duplicates ... not dealing with variation among stocks within species
    Data = Data[match(unique(Data$Genus_species),Data$Genus_species), ]

    # log-transform to simplify later syuntax
    Data = cbind( Data, "logM" = log(Data[,'M']),
                        "logK" = log(Data[,'K']),
                        "logtmax" = log(Data[,'tmax']),
                        "logLinf" = log(Data[,'Linf']) )

    # Only include modeled variables in trait-data passed to phylosem
    rownames(Data) = Data$Genus_species
    Data = Data[,c('logM','logK','logtmax','logLinf')]

    # Fit model
    psem_taxon = phylosem( data=Data,
                     tree = tree,
                     sem = sem_structure,
                     estimate_ou = TRUE,
                     control = phylosem_control(quiet = TRUE) )

    # Plot path diagram
    my_fitted_DAG = as_fitted_DAG(psem_taxon)
    plot(my_fitted_DAG, type="color")