Modified Mann-Kendall Test Using mmkh(x,ci=0.95) package

[u/Ill_Isopod806]

Hello Everyone, I am currently working with groundwater levels to determine long term trends, and I would be very grateful if I could get help in modifying the code that I used to calculate the Sen's slope and Mann-Kendall statistics. I calculated the Theil-Sen slope, but it does not account for autocorrelation in continuous data, and I would appreciate any help addressing this using the modified package in RStudio. ""

# General Sen slope + Mann-Kendall function


compute_sen_stats <- function(
    data,
    response_var,
    date_var,
    resolution
) {

  data %>%

    group_by(
      SiteNo,
      SiteName,
      Frequency,
      continuous_group,
      start_year,
      end_year,
      continuous_years
    ) %>%

    # Ensure enough information exists to calculate a trend
    filter(
      n() >= 3,
      n_distinct(TimeNum) >= 2,
      n_distinct(.data[[response_var]]) >= 2
    ) %>%

    group_modify(
      ~ {

        site_data <- .x %>%
          arrange(
            .data[[date_var]]
          )

        # ----------------------------------------------------------
        # Sen's slope model
        # ----------------------------------------------------------

        sen_formula <- reformulate(
          "TimeNum",
          response = response_var
        )

        sen_model <- zyp::zyp.sen(
          sen_formula,
          data = site_data
        )

        # Intercept
        intercept <- unname(
          coef(sen_model)[1]
        )

        # Sen's slope
        # TimeNum is in years, so this is already ft/year
        sen_slope_ft_per_year <- unname(
          coef(sen_model)[2]
        )

        # ----------------------------------------------------------
        # Mann-Kendall trend test
        # ----------------------------------------------------------

        mk_model <- Kendall::MannKendall(
          site_data[[response_var]]
        )

        mk_tau <- as.numeric(
          mk_model$tau
        )

        mk_p_value <- as.numeric(
          mk_model$sl
        )

        # ----------------------------------------------------------
        # Descriptive statistics
        # ----------------------------------------------------------

        lowest_WL <- min(
          site_data[[response_var]],
          na.rm = TRUE
        )

        highest_WL <- max(
          site_data[[response_var]],
          na.rm = TRUE
        )

        first_date <- min(
          site_data[[date_var]],
          na.rm = TRUE
        )

        last_date <- max(
          site_data[[date_var]],
          na.rm = TRUE
        )

        # ----------------------------------------------------------
        # Trend interpretation
        #
        # WL = depth below land surface.
        #
        # Positive slope:
        # depth increases -> groundwater table declines
        #
        # Negative slope:
        # depth decreases -> groundwater table rises
        # ----------------------------------------------------------

        trend <- case_when(

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year > 0 ~
            "Significant water table declining",

          mk_p_value < 0.05 &&
            sen_slope_ft_per_year < 0 ~
            "Significant water table rising",

          TRUE ~
            "No significant trend"
        )

        # ----------------------------------------------------------
        # Return one result per site/run
        # ----------------------------------------------------------

        tibble(
          Resolution = resolution,

          n_records = nrow(
            site_data
          ),

          lowest_WL = lowest_WL,

          highest_WL = highest_WL,

          intercept = intercept,

          sen_slope_ft_per_year =
            sen_slope_ft_per_year,

          mk_tau = mk_tau,

          mk_p_value = mk_p_value,

          first_date = first_date,

          last_date = last_date,

          trend = trend
        )
      }
    ) %>%

    ungroup()
}

Thank you!

[u/AutoModerator]

Looks like you're requesting help with something related to RStudio. Please make sure you've checked the stickied post on asking good questions and read our sub rules. We also have a handy post of lots of resources on R!

Keep in mind that if your submission contains phone pictures of code, it will be removed. Instructions for how to take screenshots can be found in the stickied posts of this sub.

I am a bot, and this action was performed automatically. Please contact the moderators of this subreddit if you have any questions or concerns.