0

我正在尝试使用drc包中的 Gompertz 函数来拟合纵向数据的增长模型。为了同时在多个国家/地区生成模型，我嵌套了数据并使用了tidyverse 的map函数。使用purrr::map的建模步骤运行良好，但是当我尝试使用 broom 中的glance选项获取模型摘要时，我收到一个错误，即drc 类的对象不存在glance 方法。但是，根据此网页，glance 似乎对 drc 对象有效：https ://rdrr.io/github/tidyverse/broom/man/glance.drc.html

下面是一个示例，包括涉及 2 个国家和每个国家 6 个时间点的玩具数据集。

#数据集测试

iso_code    total_cases num_days  
IND 1   1  
IND 3   10  
IND 3   20  
IND 3   30  
IND 165 50  
IND 979 60  
SGP 3   1  
SGP 18  10  
SGP 47  20  
SGP 86  30  
SGP 187 50  
SGP 455 60

#安装包tidyverse和drc

> library(tidyverse)
> library(drc)

> data = read.delim("test_data.txt", sep='\t', header=T)
> dim(data)

[1] 12  3

> str(data)

'data.frame':   12 obs. of  3 variables:
 $ iso_code   : chr  "IND" "IND" "IND" "IND" ...
 $ total_cases: int  1 3 3 3 165 979 3 18 47 86 ...
 $ num_days   : int  1 10 20 30 50 60 1 10 20 30 ...

#按国家/地区嵌套数据

> by_country = data %>% group_by(iso_code)%>% nest()

#在嵌套数据上运行 drc gompertz 模型

> by_country_g = by_country %>% mutate(model= purrr::map(data,~  drm((total_cases)/max(total_cases) ~ num_days, fct=G.4(), data=.)))

#检查模型内容

> by_country_g$model

[[1]]

A 'drc' model.

Call:
drm(formula = (total_cases)/max(total_cases) ~ num_days, data = .,     fct = G.4())

Coefficients:
b:(Intercept)  c:(Intercept)  d:(Intercept)  e:(Intercept)  
    -0.164861       0.002555       1.531310      54.838386  


[[2]]

A 'drc' model.

Call:
drm(formula = (total_cases)/max(total_cases) ~ num_days, data = ., fct = G.4())

Coefficients:
b:(Intercept)  c:(Intercept)  d:(Intercept)  e:(Intercept)  
      -0.1257         0.0845         1.4638        52.9056

#通过一目了然提取模型摘要

> summary_g = by_country_g %>% mutate(glance = purrr::map(model, broom::glance))%>% unnest(glance)


Error: Problem with mutate() input glance.  
x No glance method for objects of class drc  
ℹ Input glance is purrr::map(model, broom::glance).  
ℹ The error occured in group 1: iso_code = "IND".
Run rlang::last_error() to see where the error occurred.
4

1 回答 1

0 
library(tidyverse)
library(drc)

data %>% 
  nest(data = -iso_code) %>%
  mutate(model= map(data,~  drm((total_cases)/max(total_cases) ~ num_days, 
                    fct=G.4(), data=.)), 
         glance = map(model, broom::glance)) %>%
  unnest(glance)

#   iso_code data             model     AIC    BIC logLik    df.residual
#  <chr>    <list>           <list>  <dbl>  <dbl> <logLik>        <int>
#1 IND      <tibble [6 × 2]> <drc>  -59.8  -60.8  34.884966           2
#2 SGP      <tibble [6 × 2]> <drc>   -7.39  -8.43  8.694533           2
于 2021-05-30T06:23:22.000 回答