#' 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)")
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.