library(QFASA)
library(dplyr)
library(compositions)
Prior to starting make sure that:
.csv
file that contains a single column with a header row with the names of thee fatty acids listed below (see example file “FAset.csv”)..csv
file with the full set of FAs and then add code to subset the FAs you wish to use from that set -> this alternative is useful if you are planning to test multiple sets.data(FAset)
= as.vector(unlist(FAset)) fa.set
.csv
file should be proportions summing to 1.p.SMUFASA
BUT make sure that the same FAs appear in the predator and prey files..csv
file can contain as much tombstone data columns as you like, you must extract the predator FA signatures as separate input in order to run in p.SMUFASA
. For example: in the code below, the predator .csv
file (“predatorFAs.csv
”) has 4 tombstone columns (SampleCode, AnimalCode, SampleGroup, Biopsy). Prior to running p.SMUFASA
, the tombstone (columns 1-4) and FA data (columns 5 onward) are each extracted from the original data frame. The FA data becomes the predator.matrix
(which is passed to p.SMUFASA
) and the tombstone data is retained so that it can be recombined with the model output later on.data(predatorFAs)
= predatorFAs[,1:4]
tombstone.info = predatorFAs[,5:(ncol(predatorFAs))]
predator.matrix = nrow(predator.matrix) npredators
.csv
file should be proportions summing to 1.data(preyFAs)
= preyFAs[,-c(1,3)]
prey.matrix
# Selecting 5 prey species to include
<-c("capelin", "herring", "mackerel", "pilchard", "sandlance")
spec.red <- sort(spec.red)
spec.red <- prey.matrix %>%
prey.red filter(Species %in% spec.red)
FC - rep(1,nrow(prey.matrix))
.= preyFAs[,c(2,3)]
FC = FC %>%
FC arrange(Species)
= tapply(FC$lipid,FC$Species,mean,na.rm=TRUE)
FC.vec <- FC.vec[spec.red] FC.red
<- p.SMUFASA(predator.matrix, prey.red, FC.red, fa.set) smufasa.est
The MUFASA output is a list with 4 components:
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.
<- smufasa.est$Cal_Estimates CalEst
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.
<- smufasa.est$Diet_Estimates DietEst
This is a vector of the negative log likelihood values at each iteration of the optimizer.
<- smufasa.est$nll nll
This is the optimized diagonal values of the variance-covariance matrix of the errors. See reference in help file for details.
<- smufasa.est$Var_Epsilon VarEps
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.
= p.QFASA(predator.matrix=predator.matrix, prey.matrix=prey.matrix,
Q 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:
<- Q$'Diet Estimates' DietEst.Q
See p.QFASA
documnetation for further details.
<- p.SMUFASA(predator.matrix, prey.red, cal.est, FC.red, fa.set)
M <- M$Diet_Estimates DietEst.M
See p.MUFASA
documnetation for further details.