r语言 - 如何在组的类别上运行模型?

Question

我正在尝试使用垃圾过滤器包来模拟植物分解曲线。我的数据集有三列，时间，剩余质量和站点代码。我想将模型应用于网站代码的类别，并提取模型参数，但我有一些困难。下面是一个错误代码的例子:

library(litterfitter)
library(tidyverse)
decomp_test <- structure(list(site_code = c("CCPp1a", "CCPp1a", "CCPp1a", "CCPp1a", 
"CCPp1a", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1c", 
"CCPp1c", "CCPp1c", "CCPp1c", "CCPp1c", "CCPp1d", "CCPp1d", "CCPp1d", 
"CCPp1d", "CCPp1d", "CCPp1e", "CCPp1e", "CCPp1e", "CCPp1e", "CCPp1e", 
"CCPp1f", "CCPp1f", "CCPp1f", "CCPp1f", "CCPp1f"), days_between = c(0L, 
118L, 229L, 380L, 572L, 0L, 118L, 229L, 380L, 572L, 0L, 118L, 
229L, 380L, 572L, 0L, 118L, 229L, 380L, 572L, 0L, 118L, 229L, 
380L, 572L, 0L, 118L, 229L, 380L, 572L), mass_remaining = c(1, 
0.7587478816, 0.7366473295, 0.6038150404, 0.6339366063, 1, 0.7609346914, 
0.7487194938, 0.7336179508, 0.6595702348, 1, 0.777213425, 0.734006734, 
0.6963752241, 0.5827854154, 1, 0.7716566866, 0.7002094345, 0.6913555798, 
0.7519095328, 1, 0.7403565314, 0.6751289171, 0.6572164948, 0.620339994, 
1, 0.8126440236, 0.7272999401, 0.7223268259, 0.6805293006)), row.names = c(NA, 
-30L), class = "data.frame")
#Test data-frame with a small number of sites
discrete_paralell <-
decomp_test %>%
nest(-site_code) %>%
mutate(fit = map(decomp_test, ~ fit_litter(time=decomp_test$days_between ,mass.remaining= decomp_test$mass_remaining,
model='discrete.parallel',iters=1000)),
results = map(fit, glance)) %>% 
unnest(results)

错误:mutate()列fit有问题。ifit = map(...)。ifit的大小必须为6或1，而不是3。

#or
discrete_paralell <-
decomp_test %>%
nest(-site_code) %>%
mutate(fit = map(decomp_test, ~ fit_litter(time=decomp_test$days_between ,mass.remaining= decomp_test$mass_remaining,
model='discrete.parallel',iters=1000)),
coef = map_dbl(fit, "coef"),
actual = map_dbl(fit, "mass"),
preds = map_dbl(fit, "predicted"),
AIC = map_dbl(fit, "fitAIC"),
model = map_dbl(fit, "model"))

错误:mutate()列fit有问题。ifit = map(...)。ifit的大小必须为6或1，而不是3。

我知道不是所有的模型都适合，稍后我会检查所有的模型相对于其他模型的适合程度。

Answer 1

这个怎么样:

library(litterfitter)
library(tidyverse)
decomp_test <- structure(list(site_code = c("CCPp1a", "CCPp1a", "CCPp1a", "CCPp1a", 
"CCPp1a", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1b", "CCPp1c", 
"CCPp1c", "CCPp1c", "CCPp1c", "CCPp1c", "CCPp1d", "CCPp1d", "CCPp1d", 
"CCPp1d", "CCPp1d", "CCPp1e", "CCPp1e", "CCPp1e", "CCPp1e", "CCPp1e", 
"CCPp1f", "CCPp1f", "CCPp1f", "CCPp1f", "CCPp1f"), days_between = c(0L, 
                                          118L, 229L, 380L, 572L, 0L, 118L, 229L, 380L, 572L, 0L, 118L, 
                                          229L, 380L, 572L, 0L, 118L, 229L, 380L, 572L, 0L, 118L, 229L, 
                                          380L, 572L, 0L, 118L, 229L, 380L, 572L), mass_remaining = c(1, 
                                                                                                      0.7587478816, 0.7366473295, 0.6038150404, 0.6339366063, 1, 0.7609346914, 
                                                                                                      0.7487194938, 0.7336179508, 0.6595702348, 1, 0.777213425, 0.734006734, 
                                                                                                      0.6963752241, 0.5827854154, 1, 0.7716566866, 0.7002094345, 0.6913555798, 
                                                                                                      0.7519095328, 1, 0.7403565314, 0.6751289171, 0.6572164948, 0.620339994, 
                                                                                                      1, 0.8126440236, 0.7272999401, 0.7223268259, 0.6805293006)), row.names = c(NA, 
                                                                                                                                                                                 -30L), class = "data.frame")
#Test data-frame with a small number of sites
discrete_paralell <-
decomp_test %>%
group_by(site_code) %>% 
summarise(fit = list(fit_litter(time=days_between ,mass.remaining= mass_remaining,
model='discrete.parallel',iters=1000)$optimFit$par)) %>% 
unnest_wider(fit) %>% 
setNames(c("site_code", "par_1", "par_2", "par_3"))
#> Number of successful fits:  898  out of 1000
#> Number of successful fits:  912  out of 1000
#> Warning in fit_litter(time = days_between, mass.remaining = mass_remaining, :
#> one or more parameters fit on the boundary, check fit closely
#> Number of successful fits:  866  out of 1000
#> Number of successful fits:  981  out of 1000
#> Warning in fit_litter(time = days_between, mass.remaining = mass_remaining, :
#> one or more parameters fit on the boundary, check fit closely
#> Number of successful fits:  895  out of 1000
#> Warning in fit_litter(time = days_between, mass.remaining = mass_remaining, :
#> one or more parameters fit on the boundary, check fit closely
#> Number of successful fits:  907  out of 1000
#> Warning in fit_litter(time = days_between, mass.remaining = mass_remaining, :
#> one or more parameters fit on the boundary, check fit closely
#> New names:
#> • `` -> `...1`
#> • `` -> `...2`
#> • `` -> `...3`
#> New names:
#> New names:
#> New names:
#> New names:
#> New names:
#> • `` -> `...1`
#> • `` -> `...2`
#> • `` -> `...3`
discrete_paralell
#> # A tibble: 6 × 4
#>   site_code par_1    par_2    par_3
#>   <chr>     <dbl>    <dbl>    <dbl>
#> 1 CCPp1a    0.819 0.000512 0.596   
#> 2 CCPp1b    0.738 0.0001   0.0161  
#> 3 CCPp1c    0.120 0.257    0.000685
#> 4 CCPp1d    0.256 0.0187   0.0001  
#> 5 CCPp1e    0.332 0.0119   0.0001  
#> 6 CCPp1f    0.273 0.00954  0.0001

^{在2023-01-17由reprex包(v2.0.1)创建}