Stomach Data plots

[lentinj]

Add plots for __cons stomach data. @vbartolino already has some examples for doing so.

[vbartolino]

These are simple plotting functions that I made time ago to visualise the fitting of stomach data in gadget model.
To keep in mind that they are based on:

• g2 models and output from Rgadget::gadget.fit (although that's pretty similar structure to gadgetutils::g3_fit)
• stomach likelihood data are usually broken into two components, (LH1) prey-species proportion in the stomachs and (LH2) prey length proportions in the stomachs

ggPropPreyStomach <- function(Gfit, compName, yrs=NULL, q=NULL, area=NULL, pred=NULL, preyList=NULL, ...){

    Gfit$stomachcontent$species <- substring(Gfit$stomachcontent$prey,1,3)
    tmp <- filter(Gfit$stomachcontent, component==compName)
    if(!is.null(yrs)==TRUE){
        tmp <- filter(tmp, year %in% yrs)
    } else { NULL }
    if(!is.null(q)==TRUE){
        tmp <- filter(tmp, step %in% q)
    } else { NULL }
    if(!is.null(area)==TRUE){
        tmp <- filter(tmp, area %in% area)
    } else { NULL }
    if(!is.null(pred)==TRUE){
        tmp <- filter(tmp, predator %in% pred)
    } else { NULL }
    if(!is.null(preyList)==TRUE){
        tmp <- filter(tmp, species %in% preyList)
    } else { NULL }

    if(is.null(tmp$model)){
    ggplot(tmp, aes(year,predicted)) +
        geom_line(aes(color=species)) +
        geom_point(data=tmp, aes(year,observed,color=species)) +
        facet_wrap(~predator+step) # at the moment only single area
    } else {
     ggplot(tmp, aes(year,predicted)) +
        geom_line(aes(color=species,lty=model)) +
        geom_point(data=tmp, aes(year,observed,color=species)) +
        facet_wrap(~predator+step) # at the moment only single area
    }
}

Example

ggPropPreyStomach(fit, compName="spp.cod.stom") + facet_grid(predator~step)

---

ggLenStomach <- function(Gfit, compName, yrs=NULL, q=NULL, area=NULL, predList, preyList=NULL, ...){

    ## Gfit:     gadget.fit obj
    ## compName: name of the likelihood component with data on prey size
    ## yrs:      vector of years (default all)
    ## q:        vector of timesteps (default all)
    ## area:     area (default all)
    ## predList: vector of predators (depending on the name of predator length aggregations, ie "len20", "len40")
    ## predList: vector of preys (3 first letters of the name of the prey length aggregation, ie "spr")

    Gfit$stomachcontent$species <- substring(Gfit$stomachcontent$prey,1,3)
    tmp <- filter(Gfit$stomachcontent, component==compName &
                                       predator %in% predList &
                                       prey.length!=0)
    if(!is.null(yrs)==TRUE){
        tmp <- filter(tmp, year %in% yrs)
    } else { NULL }
    if(!is.null(q)==TRUE){
        tmp <- filter(tmp, step %in% q)
    } else { NULL }
    if(!is.null(area)==TRUE){
        tmp <- filter(tmp, area %in% area)
    } else { NULL }
    if(!is.null(preyList)==TRUE){
        tmp <- filter(tmp, species %in% preyList)
    } else { NULL }

    if(is.null(tmp$model)){
     ggplot(tmp, aes(prey.length,predicted)) +
        geom_line(aes(color=species)) +
        geom_point(data=tmp, aes(prey.length,observed,color=species)) +
        facet_wrap(~year+step+area+predator) # at the moment only single area
    } else {
     ggplot(tmp, aes(prey.length,predicted)) +
        geom_line(aes(color=species,lty=model)) +
        geom_point(data=tmp, aes(prey.length,observed,color=species)) +
        facet_wrap(~year+step+area+predator) # at the moment only single area
     }
}

Example

ggLenStomach(fit, compName="len.cod.stom",
               yrs=2010:2014, q=c(1,4),
               predList=c("cod36"), preyList=c("spr","her")) +
      facet_grid(step~year)

In addition, I've also been using simple scatterplots of deviations and obs vs pred like these below

ggplot(fit$stomachcontent %>% filter(component=="len.cod.stom")) +
    geom_point(aes(prey.length, observed-predicted, group=pred.length,col=as.factor(pred.length))) +
    facet_grid(step~.)

ggplot(fit$stomachcontent %>% filter(component=="len.cod.stom" & year %in% 2010:2014)) +
    geom_point(aes(observed,predicted,pch=substring(prey,1,3))) +
    geom_abline(slope=1, intercept=0, col="red", lwd=0.5) +
    facet_wrap(~year) + xlim(0,1) + ylim(0,1)

[vbartolino]

in my view, there are few additional plots of general interest from the model predictions that I've been using when predation is switched on:

• a time series plot of predation mortality (see below ggPreyCompDiet). It's not uncommon to have it with all the other sources of mortality such as fishing mortality and background M
• a time series plot of prey biomass by prey species (see below ggPredMortality)
• a plot with the functional form of the predator-prey suitability function (usually hand made from the parameters and suitability function)

ggPreyCompDiet <- function(Gfit, predList, yrs=NULL, q=NULL, area=NULL, yearly=TRUE, biomass=TRUE, ...){

    tmp <- Gfit$predator.prey %>%
        filter(predator %in% predList)

# time-area selection    
    if(!is.null(yrs)==TRUE){
        tmp <- filter(tmp, year %in% yrs)
    } else { NULL }
    if(!is.null(q)==TRUE){
        tmp <- filter(tmp, step %in% q)
    } else { NULL }
    if(!is.null(area)==TRUE){
        tmp <- filter(tmp, area %in% area)
    } else { NULL }

# biomass or number    
    if(biomass==TRUE){
        tmp <- tmp %>% mutate(y = biomass_consumed)
        ylabel <-  "Biomass"
    }
    if(biomass==FALSE){
        tmp <- tmp %>% mutate(y = number_consumed)
        ylabel <-  "Numbers"
    } else {NULL}


# timestep or annual aggregation and plotting
    if(yearly==TRUE){
        tmp <- tmp %>%
            group_by(year,area,prey) %>%
            summarise(y=sum(y)) %>%
            filter(y>0)
        p1 <- ggplot(tmp) +
            ylab(ylabel) +
            geom_bar(aes(year,y,fill=prey), stat="identity", position="stack")
    }
    if(yearly==FALSE){
        tmp <- tmp %>%
            group_by(year,step,area,prey) %>%
            summarise(y=sum(y)) %>%
            filter(y>0)
        p1 <- ggplot(tmp) +
            geom_bar(aes(year,y,fill=prey), stat="identity", position="stack") +
            ylab(ylabel) +
            facet_wrap(~step)
    } else {NULL}
    p1
}

ggPredMortality <- function(Gfit, predList, preyList=NULL, preyLen=NULL, yrs=NULL, q=NULL, area=NULL, yearly=TRUE, ...){

    tmp <- Gfit$predator.prey %>%
        filter(predator %in% predList)
    
# time-area selection    
    if(!is.null(yrs)==TRUE){
        tmp <- filter(tmp, year %in% yrs)
    } else { NULL }
    if(!is.null(q)==TRUE){
        tmp <- filter(tmp, step %in% q)
    } else { NULL }
    if(!is.null(area)==TRUE){
        tmp <- filter(tmp, area %in% area)
    } else { NULL }

# prey selection    
    if(!is.null(preyList)==TRUE){
        tmp <- filter(tmp, prey %in% preyList)
    } else { NULL }

# prey length    
    if(!is.null(preyLen)==TRUE){
        tmp <- filter(tmp, length >= preyLen[1] & length <= preyLen[2])
    } else { NULL }

# timestep or annual aggregation and plotting
    if(yearly==TRUE){
        tmp <- tmp %>%
            group_by(prey,year,step,length) %>%
            summarise(mortality=sum(mortality)) %>% # sum contribution of different predators to M
            group_by(prey,year,step) %>%
            summarise(mortality=mean(mortality)) %>% # avg contribution over different lengths to M
            group_by(prey,year) %>%
            summarise(mortality=mean(mortality)) # sum contribution of different timesteps to M
        p1 <- ggplot(tmp) +
            geom_line(aes(year,mortality,col=prey))
    }
    if(yearly==FALSE){
        tmp <- tmp %>%
            group_by(prey,year,step,length) %>%
            summarise(mortality=sum(mortality)) %>% # sum contribution of different predators to M
            group_by(prey,year,step) %>%
            summarise(mortality=mean(mortality)) # avg contribution over different lengths to M
        p1 <- ggplot(tmp) +
            geom_line(aes(year,mortality,col=prey)) +
            facet_wrap(~step)
    } else {NULL}
    p1
}