【问题标题】:Using purrr and modelr to make predictions from dose-response models (drm)使用 purrr 和 modelr 从剂量反应模型 (drm) 进行预测
【发布时间】:2020-09-24 22:06:14
【问题描述】:

我有一个缩写形式的数据集,如下所示:

library(tidyverse)    
dat_s<-tibble(
      type=c(rep("A", 9), rep("B", 8), rep("C", 10)),
      ref=c("ref3", "ref3", "ref1",  "ref2", "ref2", "ref1", "ref2", "ref2", "ref2", "ref2", 
            "ref1", "ref2", "ref2", "ref3", "ref2", "ref3", "ref1", "ref3", 
            "ref2", "ref3", "ref1", "ref1", "ref3", "ref1", "ref1", "ref2", "ref2"),
      info=as.character(sample(100, 27)),
      liv=c(3.0e-05, 2.9e-07, 2.2e-07, 2.7e-07, 2.6e-06, 4.8e-07, 1.4e-05, 2.6e-06, 7.7e-06, 2.2e-06, 
            1.5e-07, 1.6e-07, 1.8e-06, 6.1e-08, 4.9e-06, 4.9e-06, 1.8e-06, 1.5e-07,
            4.3e-08, 1.8e-06, 1.0e-07, 1.6e-07, 9.7e-07, 1.0e-06, 6.4e-07, 1.2e-07, 5.7e-06),
      prod=c(0.00, 2, 3, 4.80, 2.10, 5.10, 0.00, 0.13, 2.00, 0.13, 0.00, 4.10, 4.60, 2.10, 0.26, 0.00, 
             4.60, 0.00, 4.60, 2.10, 4.80, 0.00, 0.00, 1.80, 3.60, 4.10, 0.00)
    )%>%
  mutate(livp1=liv+1)

我想计算type 和ref 的每种组合的剂量反应关系,进行预测以绘制曲线,并计算残差。 info 列是为了反映我在此数据框中还有其他列需要保留,但在剂量反应分析中并不重要。

我首先使用函数和嵌套数据框创建模型:

dr_s<-function(df){drc::drm(data=df, prod~log(livp1), fct=drc::LL.3())}

dat_mods<-
  dat_s%>%
  group_by(type, ref)%>%
  nest() %>%
  mutate(dr_mod=map(data, dr_s))

这有助于创建模型并将它们放入数据框中。要将add_predictions 与drm 类型的模型一起使用,输入必须是data.frame(而不是tibble)。当我尝试为每个 livp1 变量添加预测时(根据下面的 cmets):

dat_mods%>%

mutate(mod_preds=map2(data, dr_mod, ~add_predictions(data=as.data.frame(.x), model=.y))))

我收到 non-numeric argument to binary operator 错误消息。当info 列不是character 类时,此代码可以正常工作。但是,我需要将这些信息与预测数据一起保留,并希望尽可能避免将其从数据框中提取出来。

感谢任何指导!

【问题讨论】:

    标签: r dplyr purrr drc modelr


    【解决方案1】:

    这个很傻,drm 类模型的predict 方法不适用于tibble 类的对象。所以你必须将.x 转换为data.frame。

    dat_s%>%
      group_by(type, ref)%>%
      nest() %>%
      mutate(dr_mod=map(data, dr_s),
             mod_preds=map2(data, dr_mod, 
                            ~modelr::add_predictions(data = as.data.frame(.x[,"livp1"]),
                                                     model = .y)))
    ## A tibble: 9 x 5
    ## Groups:   type, ref [9]
    #  type  ref   data             dr_mod mod_preds       
    #  <chr> <chr> <list>           <list> <list>          
    #1 A     ref3  <tibble [2 × 4]> <drc>  <df[,2] [2 × 2]>
    #2 A     ref1  <tibble [2 × 4]> <drc>  <df[,2] [2 × 2]>
    #3 A     ref2  <tibble [5 × 4]> <drc>  <df[,2] [5 × 2]>
    #4 B     ref2  <tibble [4 × 4]> <drc>  <df[,2] [4 × 2]>
    #5 B     ref1  <tibble [2 × 4]> <drc>  <df[,2] [2 × 2]>
    #6 B     ref3  <tibble [2 × 4]> <drc>  <df[,2] [2 × 2]>
    #7 C     ref3  <tibble [3 × 4]> <drc>  <df[,2] [3 × 2]>
    #8 C     ref2  <tibble [3 × 4]> <drc>  <df[,2] [3 × 2]>
    #9 C     ref1  <tibble [4 × 4]> <drc>  <df[,2] [4 × 2]>
    

    【讨论】:

    • 有道理!但是,当我仔细查看结果时,它们并不是我所期望的预测输出——看起来模型将每一列都应用为可能的输入。 dat_mods&lt;- dat_s%&gt;% group_by(type, ref)%&gt;% nest() %&gt;% mutate(dr_mod=map(data, dr_s))%&gt;% mutate(mod_preds=map2(data, dr_mod, ~modelr::add_predictions(data - .x, model = .y, type = "data.frame")))%&gt;% unnest(mod_preds) 难道是我必须指定别的东西为.x 和.y?
    • 这非常适合我提供的数据!但是,我显然使示例数据不能完全反映我的原始数据——info 在技术上应该是几个字符列。当时我认为差异并不重要(而且我认为列数并不重要),但显然调用as.data.frame(.x) 会影响livp1 的数据类型,因为我得到了很好的结果non-numeric argument to binary operator 错误。我正在编辑我原来的问题,如果您能再看一下,将不胜感激?
    • 问题在于子集的位置。如果我在add_predictions 内执行此操作(例如data=as.data.frame(.x)$livp1),它将使其成为一个向量,predict 以任何形式都不喜欢。除非我们能找到更优雅的解决方案,否则我认为我将不得不拆分并重新组合表格。再次抱歉弄错了原始问题代码!
    • 看看我的编辑,如果你使用[ 操作符对.x 进行子集化,它会返回一个tibble,然后可以将其转换为data.frame。 [ 和 $ 之间存在令人讨厌的细微差别。
    • 多田!我还刚刚发现,如果您在预测 mutate(dr_mod=map(data, dr_s))%&gt;% unnest(data)%&gt;% mutate(mod_preds=map2(livp1, dr_mod, ~add_predictions(data=as.data.frame(.x), model=.y))) 之前 unnest() 也可以使用,但我完全赞成节省代码行 :)
    猜你喜欢
    • 1970-01-01
    • 2020-08-24
    • 2016-02-16
    • 1970-01-01
    • 1970-01-01
    • 2023-04-03
    • 1970-01-01
    • 2018-03-02
    • 1970-01-01
    相关资源
    最近更新 更多