IFNdyn: IFN forest dynamics and scenario analysis

Documented in drawSubmodelEffects

#' Draw submodel effects
#'
#' Draws the effect of individual variables on the predictions of a specific submodel of forest dynamics
#'
#' @param species Character vector of species codes to be studied.
#' @param submodel Submodel to be studied. Either \code{"growth"}, \code{"survival"}, \code{"survivalPG"},
#' \code{"ingrowth"},  \code{"ingrowthdisp"},
#' \code{"ingrowthB"},  \code{"ingrowthBdisp"},
#' \code{"ingrowthN"},  \code{"ingrowthNdisp"},
#' \code{"height"} or \code{"heightgrowth"}.
#' @param predictor String with the predictor whose effect is to be studied:
#' \itemize{
#' \item{\code{DBH}: Tree diameter.}
#' \item{\code{H}: Tree height.}
#' \item{\code{H/D}: Tree height to diameter ratio.}
#' \item{\code{N}: Stand density.}
#' \item{\code{Nsp.N}: Proportion of the target species.}
#' \item{\code{G}: Stand's basal area.}
#' \item{\code{DBHincprev}: DBH increment of the previous 10-yr timestep.}
#' \item{\code{BAL}: Basal area of larger trees.}
#' \item{\code{SWHC}: Soil water holding capacity.}
#' \item{\code{Rad}: Annual radiation (daily mean).}
#' \item{\code{Temp}: Mean annual temperature.}
#' \item{\code{Prec}: Mean annual precipitation (mm).}
#' \item{\code{PET}: Mean annual potential evapotranspiration (mm).}
#' \item{\code{P/PET}: Moisture index (P/PET).}
#' \item{\code{slope}: Slope (degrees).}
#' \item{\code{elevation}: Elevation (m).}
#' }
#' @param rem.flat A flag to indicate that species with no response to the chosen predictor should be excluded from the plot
#'
#' @details The function uses species mean values for the non-target predictors. Ingrowth is modelled in two stages. Hence, the effect of particular predictors
#' can be drawn for: (a) the binomial model to predict presence of ingrowth ('ingrowthB' or 'ingrowthBdisp');
#' (b) the density of ingrowth conditional to presence ('ingrowthN' or 'ingrowthNdisp');
#' (c) the density of ingrowth resulting from multiplying the two previous predictions ('ingrowth' or 'ingrowthdist').
#' In addition, two sets of ingrowth models are provided, either without or including dispersal effects.
#'
#' @examples
#' pines = c("21", "22", "23", "24", "25", "26", "27","28","31","17")
#' oaks = c("45", "46", "44", "41", "43", "71","72")
#' other = c("51", "55","60", "73", "3","37", "38", "36", "83", "80")
#' drawSubmodelEffects(species = pines, predictor = "DBH", submodel = "growth")
#' drawSubmodelEffects(species = oaks, predictor = "DBH", submodel = "growth")
#' drawSubmodelEffects(species = other, predictor = "DBH", submodel = "growth")
#'
drawSubmodelEffects<-function(species, submodel="growth", predictor = "DBH", rem.flat = TRUE) {

  submodel = match.arg(submodel, c("growth", "survival", "survivalPG",
                                   "ingrowth", "ingrowthdisp",
                                   "ingrowthB", "ingrowthBdisp",
                                   "ingrowthN", "ingrowthNdisp",
                                   "height", "heightgrowth"))
  predictor = match.arg(predictor, c("DBH","H", "H/D", "N", "Nsp.N","G", "DBHincprev","DBHincprevrel","BAL", "SWHC", "Rad", "Temp",
                                     "Prec", "PET","P/PET", "slope", "elevation"))
  getParamMeanDF<-function(sp, nrow=50) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(paramMeans[rep(i,nrow),-c(1,2)])
    return(NULL)
  }
  getParamMin<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(paramMins[i,-c(1,2)])
    return(NULL)
  }
  getParamMax<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(paramMaxs[i,-c(1,2)])
    return(NULL)
  }
  getParamList<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params[i,-c(1,2)]))
    return(NULL)
  }
  getParam2List<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.list(params2[i,-c(1,2)]))
    return(NULL)
  }
  getParamSpeciesName<-function(sp) {
    for(i in 1:length(s)) if(as.character(sp) %in% s[[i]]) return(as.character(params[i,2]))
    stop("Species not found")
  }

  if(predictor=="DBH") xlab = "Diameter (cm)"
  else if(predictor=="H") xlab = "Tree height (m)"
  else if(predictor=="H/D") xlab = "Tree height to diameter ratio (m/cm)"
  else if(predictor=="N") xlab = "Tree density (ind/ha)"
  else if(predictor=="Nsp.N") xlab = "Proportion of target species [0-1]"
  else if(predictor=="G") xlab = "Basal area (m2/ha)"
  else if(predictor=="DBHincprev") xlab = "Diameter increment in previous 10-yr (cm)"
  else if(predictor=="DBHincprevrel") xlab = "Relative diameter increment in previous 10-yr [0-1]"
  else if(predictor=="BAL") xlab = "Basal area of larger trees (m2/ha)"
  else if(predictor=="SWHC") xlab = "Soil water holding capacity (mm)"
  else if(predictor=="Rad") xlab = "Mean annual radiation (MJ)"
  else if(predictor=="Temp") xlab = "Mean annual temperature (degrees C)"
  else if(predictor=="Prec") xlab = "Mean annual precipitation (mm)"
  else if(predictor=="PET") xlab = "Mean annual potential evapotranspiration (mm)"
  else if(predictor=="P/PET") xlab = "Moisture index (P/PET)"
  else if(predictor=="slope") xlab = "Slope (degrees)"
  else if(predictor=="elevation") xlab = "Elevation (m)"


  x =numeric(0)
  y = numeric(0)
  sp = character(0)
  if(submodel=="growth") {
    params = get("defaultGrowthParams")
    paramMeans = get("defaultGrowthMeans")
    paramMins = get("defaultGrowthMins")
    paramMaxs = get("defaultGrowthMaxs")
    vars = .getVarsGrowth()
    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p = getParamList(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p)) {
        if(predictor=="DBH") {
          df$d =  seq(minvec[1,"d"], maxvec[1,"d"], length.out = nrow(df))
          df$inv.d = 1/df$d
          df$ln.d = log(df$d+1)
          df$sq.d = df$d^2
          df$sqrt.d = sqrt(df$d)
          df$BAL.ln.d = df$BAL/log(df$d+1)
          x1 = df$d
        }
        else if(predictor=="G") {
          df$G =  seq(minvec[1,"G"], maxvec[1,"G"], length.out = nrow(df))
          x1 = df$G
        }
        else if(predictor=="BAL") {
          df$BAL =  seq(minvec[1,"BAL"], maxvec[1,"BAL"], length.out = nrow(df))
          df$BAL.ln.d = df$BAL/log(df$d+1)
          x1 = df$BAL
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Temp") {
          df$Temp =  seq(minvec[1,"Temp"], maxvec[1,"Temp"], length.out = nrow(df))
          x1 = df$Temp
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="slope") {
          df$slope =  seq(minvec[1,"slope"], maxvec[1,"slope"], length.out = nrow(df))
          x1 = df$slope
        }
        df$ln.G = log(df$G)
        df$Temp15sq = (df$Temp-15)^2

        # print(head(df))
        ss = as.numeric(rowSums(sweep(as.matrix(df[,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
        y1 = exp(p[["Intercept"]]+ss)
        # print(y1)
        include = TRUE
        if(rem.flat) {
          if(sd(y1, na.rm=T)==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    ggplot(res)+
      geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
      xlab(xlab)+
      ylab("10-yr diameter increment (cm)")
  }
  else if(submodel=="survival") {
    params = get("defaultSurvivalParams")
    paramMeans = get("defaultSurvivalMeans")
    paramMins = get("defaultSurvivalMins")
    paramMaxs = get("defaultSurvivalMaxs")
    vars = .getVarsSurvival()
    # print(vars)
    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p = getParamList(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p)) {
        if(predictor=="DBH") {
          df$d =  seq(minvec[1,"d"], maxvec[1,"d"], length.out = nrow(df))
          df$ln.d = log(df$d+1)
          df$inv.d = 1/df$d
          df$sq.d = df$d^2
          df$sqrt.d = sqrt(df$d)
          x1 = df$d
        }
        else if(predictor=="H") {
          df$h =  seq(minvec[1,"h"], maxvec[1,"h"], length.out = nrow(df))
          x1 = df$h
        }
        else if(predictor=="G") {
          df$G =  seq(minvec[1,"G"], maxvec[1,"G"], length.out = nrow(df))
          x1 = df$G
        }
        else if(predictor=="BAL") {
          df$BAL =  seq(minvec[1,"BAL"], maxvec[1,"BAL"], length.out = nrow(df))
          x1 = df$BAL
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Temp") {
          df$Temp =  seq(minvec[1,"Temp"], maxvec[1,"Temp"], length.out = nrow(df))
          x1 = df$Temp
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="slope") {
          df$slope =  seq(minvec[1,"slope"], maxvec[1,"slope"], length.out = nrow(df))
          x1 = df$slope
        }
        df$ln.G = log(df$G)
        df$Temp15sq = (df$Temp-15)^2
        # print(head(df))
        ss = as.numeric(rowSums(sweep(as.matrix(df[,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
        y1 = 1/(1+exp(-1*(p[["Intercept"]]+ss)))
        include = TRUE
        if(rem.flat) {
          if(sd(y1, na.rm=T)==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    ggplot(res)+
      geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
      scale_y_continuous(limits=c(0,1))+
      xlab(xlab)+
      ylab("Probability of survival [0-1]")
  }
  else if(submodel=="survivalPG") {
    params = get("defaultSurvivalPGParams")
    paramMeans = get("defaultSurvivalPGMeans")
    paramMins = get("defaultSurvivalPGMins")
    paramMaxs = get("defaultSurvivalPGMaxs")
    vars = .getVarsSurvivalPG()
    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p = getParamList(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p)) {
        if(predictor=="DBH") {
          df$d =  seq(minvec[1,"d"], maxvec[1,"d"], length.out = nrow(df))
          df$ln.d = log(df$d+1)
          df$inv.d = 1/df$d
          df$sq.d = df$d^2
          df$sqrt.d = sqrt(df$d)
          df$ln.d.inc.prev.ln.d = df$ln.d.inc.prev/log(df$d+1)
          x1 = df$d
        }
        else if(predictor=="H") {
          df$h =  seq(minvec[1,"h"], maxvec[1,"h"], length.out = nrow(df))
          x1 = df$h
        }
        else if(predictor=="G") {
          df$G =  seq(minvec[1,"G"], maxvec[1,"G"], length.out = nrow(df))
          df$ln.G = log(df$G)
          x1 = df$G
        }
        else if(predictor=="DBHincprev") {
          x1 = seq(exp(minvec[1,"ln.d.inc.prev"]), exp(maxvec[1,"ln.d.inc.prev"]), length.out = nrow(df))
          df$ln.d.inc.prev =  log(x1+1)
          df$ln.d.inc.prev.ln.d = df$ln.d.inc.prev/log(df$d+1)
        }
        else if(predictor=="DBHincprevrel") {
          x1 = seq(0, 1, length.out = nrow(df))
          df$ln.d.inc.prev.ln.d = x1
          df$ln.d.inc.prev = x1*log(df$d+1)
        }
        else if(predictor=="BAL") {
          df$BAL =  seq(minvec[1,"BAL"], maxvec[1,"BAL"], length.out = nrow(df))
          x1 = df$BAL
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Temp") {
          df$Temp =  seq(minvec[1,"Temp"], maxvec[1,"Temp"], length.out = nrow(df))
          x1 = df$Temp
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="slope") {
          df$slope =  seq(minvec[1,"slope"], maxvec[1,"slope"], length.out = nrow(df))
          x1 = df$slope
        }
        df$ln.G = log(df$G)
        df$Temp15sq = (df$Temp-15)^2
        ss = as.numeric(rowSums(sweep(as.matrix(df[,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
        y1 = 1/(1+exp(-1*(p[["Intercept"]]+ss)))
        include = TRUE
        if(rem.flat) {
          if(sd(y1, na.rm=T)==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    ggplot(res)+
      geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
      scale_y_continuous(limits=c(0,1))+
      xlab(xlab)+
      ylab("Probability of survival [0-1]")
  }
  else if(submodel %in% c("ingrowth", "ingrowthdisp",
                          "ingrowthB", "ingrowthBdisp",
                          "ingrowthN", "ingrowthNdisp")) {
    if(submodel %in% c("ingrowth", "ingrowthB", "ingrowthN")) {
      params = get("defaultIngrowthBParams")
      params2 = get("defaultIngrowthNParams")
      paramMeans = get("defaultIngrowthNMeans")
      paramMins = get("defaultIngrowthNMins")
      paramMaxs = get("defaultIngrowthNMaxs")
      varsB = .getVarsIngrowthB()
      varsN = .getVarsIngrowthN()
    } else if(submodel %in% c("ingrowthdisp", "ingrowthBdisp", "ingrowthNdisp")) {
      params = get("defaultIngrowthBParams_disp")
      params2 = get("defaultIngrowthNParams_disp")
      paramMeans = get("defaultIngrowthNMeans_disp")
      paramMins = get("defaultIngrowthNMins_disp")
      paramMaxs = get("defaultIngrowthNMaxs_disp")
      varsB = .getVarsIngrowthB_disp()
      varsN = .getVarsIngrowthN_disp()
    }

    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p_B = getParamList(spi)
      p_N = getParam2List(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p_B)) {
        if(predictor=="N") {
          df$ln.N =  log(seq(exp(minvec[1,"ln.N"]), exp(maxvec[1,"ln.N"]), length.out = nrow(df)))
          x1 = df$ln.N
        }
        else if(predictor=="Nsp.N") {
          df$Nsp.N =  seq(0.01, 1, length.out = nrow(df))
          df$Nsp.N_s0.5 =  seq(0.01, 1, length.out = nrow(df))
          df$Nsp.N_s1 =  seq(0.01, 1, length.out = nrow(df))
          df$Nsp.N_s2 =  seq(0.01, 1, length.out = nrow(df))
          df$Nsp.N_s4 =  seq(0.01, 1, length.out = nrow(df))
          x1 = df$Nsp.N
        }
        else if(predictor=="G") {
          df$G =  seq(minvec[1,"G"], maxvec[1,"G"], length.out = nrow(df))
          x1 = df$G
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Temp") {
          df$Temp =  seq(minvec[1,"Temp"], maxvec[1,"Temp"], length.out = nrow(df))
          x1 = df$Temp
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="slope") {
          df$slope =  seq(minvec[1,"slope"], maxvec[1,"slope"], length.out = nrow(df))
          x1 = df$slope
        }
        df$ln.G = log(df$G)
        df$Temp15sq = (df$Temp-15)^2

        ss_B = as.numeric(rowSums(sweep(as.matrix(df[,varsB]),2,as.numeric(p_B[varsB]), FUN="*"), na.rm=T))
        B =1/(1+exp(-1*(as.numeric(p_B["Intercept"])+ss_B)))
        ss_N = as.numeric(rowSums(sweep(as.matrix(df[,varsN]),2,as.numeric(p_N[varsN]), FUN="*"), na.rm=T))
        N = exp(p_N[["Intercept"]]+ss_N)
        if(submodel %in% c("ingrowth", "ingrowthdisp")) {y1 = N*B}
        else if(submodel %in% c("ingrowthN", "ingrowthNdisp")) {y1 = N}
        else {y1 = B}
        include = TRUE
        if(rem.flat) {
          if(sd(y1, na.rm=T)==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    if(submodel %in% c("ingrowth", "ingrowthdisp")) {
      ggplot(res)+
        geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
        xlab(xlab)+
        scale_y_log10()+
        ylab("Ingrowth density (ind/ha)")
    } else if(submodel %in% c("ingrowthN", "ingrowthNdisp")) {
      ggplot(res)+
        geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
        xlab(xlab)+
        scale_y_log10()+
        ylab("Ingrowth conditional density (ind/ha)")
    } else {
      ggplot(res)+
        geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
        xlab(xlab)+
        ylab("Ingrowth probability [0-1]")
    }
  }
  else if(submodel=="height") {
    params = get("defaultHeightParams")
    paramMeans = get("defaultHeightMeans")
    paramMins = get("defaultHeightMins")
    paramMaxs = get("defaultHeightMaxs")
    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p = getParamList(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p)) {
        if(predictor=="DBH") {
          df$d =  seq(minvec[1,"d"], maxvec[1,"d"], length.out = nrow(df))
          x1 = df$d
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="elevation") {
          df$elevation =  seq(minvec[1,"elevation"], maxvec[1,"elevation"], length.out = nrow(df))
          x1 = df$elevation
        }
        y1 = (p[["A0"]] + p[["B1"]]*df$PPET+p[["B2"]]*df$SWHC + p[["B3"]]*df$Rad + p[["B4"]]*df$elevation)/(1+(p[["C1"]]/df$d)+(p[["C2"]]/(df$d^2)))
        y1 = pmax(0,y1)
        include = TRUE
        if(rem.flat) {
          if(sd(y1, na.rm=T)==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    ggplot(res)+
      geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
      xlab(xlab)+
      ylab("Tree height (m)")
  }
  else if(submodel=="heightgrowth") {
    params = get("defaultHeightGrowthParams")
    paramMeans = get("defaultHeightGrowthMeans")
    paramMins = get("defaultHeightGrowthMins")
    paramMaxs = get("defaultHeightGrowthMaxs")
    vars = .getVarsHeightGrowth()
    sp_vec = as.character(params$species)
    for(spi in species) {
      s = strsplit(sp_vec, split=",")
      p = getParamList(spi)
      df = getParamMeanDF(spi)
      minvec = getParamMin(spi)
      maxvec = getParamMax(spi)
      if(!is.null(p)) {
        if(predictor=="H") {
          df$h =  seq(minvec[1,"h"], maxvec[1,"h"], length.out = nrow(df))
          df$inv.h = 1/df$h
          df$ln.h = log(df$h+1)
          x1 = df$h
        }
        else if(predictor=="H/D") {
          df$h.d =  seq(minvec[1,"h.d"], maxvec[1,"h.d"], length.out = nrow(df))
          x1 = df$h.d
        }
        else if(predictor=="G") {
          df$G =  seq(minvec[1,"G"], maxvec[1,"G"], length.out = nrow(df))
          x1 = df$G
        }
        else if(predictor=="SWHC") {
          df$SWHC =  seq(minvec[1,"SWHC"], maxvec[1,"SWHC"], length.out = nrow(df))
          x1 = df$SWHC
        }
        else if(predictor=="Rad") {
          df$Rad =  seq(max(10,minvec[1,"Rad"]), maxvec[1,"Rad"], length.out = nrow(df))
          x1 = df$Rad
        }
        else if(predictor=="Temp") {
          df$Temp =  seq(minvec[1,"Temp"], maxvec[1,"Temp"], length.out = nrow(df))
          x1 = df$Temp
        }
        else if(predictor=="Prec") {
          df$Prec =  seq(minvec[1,"Prec"], maxvec[1,"Prec"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$Prec
        }
        else if(predictor=="PET") {
          df$PET =  seq(minvec[1,"PET"], maxvec[1,"PET"], length.out = nrow(df))
          df$PPET = pmin(df$Prec/df$PET, 1.3)
          x1 = df$PET
        }
        else if(predictor=="P/PET") {
          df$PPET =  seq(minvec[1,"PPET"], maxvec[1,"PPET"], length.out = nrow(df))
          x1 = df$PPET
        }
        else if(predictor=="slope") {
          df$slope =  seq(minvec[1,"slope"], maxvec[1,"slope"], length.out = nrow(df))
          x1 = df$slope
        }
        df$ln.G = log(df$G)
        df$Temp15sq = (df$Temp-15)^2

        # print(head(df))
        # print(p[vars])
        ss = as.numeric(rowSums(sweep(as.matrix(df[,vars]),2,as.numeric(p[vars]), FUN="*"), na.rm=T))
        # print(ss)
        y1 = exp(p[["Intercept"]]+ss)
        # print(cbind(x1, y1))
        include = TRUE
        if(rem.flat) {
          sdval = sd(y1, na.rm=T)
          if(is.na(sdval)) include = FALSE
          else if(sdval==0) include = FALSE
        }
        if(include) {
          x = c(x, x1)
          y = c(y, y1)
          sp = c(sp, rep(getParamSpeciesName(spi), length(x1)))
        }
      }
    }
    res = data.frame(x = x, y = y, Species = sp, stringsAsFactors = F)
    ggplot(res)+
      geom_line(mapping=aes(x, y, colour=Species), size=1.5)+
      xlab(xlab)+
      ylab("10-yr height increment (m)")
  }
}