library(QFASA) library(plyr)
Prior to starting make sure that:
Choose from one of three distance measures:
1) KL (Kullback-Leibler) 2) AIT (Aitchison) 3) CSD (Chi-Squared)
dist.meas=1
data(FAset) fa.set = as.vector(unlist(FAset))
The FA signatures in the originating .csv file should sum to 100 or 1.
Each predator signature is a row with the FAs in columns (see example file "predatorFAs.csv").
the FA signatures are subsetted for the chosen FA set (created above) and renormalized during the modelling so there is no need to subset and/or renormalize prior to loading the .csv file or running p.QFASA BUT make sure that the the same FAs appear in the predator and prey files (if a FA appears in one but not the other the code will give you an error).
Unlike the original QFASApack code the predator FA .csv file can contain as much tombstone data in columns as you wish but the predator FA signatures must be extracted as a separate input in order to run in p.QFASA. For example: in the code below the predator .csv file ("predatorFAs.csv") has 4 tombstone columns (SampleCode, AnimalCode, SampleGroup, Biopsy). Prior to running QFASA the tombstone (columns 1-4) and FA data (columns 5 onward) are each extracted from the original data frame. The FA data become the the predator.matrix (which is passed to p.QFASA) and the tombstone data is retained so that it can be recombined with the model output later on.
data(predatorFAs) tombstone.info = predatorFAs[,1:4] predator.matrix = predatorFAs[,5:(ncol(predatorFAs))] # number of predator FA signatures this is used to create the matrix of CC values (see section 6 below) npredators = nrow(predator.matrix)
The FA signatures in the originating .csv file should sum to 100 or 1.
The prey file should contain all of the individual FA signatures of the prey and their lipid contents (where appropriate) - a matrix of the mean values for the FAs (prey.matrix) by the designated prey modelling group is then calculated using the MEANmeth function loaded above.
Like the predator .csv file you can have as many tombstone data columns as required but there must be at least one column that identifies the modelling group to be used (in the example file used below "preyFAs.csv" it is the "Species" column).
Unlike the predator data, the prey data is not subsetted and renormalized during the modelling so the prey file needs to be subsetted for the desired FA set (created above) and renormalized to sum to 1 prior to calculating the mean values (see code below). Example: in the code below the "preyFAs.csv" file has 3 tombstone columns. The full FA set is extracted from the data frame (columns 4 onward), subsetted for the FA set in use and then renormalized over 1. The modelling group names (the "Species" column in this case) is then added back to the subsetted and renormalized data (as the first column) and the average values calculated using the MEANmeth function. Note that for the MEANmeth function to work the modelling group name must be in the first column.
#full file data(preyFAs) #extract prey FA only from data frame and subset them for the FA set designated above prey.sub=(preyFAs[,4:(ncol(preyFAs))])[fa.set] #renormalize over 1 prey.sub=prey.sub/apply(prey.sub,1,sum) #extract the modelling group names from the full file group=as.vector(preyFAs$Species) #add modelling group names to the subsetted and renormalized FAs prey.matrix=cbind(group,prey.sub) #create an average value for the FA signature for each designated modelling group prey.matrix=MEANmeth(prey.matrix)
#numbers are the column which identifies the modelling group, and the column which contains the lipid contents FC = preyFAs[,c(2,3)] FC = as.vector(tapply(FC$lipid,FC$Species,mean,na.rm=TRUE))
data(CC) cal.vec = CC[,2] cal.mat = replicate(npredators, cal.vec)
Q = p.QFASA(predator.matrix, prey.matrix, cal.mat, dist.meas, gamma=1, FC, start.val=rep(1,nrow(prey.matrix)), fa.set)
The QFASA output is a list with 2 components:
This is a matrix of the diet estimate for each predator (by rows, in the same order as the input file) by the modelling groups (by column, in the same order as the prey.matrix file). The estimates are expressed as a proportion (they will sum to 1). In the code below the Diet Estimate matrix is extracted from the QFASA output and the modelling group identities and predator tombstone data (created above) are added to the matrix:
DietEst = Q$'Diet Estimates' #estimates changed from proportions to percentages DietEst = round(DietEst*100,digits=2) DietEst = cbind(tombstone.info,DietEst)
This is a list of lists where each list (one per predator) is itself a list of four outputs:
ModFAS: the value of the modelled FA. These are expressed as proportions (they will sum to 1).
DistCont: the contribution of each FA to the final minimized distance.
PropDistCont: the contribution of each FA to the final minimized distance as a proportion of the total.
MinDist: the final minimized distance in the code below the 'ldply' function from the plyr package is used to compile the lists within 'Additional Measures' into a data frame with one row per predator (in the same order as the input predator matrix) and the values for each of the 4 lists arranged into columns. The 'ldply' function automatically names the columns of the data frame with a concatenation of the originating list name and the FA name so that the 4 sets of outputs can be easily identified within the data frame.
Add.meas = ldply(Q$'Additional Measures', data.frame)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.