TapeS (version 0.13.2)

NSURvar: estimate variance components for component biomass functions

Description

estimate variance components for component biomass functions

Usage

NSURvar(
  data,
  estBM = NULL,
  comp = NULL,
  interval = "confidence",
  level = 0.95,
  adjVarPar = TRUE,
  as.list = TRUE
)

Value

a data.frame with information on lower and upper bound of required interval as well as the (given) estimate and the respective mean squared error

Arguments

data

data / predictors given for prediction by nsur incl. species code for component biomass function, see BaMap.

estBM

estimated biomass components for which variance information is required, given as data.frame, possibly use df[,,drop=FALSE]

comp

which components are required, see tprBiomass

interval

either none, confidence or prediction

level

Tolerance / confidence level, defaults 0.95

adjVarPar

should the variance information be taken from stable models? defaults to TRUE

as.list

Should the return value be a list or rbind to a data.frame? Defaults to TRUE.

Details

Estimates confidence and prediction intervals according to the methods presented in Parresol (2001).

In case, adjVarPar = TRUE, the models with instable variance estimates like Silver fir, Scots pine, Maple and Ash are, firstly, fitted by Norway spruce and European beech, respectively, and, secondly, adjusted to the expected value of the species specific model by substracting the difference to the first model. With that, more stable and imho more realistic confidence and prediction intervals are given. True, this assumes comparability of the variances between species.

Examples

d1 <- seq(42, 56, 2)
h <- estHeight(d1, 1)
data <- data.frame(spp = 1:8, # from BaMap(1, 7)
                   dbh = d1,
                   ht = h,
                   sth = 0.01*h,
                   D03 = 0.8 * d1,
                   kl = 0.7 * h)
estBM <- as.data.frame( nsur(spp = data$spp,
                             dbh = data$dbh,
                             ht = data$ht,
                             sth = data$sth,
                             d03 = data$D03,
                             kl = data$kl) )
estBM$agb <- rowSums(estBM[, -which(colnames(estBM)=="id")])
comp = c("sw", "agb")
interval = "confidence"
level = 0.95
adjVarPar = TRUE
e1 <- TapeS:::NSURvar(data, estBM, comp, interval="confidence", level=0.95, adjVarPar = TRUE)
e2 <- TapeS:::NSURvar(data, estBM, comp, interval="confidence", level=0.95, adjVarPar = FALSE)

