SMUFASA Workflow Example

Jennifer McNichol

2023-11-15

Load Packages

library(QFASA)
library(dplyr)
library(compositions)

Modeling Inputs

Prior to starting make sure that:

Fatty Acid Set

data(FAset)
fa.set = as.vector(unlist(FAset))

Matrix of Predator FA Signatures

data(predatorFAs)
tombstone.info = predatorFAs[,1:4]
predator.matrix = predatorFAs[,5:(ncol(predatorFAs))]
npredators = nrow(predator.matrix)

Matrix of Prey FA Signatures

data(preyFAs)
prey.matrix = preyFAs[,-c(1,3)]

# Selecting 5 prey species to include
spec.red <-c("capelin", "herring", "mackerel", "pilchard", "sandlance")
spec.red <- sort(spec.red)
prey.red <- prey.matrix %>%
  filter(Species %in% spec.red)

Prey Lipid Content

FC = preyFAs[,c(2,3)] 
FC = FC %>%
  arrange(Species)
FC.vec = tapply(FC$lipid,FC$Species,mean,na.rm=TRUE)
FC.red <- FC.vec[spec.red]

Running SMUFASA

smufasa.est <- p.SMUFASA(predator.matrix, prey.red, FC.red, fa.set)

p.SMUFASA Output

The MUFASA output is a list with 4 components:

Calibration Coefficient Estimates

A vector of length equal to the number of FAs used and whose sum is the total number of FAs used. Thos is a set of calibration coefficients common to all predators used.

CalEst <- smufasa.est$Cal_Estimates

Diet Estimates

The diet estimate vector returned by p.SMUFASA represents an overall common diet for all predators in predator.matrix . Note: If you wish to obtain diet estimates for each individual predator see the steps below.

DietEst <- smufasa.est$Diet_Estimates

nll

This is a vector of the negative log likelihood values at each iteration of the optimizer.

nll <- smufasa.est$nll

Var_Epsilon

This is the optimized diagonal values of the variance-covariance matrix of the errors. See reference in help file for details.

VarEps <- smufasa.est$Var_Epsilon

Obtaining Diet Estimates

Once a vector of calibration coefficients is obtained via p.SMUFASA you can pass this vector to p.QFASA or p.MUFASA to obtain individual diet estimates.

QFASA

Q = p.QFASA(predator.matrix=predator.matrix, prey.matrix=prey.matrix, 
            cal.mat=CalEst, dist.meas=2, gamma=1, FC=FC.red, 
            start.val=rep(1,nrow(prey.red)), ext.fa=fa.set)

QFASA Diet Estimates:

DietEst.Q <- Q$'Diet Estimates'

See p.QFASA documnetation for further details.

MUFASA

M <- p.SMUFASA(predator.matrix, prey.red, cal.est, FC.red, fa.set)
DietEst.M <- M$Diet_Estimates

See p.MUFASA documnetation for further details.