MUFASA Workflow Example

Holly Steeves

2021-07-20

Load Package

library(QFASA)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(compositions)
## Welcome to compositions, a package for compositional data analysis.
## Find an intro with "? compositions"
## 
## Attaching package: 'compositions'
## The following objects are masked from 'package:stats':
## 
##     cor, cov, dist, var
## The following objects are masked from 'package:base':
## 
##     %*%, norm, scale, scale.default

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 (Fat 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]

Calibration Coefficients

data(CC)
cal.vec = CC[,2]
cal.m = replicate(npredators, cal.vec)
rownames(cal.m) <- CC$FA

Running MUFASA

M <- p.MUFASA(predator.matrix, prey.red, cal.m, FC.red, fa.set)

p.MUFASA Output

The MUFASA output is a list with 3 components:

Diet Estimates

This is a matrix of the diet estimate for each predator (by rows, in the same order as the input file) by the species groups (by column, in the same order as the prey.red file). The estimates are expressed as a proportion (they will sum to 1).

Diet_Estimates <- M$Diet_Estimates

nll

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

nll <- M$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 <- M$Var_Epsilon