rpipetibblesapplyrstatix

error with rstatix::tukey_hsd loop: no applicable method for 'tbl_vars' applied to an object of class "NULL"


given the dataset df

df<-structure(list(stream = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L), levels = c("blendziava", "smeltaite", "sventoji"), class = "factor"), 
    season = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 4L, 4L, 4L, 4L), levels = c("fall", "winter", "spring", 
    "summer"), class = c("ordered", "factor")), event = structure(c(1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = c("drought", 
    "flood"), class = "factor"), c = c(533, 596, 622, 590, 282, 
    313, 310, 293, 285, 292, 280, 262, 357, 360, 348, 337, 424, 
    464, 441, 427, 506, 703, 606, 579, 506, 605, 587, 526, 307, 
    366, 301, 330, 385, 346, 396, 210, 225, 218, 217, 228, 226, 
    270, 130, 136, 135, 143, 150, 153, 151, 222, 232, 234, 235, 
    238, 238, 240, 299, 296, 296, 297, 308, 305, 293, 360, 384, 
    382, 368, 367, 360, 362, 387, 417, 404, 389, 386, 382, 378, 
    386, 390, 397, 360, 190, 195, 191, 203, 175, 176, 173, 177, 
    235, 226, 220, 217, 316, 303, 291, 292, 366, 355, 343, 350, 
    378, 368, 341, 384), ph = c(8.31, 8.04, 8.28, 8.2, 8.35, 
    8.06, 8.02, 8.19, 8.14, 7.98, 7.91, 8.04, 7.91, 7.71, 7.82, 
    7.98, 8.3, 7.97, 8.17, 8.21, 8.05, 7.5, 7.86, 8.06, 7.86, 
    7.68, 7.88, 7.92, 9.34, 9.18, 9.34, 9.27, 9.16, 9.18, 8.97, 
    8.08, 8.05, 8, 8.03, 7, 8.07, 8.38, 7.97, 7.95, 7.97, 8.03, 
    8.05, 8.23, 8.65, 8.07, 8.17, 8.22, 8.21, 8.23, 8.27, 8.34, 
    8.3, 8.54, 8.53, 8.46, 8.49, 8.51, 8.59, 8.2, 8.18, 8.27, 
    8.3, 8.44, 8.46, 8.47, 8.04, 7.94, 8.15, 8.17, 8.29, 8.38, 
    8.39, 8.34, 8.33, 8.38, 8.26, 8.12, 8.09, 8.07, 8.05, 8.1, 
    8.1, 6.03, 8.01, 7.94, 7.94, 8, 8.16, 8.23, 8.08, 8.17, 8.16, 
    7.93, 7.87, 7.87, 7.7, 7.65, 7.71, 7.72, 7.81), t = c(12.9, 
    12.2, 12.1, 12.1, 8.2, 8.3, 8.2, 7.9, 2.5, 2.3, 2.3, 1.9, 
    4.7, 3.6, 3.3, 4, 11.5, 11.9, 10.5, 10.5, 18, 19, 17.6, 16, 
    17.3, 19.1, 17.9, 17, 10.7, 11.3, 11.3, 10.9, 11.3, 11.3, 
    11.2, 4.3, 5, 5.1, 5.2, 5.4, 5.4, 5.4, 1.6, 1.8, 1.9, 2.1, 
    2.1, 2.2, 2.4, 4.9, 4.7, 4.9, 4.2, 4.2, 4.3, 4.6, 10.5, 11.8, 
    11.5, 10.7, 12, 12.2, 11.3, 15.2, 17.5, 17.3, 16.1, 17.2, 
    16.8, 17.6, 17, 19.5, 19, 17.9, 18.6, 18.4, 19.4, 11.3, 11.1, 
    11.5, 11.5, 5.3, 5.3, 5.4, 5.6, 2.5, 2.3, 2.4, 2.4, 4.2, 
    3.6, 3.3, 3.5, 11, 10.6, 9.7, 9.8, 18, 18, 17, 16.6, 19.9, 
    20, 19, 19.6), do = c(9.2, 5.8, 7.8, 9.2, 7.9, 7.5, 8, 9.1, 
    11.7, 11.3, 11.7, 13.5, 13.3, 11.7, 13.2, 13.7, 14.9, 11.4, 
    13.8, 11.6, 7.9, 5, 5, 7.9, 8.2, 4.3, 7, 8.3, 9.9, 9.8, 10.5, 
    10.1, 10.8, 11, 10.6, 11.2, 11.5, 11.6, 11.4, 11.3, 11.9, 
    12.6, 12.4, 12.9, 12.5, 12.2, 12.3, 12.6, 14.7, 13.7, 14.9, 
    14.3, 14.3, 14.5, 14.5, 14.2, 14.2, 15.7, 15, 13.5, 13.4, 
    13.2, 13.9, 10.5, 8.6, 9.9, 10.6, 11.6, 10.6, 9.7, 10, 7.3, 
    9.7, 9.9, 10.6, 10.9, 10.9, 8.7, 9.8, 10.1, 9.8, 9.4, 9, 
    9.5, 9.8, 13, 12.5, 12.1, 12.7, 13.3, 14.4, 14.8, 15.5, 11.8, 
    12.1, 13, 14.4, 8, 8.8, 8.9, 7.3, 5.4, 8.4, 7.8, 6.1)), row.names = c(NA, 
-105L), class = c("tbl_df", "tbl", "data.frame"))

with 4 numeric environmental variables, c,ph,t,do and three factors stream, season, event, I'm running different loops using this formula:

env_vars <- c("c","ph","t","do")
proc_model<-function(sub_df) {
   sapply(env_vars, function(env) {
      model <-tukey_hsd(reformulate("stream",env), data = sub_df)%>% 
         add_significance()
      list(
         mytukey_result = model
      )
   }, simplify=FALSE)
}

The model iterate the particular test on each numeric column, represented by the env_vars vector. Then, I run the function on my df that I grouped in respect to the other two factors season and event with:

 stat.test<-df%>%
   group_by(season,event) %>% 
    proc_model()
stat.test

grouping everything in the stat_test.

I do this procedure for different pipe friendly tests from rstatix::package like anova_test(), kruskal_test() and their post hoc comparisons tukey_hsd() and dunn_test(), respectively.

The problem comes when I particularly do it for tukey_hsd , giving me the error:

Error in UseMethod("tbl_vars") : 
no applicable method for 'tbl_vars' applied to an object of class "NULL" 

I suspect is because while the other tests are returning just a data.frame , the tukey_hsd returns a tibble data.frame . Is there a convenient method to export this output anyway, bypassing the above mentioned error? Thanks in advance!


Solution

  • rstatix::tukey_hsd is a wrapper function for stats::TukeyHSD which according to the description needs "a fitted model object, usually an aov fit". Therefore, rstatix::tukey_hsd works with lm() and aov().

    If you include lm() into your code, the error disappears and you get an output:

    env_vars <- c("c","ph","t","do")
    proc_model<-function(sub_df) {
     sapply(env_vars, function(env) {
       model <-tukey_hsd(lm(reformulate("stream",env), data = sub_df))%>% 
         add_significance()
       list(
         mytukey_result = model
       )
     }, simplify=FALSE)
    }
    

    However, in my opinion this approach requires a much more complex analysis including appropriate GLM or GLMM fits followed by meaningful comparisons.