NSURvar: estimate variance components for component biomass functions

In TapeS: Tree Taper Curves and Sorting Based on 'TapeR'

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
)

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.

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

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)

## Not run: 
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)
}))

## End(Not run)

TapeS documentation built on April 3, 2025, 11:09 p.m.