if (FALSE) {
par(mfrow=c(1, 2))
plot(x = data$dbh, y = e1$agb_ECBM, main="adjusted Var-Parameter", pch=data$spp,
     ylim=c(0.5*min(e1$agb_ECBM), 1.2*max(e1$agb_ECBM)), las=1,
     ylab="estimated AGB", xlab = "DBH [cm]")
invisible(sapply(1:nrow(e1), function(a){
  # a <- 1
  # lines(x = rep(data$dbh[a], 2), y = c(e2$agb_lwr[a], e2$agb_upr[a]),
  #       col="blue", lwd=2)
  rect(xleft = data$dbh[a] - 0.1, xright = data$dbh[a] + 0.1,
       ybottom = e2$agb_lwr[a], ytop = e2$agb_upr[a], border = "blue")
  lines(x = rep(data$dbh[a], 2), y = c(e1$agb_lwr[a], e1$agb_upr[a]),
        col="red", lwd=2)
}))
legend("bottomright", legend=c("Fi", "Ta", "Kie", "Dgl", "Bu", "Ei", "BAh", "Es"), pch=1:8)

## prediction intervals
e1 <- TapeS:::NSURvar(data, estBM, comp, interval="prediction", level=0.95, adjVarPar = TRUE)
e2 <- TapeS:::NSURvar(data, estBM, comp, interval="prediction", level=0.95, adjVarPar = FALSE)

plot(x = data$dbh, y = e1$agb_ECBM, main="adjusted Var-Parameter", pch=data$spp,
     ylim=c(0, 2*max(e1$agb_ECBM)), las=1,
     ylab="estimated AGB", xlab = "DBH [cm]")
invisible(sapply(1:nrow(e1), function(a){
  # a <- 1
  # lines(x = rep(data$dbh[a], 2), y = c(e2$agb_lwr[a], e2$agb_upr[a]),
  #       col="blue", lwd=2)
  rect(xleft = data$dbh[a] - 0.1, xright = data$dbh[a] + 0.1,
       ybottom = e2$agb_lwr[a], ytop = e2$agb_upr[a], border = "blue")
  lines(x = rep(data$dbh[a], 2), y = c(e1$agb_lwr[a], e1$agb_upr[a]),
        col="red", lwd=2)
}))
legend("topleft", legend=c("Fi", "Ta", "Kie", "Dgl", "Bu", "Ei", "BAh", "Es"), pch=1:8)

## one species, large diameter range
spp <- 1 # spruce
spp <- 5 # beech
spp <- 2 # silver fir
spp <- 8 # ash
d1 <- seq(7, 80, 2)
h <- estHeight(d1, spp)
data <- data.frame(spp = spp,
                   dbh = d1,
                   ht = h,
                   sth = 0.01*h,
                   D03 = 0.8 * d1,
                   kl = 0.7 * h)
estBM <- as.data.frame( nsur(spp = data$spp,
                             dbh = data$dbh,
                             ht = data$ht,
                             sth = data$sth,
                             d03 = data$D03,
                             kl = data$kl) )
estBM$agb <- rowSums(estBM[, -which(colnames(estBM)=="id")])
comp = c("sw", "agb")
interval = "confidence"
level = 0.95
adjVarPar = TRUE
e1 <- TapeS:::NSURvar(data, estBM, comp, interval="confidence", level=0.95, adjVarPar = TRUE)
e2 <- TapeS:::NSURvar(data, estBM, comp, interval="confidence", level=0.95, adjVarPar = FALSE)

par(mfrow=c(1, 2))
plot(x = data$dbh, y = e1$agb_ECBM, main="adjusted Var-Parameter", pch=data$spp,
     ylim=c(0.5*min(e1$agb_ECBM), 1.2*max(e1$agb_ECBM)), las=1,
     ylab="estimated AGB", xlab = "DBH [cm]")
invisible(sapply(1:nrow(e1), function(a){
  # a <- 1
  # lines(x = rep(data$dbh[a], 2), y = c(e2$agb_lwr[a], e2$agb_upr[a]),
  #       col="blue", lwd=2)
  rect(xleft = data$dbh[a] - 0.1, xright = data$dbh[a] + 0.1,
       ybottom = e2$agb_lwr[a], ytop = e2$agb_upr[a], border = "blue")
  lines(x = rep(data$dbh[a], 2), y = c(e1$agb_lwr[a], e1$agb_upr[a]),
        col="red", lwd=2)
}))


## prediction intervals
e1 <- TapeS:::NSURvar(data, estBM, comp, interval="prediction", level=0.95, adjVarPar = TRUE)
e2 <- TapeS:::NSURvar(data, estBM, comp, interval="prediction", level=0.95, adjVarPar = FALSE)

plot(x = data$dbh, y = e1$agb_ECBM, main="adjusted Var-Parameter", pch=data$spp,
     ylim=c(0, 2*max(e1$agb_ECBM)), las=1,
     ylab="estimated biomass", xlab = "DBH [cm]")
invisible(sapply(1:nrow(e1), function(a){
  # a <- 1
  # lines(x = rep(data$dbh[a], 2), y = c(e2$agb_lwr[a], e2$agb_upr[a]),
  #       col="blue", lwd=2)
  rect(xleft = data$dbh[a] - 0.1, xright = data$dbh[a] + 0.1,
       ybottom = e2$agb_lwr[a], ytop = e2$agb_upr[a], border = "blue")
  lines(x = rep(data$dbh[a], 2), y = c(e1$agb_lwr[a], e1$agb_upr[a]),
        col="red", lwd=2)
}))
}