R/rmse_basis.R

In inbo/dhcurve: Automated Modelling of Diameter Height Curves for Trees

#' @title Berekent RMSE van basismodel
#'
#' @description
#' Deze functie berekent de RMSE door cross-validatie op basis van 6 subsets.
#' Deze functie kan ook gebruikt worden voor het lokaal model (ze bepaalt het
#' verschil tussen de datasets op basis van het al dan niet aanwezig zijn van
#' een veld `DOMEIN_ID` in de dataset).  Opgelet!  In tegenstelling tot de
#' meeste functies van dit package werkt deze functie op basis van de
#' meetgegevens van 1 model.
#' Zie voorbeelden voor een methode om deze functie te kunnen
#' toepassen vertrekkend van meetgegevens (bv. `Data.lokaal`) of vertrekkend van
#' een model waar de meetgegevens uit gehaald kunnen worden (bv. `Basismodel`).
#'
#' Deze functie berekent de RMSE op basis van testgroepen en omvat de volgende
#' deelstappen:
#'
#' - metingen opdelen in 6 testgroepen (veld testgroep)
#'
#' - modellen fitten voor testgroepen, waarbij ze de functie `fit.basis()`
#' 6 keer oproept
#'
#' - RMSE berekenen voor domeinmodellen en Vlaams model op basis van
#' de gemeten waarden en schattingen voor de testgroepen
#'
#' @param Data Meetgegevens van één boomsoort-domeincombinatie (dataframe
#' zoals de dataframes die in de list teruggegeven worden door de functie
#' `initiatie()`)
#' @param Typemodel "Basis" of "Lokaal"?
#' @param BMS Boomsoort
#'
#' @return Dataframe met RMSE_domein (`rmseD`), RMSE_Vlaams (`rmseVL`, niet voor
#' lokaal model) en `maxResid`
#'
#' @examples
#' library(dplyr)
#'
#' #Dataset inladen voor het basismodel
#' Data <- testdataset()
#' Datalijst <- initiatie(Data)
#' Data.basis <- Datalijst[["Basis"]]
#'
#' #De RMSE berekenen voor een basismodel op basis van de dataset
#' Data.basis %>%
#'   group_by(
#'     BMS
#'   ) %>%
#'   do(
#'     rmse.basis(., "Basis", .data$BMS)
#'   ) %>%
#'   ungroup()
#'
#' #Dataset inladen voor het lokaal model
#' Data.lokaal <- Data.basis %>%
#'   filter(DOMEIN_ID == "A")
#'
#' #De rmse berekenen voor een lokaal model
#' Data.lokaal %>%
#'   group_by(
#'     BMS,
#'     DOMEIN_ID
#'   ) %>%
#'   do(
#'     rmse.basis(., "Lokaal", .data$BMS)
#'   ) %>%
#'   ungroup()
#'
#' @export
#'
#' @importFrom dplyr %>% group_by ungroup transmute mutate bind_rows
#' summarise arrange row_number
#' @importFrom plyr .
#' @importFrom rlang .data
#' @importFrom nlme fixef
#' @importFrom stats predict
#' @importFrom assertthat assert_that
#'

rmse.basis <- function(Data, Typemodel, BMS) {

  #controle
  assert_that(is.character(Typemodel))
  Typemodel <- tolower(Typemodel)
  assert_that(Typemodel %in% c("basis", "lokaal"))

  invoercontrole(Data, "fit")
  if (Typemodel == "lokaal") {
    assert_that(length(unique(Data$DOMEIN_ID)) == 1,
                msg = "Voor een lokaal model mag de dataset Data maar 1
                domein bevatten")
  }

  #testgroepen aanmaken in dataset
  Soortdata <- Data %>%
    arrange(.data$DOMEIN_ID, .data$Omtrek, .data$HOOGTE) %>%
    mutate(
      Testgroep = (row_number(.data$DOMEIN_ID) - 1) %% 6 + 1
    ) %>%
    mutate(BMS = BMS)


  #model fitten voor de 6 testgroepen
  Data_result <- data.frame(NULL)
  for (i in 1:6) {
    Data_test <- Soortdata[Soortdata$Testgroep == i, ]
    Data_model <- Soortdata[Soortdata$Testgroep != i, ]

    if (grepl(Typemodel, "lokaal")) {
      Model <- fit.lokaal(Data_model)$Model[[1]]  #nolint
    } else {
      Model <- fit.basis(Data_model)$Model[[1]]   #nolint
    }

    Data_Boomsoort <- Data_test %>%
      mutate(
        H_Dmodel = predict(Model, newdata = .),
        ResidD = .data$HOOGTE - .data$H_Dmodel,
        ResidD2 = .data$ResidD ^ 2,
        ResidVL2 = 0
      )

    if (grepl(Typemodel, "basis")) {
      Data_Boomsoort <- Data_Boomsoort %>%
        mutate(
          H_VLmodel = as.numeric(fixef(Model)[1]) +
            as.numeric(fixef(Model)[2]) * .data$logOmtrek +
            as.numeric(fixef(Model)[3]) * .data$logOmtrek2,
          ResidVL = .data$HOOGTE - .data$H_VLmodel,
          ResidVL2 = .data$ResidVL ^ 2
        )
    }

    Data_result <- Data_result %>%
      bind_rows(Data_Boomsoort)
  }


  #rmse berekenen
  Rmse.soort <- Data_result[Data_result$Omtrek > 0.50, ] %>%
    group_by(
      .data$BMS,
      .data$DOMEIN_ID,
      .data$nBomen,
      .data$nBomenOmtrek05,
      .data$nBomenInterval,
      .data$nBomenIntervalOmtrek05,
      .data$Q5k,
      .data$Q95k
    ) %>%
    summarise(
      sseD = sum(c(.data$ResidD2)),
      sseVL = sum(c(.data$ResidVL2)),
      maxResid = max(c(.data$ResidD))
    ) %>%
    ungroup() %>%
    transmute(
      .data$BMS,
      .data$DOMEIN_ID,
      .data$nBomen,
      .data$nBomenOmtrek05,
      .data$nBomenInterval,
      .data$nBomenIntervalOmtrek05,
      .data$Q5k,
      .data$Q95k,
      rmseD = sqrt(.data$sseD / (.data$nBomenIntervalOmtrek05 - 2)),
      rmseVL = sqrt(.data$sseVL / (.data$nBomenIntervalOmtrek05 - 2)),
      .data$maxResid
    )

  #voor lokaal model het Vlaams model verwijderen (is gelijkgesteld aan 0)
  if (grepl(Typemodel, "lokaal")) {
    Rmse.soort$rmseVL <- NULL
  }

  return(Rmse.soort)
}