R/sensitivity_functions.R

In roliveros-ramos/indiseas: Analysis of Indiseas' model-based indicators

calculateFacets = function(indicators, Fmult, smooth=TRUE, dF=0.01) {
  UseMethod("calculateFacets")
}

calculateFacets.default = function(indicators, Fmult, smooth=TRUE, dF=0.01) {
  
  input = list(indicators=indicators, Fmult=Fmult)
  
  if(all(is.na(indicators))) {
    return(list(output=rep(NA, 11), input=input))
  }
  
  smoothInd = splinefun(x = Fmult, y=indicators)
  
  if(isTRUE(smooth)) {  
    Fmult = seq(from=min(Fmult, na.rm=TRUE), to=max(Fmult, na.rm=TRUE), by=dF)
    indicators = smoothInd(Fmult)
  }
  
  #   indicators = (indicators - indicators[1])/sd(indicators, na.rm=TRUE)
  
  Dind = diff(indicators)/diff(Fmult)
  
  # Facets
  # 1) degree of response
  #--------------------------
  # mean rate of change
  mrc = mean(Dind, na.rm=TRUE)
  # area under the curve
  Fmids = seq(from=min(Fmult, na.rm=TRUE), to=max(Fmult, na.rm=TRUE), by=0.005)
  yValue = smoothInd(0.5*(head(Fmult, -1) + tail(Fmult, -1)))
  auc   = sum(abs(yValue)*diff(Fmids), na.rm=TRUE)
  # 90% quantile of the derivatives
  q90 = quantile(Dind, prob=0.9, na.rm=TRUE)
  
  # 2) shape of response
  #--------------------------
  
  spearman = cor(Fmult, indicators, method = "spearman", use="complete")
  pearson  = cor(Fmult, indicators, method = "pearson", use="complete")
  kendall  = cor(Fmult, indicators, method = "kendall", use="complete")
  
  # absolute value of pearson correlation
  avp = abs(pearson)
  # absolute value Spearman correlation
  avs = abs(spearman)
  # absolute value of kendall's tau
  avk = abs(kendall)
  # absolute value of the sign of the derivative
  avsd = abs(mean(sign(Dind), na.rm=TRUE))
  # Maximum number of consecutive equal signs of the derivative 
  # over the number of derivative 
  running = rle(sign(Dind))
  mnc = max(running$lengths)/length(Dind)
  
  output = c(mrc=mrc, auc=auc, q90=q90, avp=avp, avs=avs, avk=avk, avsd=avsd, 
             mnc=mnc)
  
  
  return(list(output=output, input=input))
  
}

calculateFacets.old = function(indicators, Fmult, smooth=TRUE, dF=0.01) {
  
  input = list(indicators=indicators, Fmult=Fmult)
  
  if(all(is.na(indicators))) {
    return(list(output=rep(NA, 11), input=input))
  }

  smoothInd = splinefun(x = Fmult, y=indicators)
  
  if(isTRUE(smooth)) {  
    Fmult = seq(from=min(Fmult, na.rm=TRUE), to=max(Fmult, na.rm=TRUE), by=dF)
    indicators = smoothInd(Fmult)
  }
  
#   indicators = (indicators - indicators[1])/sd(indicators, na.rm=TRUE)
  
  Dind = diff(indicators)/diff(Fmult)
  
  # Facets
  # 1) degree of response
  #--------------------------
  # maximum value of derivative
  
  mvd = max(Dind)
  
  # area under the absolute value of the curve
  
  Fmids = seq(from=min(Fmult, na.rm=TRUE), to=max(Fmult, na.rm=TRUE), by=0.005)
  yValue = smoothInd(0.5*(head(Fmult, -1) + tail(Fmult, -1)))
  auc   = sum(abs(yValue)*diff(Fmids), na.rm=TRUE)
#   auc = integrate(smoothInd)
  # mean rate of change
  
  mrc = mean(Dind, na.rm=TRUE)
  
  # 2) shape of response
  #--------------------------
  
  spearman = cor(Fmult, indicators, method = "spearman", use="complete")
  pearson  = cor(Fmult, indicators, method = "pearson", use="complete")
  kendall  = cor(Fmult, indicators, method = "kendall", use="complete")
  
  # absolute value Spearman correlation
  avs = abs(spearman)
  # absolute value of pearson correlation
  avp = abs(pearson)
  # absolute value of kendall's tau
  avk = abs(kendall)
  # sd of derivative
  sdd = sd(Dind, na.rm=TRUE)
  
  # 3) direction of change
  #--------------------------
  # sign of spearman correlation and of pearson r
  spc = sign(pearson)
  ssc = sign(spearman)
  skc = sign(kendall)
  # sum/average of signs of derivatives
  
  mds = mean(sign(Dind), na.rm=TRUE)
  
  output = c(mvd=mvd, auc=auc, mrc=mrc, avs=avs, avp=avp, avk=avk, sdd=sdd, 
             spc=spc, ssc=ssc, skc=skc, mds=mds)
  
  input  = list(indicators=indicators, Fmult=Fmult)
  
  return(list(output=output, input=input))
  
}

calculateFacets.data.frame = function(indicator, Fmult, smooth=TRUE, dF=0.01) {

  ind = c("ecosystem", "strategy", "value", "indicator") %in% names(indicator)
  if(!all(ind)) stop("data frame does not match sensitivity data frame variables.")
  
  indValue = apply(tapply(indicators$value, 
                          INDEX = list(DF$ecosystem, DF$fishingStrategy, DF$indicator), 
                          FUN=identity), 1:3, FUN = unlist)
  
  Fmult    = unique(indicators$Fmult) # to check
  
  facets = apply(indValue, 2:4, .calculateFacets, Fmult=Fmult)
  
  rownames(facets) = c("mvd", "auc", "mrc", "avs", "avp", "avk", "sdd", 
                       "spc", "ssc", "skc", "mds")
  
  return(facets)
  
}

.calculateFacets = function(indicator, Fmult, smooth=TRUE, dF=0.01) {
  
  x = calculateFacets(indicator=indicator, Fmult=Fmult, smooth=smooth, dF=dF)
  return(x$output)
  
}