0

我正在尝试提取为我的工作设计的模型的参数，但在这样做时遇到了麻烦。我将使用 R 提供的虹膜数据作为示例数据集。

我拉入数据：

library(dplyr)
df <- iris

我设计了一个采用 Sepal Width 和 Length 的模型，并确定给定 Sepal Length 的最佳模型拟合来确定 Sepal Width （这可能是数据集的一个坏例子，我没有太注意起始参数，但我几乎可以肯定这不应该影响回答这个问题的难度）。为每个物种生产一个模型：

sepalmodel <- df %>% 
  group_by(Species) %>% 
  do(model = nls(Sepal.Width ~ a*exp(Sepal.Length*b), start = list(a = 0.5, b = 0.9), data = .)) %>% 
  ungroup()

这会生成一个数据框，其中一列是 Species，另一列是模型。在模型列中，每个模型在每一行中都由一个大列表表示，我无法从中提取其基本数据。我想直接使用参数，作为使用这些模型表示和生成数据的更干净和可重复的方式。如何从这些模型中提取 a 和 b 参数？有没有这样做的功能，还是我需要深入研究每个模型的代码？此外，该列表中是否有任何数据将其与特定物种相关联，或者它是否仅与在同一行中发现的物种直接相关？

4

3 回答 3

2

一种选择是使用broom::tidywith tidyr::unnest。

library(tidyverse)

iris %>% 
  group_by(Species) %>% 
  do(model = nls(Sepal.Width ~ a*exp(Sepal.Length*b), start = list(a = 0.5, b = 0.9), data = .)) %>% 
  ungroup() %>%
  mutate(tidy_nls = lapply(model, broom::tidy)) %>%
  unnest(tidy_nls)

#------
# A tibble: 6 x 7
  Species    model  term  estimate std.error statistic  p.value
  <fct>      <list> <chr>    <dbl>     <dbl>     <dbl>    <dbl>
1 setosa     <nls>  a       1.07      0.162       6.58 3.23e- 8
2 setosa     <nls>  b       0.232     0.0299      7.76 5.12e-10
3 versicolor <nls>  a       1.41      0.228       6.18 1.31e- 7
4 versicolor <nls>  b       0.113     0.0269      4.21 1.11e- 4
5 virginica  <nls>  a       1.80      0.261       6.87 1.17e- 8
6 virginica  <nls>  b       0.0764    0.0218      3.51 9.97e- 4
于 2021-01-27T22:06:29.417 回答
2

这里有两种方法。

1) nlme这个包是 R 自带的，所以不需要安装。可以像往常一样指定公式，后跟竖线和允许一次指定所有模型的分组变量。

library(nlme)
fm <- nlsList(Sepal.Width ~ a * exp(Sepal.Length*b) | Species, data = iris, 
  start = list(a = 0.5, b = 0.5))
co <- coef(summary(fm))
co

给出这个数组，例如，co["setosa",,], co[, "Estimate", ], co[,,"a"] 给出 setosa、所有估计值和所有“a”系数的矩阵. 整个数组是：

    , , a

               Estimate Std. Error  t value     Pr(>|t|)
    setosa     1.068204  0.1729358 6.176882 3.226364e-08
    versicolor 1.412344  0.2302122 6.134967 1.314458e-07
    virginica  1.795394  0.2453925 7.316418 1.170751e-08
    
    , , b
    
                 Estimate Std. Error  t value     Pr(>|t|)
    setosa     0.23226126 0.03189911 7.281121 5.122595e-10
    versicolor 0.11319887 0.02708865 4.178831 1.107640e-04
    virginica  0.07643349 0.02046418 3.734989 9.965872e-04

ftable可用于以各种 2d 方式查看此数组。此链接还提供了一个函数 ftable2df，可以将 ftable 转换为数据框。使用 ftable 的示例：

ftable(co, row.vars = 3:2)

给予：

                    setosa   versicolor    virginica
                                                    
a Estimate    1.068204e+00 1.412344e+00 1.795394e+00
  Std. Error  1.729358e-01 2.302122e-01 2.453925e-01
  t value     6.176882e+00 6.134967e+00 7.316418e+00
  Pr(>|t|)    3.226364e-08 1.314458e-07 1.170751e-08
b Estimate    2.322613e-01 1.131989e-01 7.643349e-02
  Std. Error  3.189911e-02 2.708865e-02 2.046418e-02
  t value     7.281121e+00 4.178831e+00 3.734989e+00
  Pr(>|t|)    5.122595e-10 1.107640e-04 9.965872e-04

作为另一个例子

ftable(co, row.vars = c(1, 3))

给予：

                  Estimate   Std. Error      t value     Pr(>|t|)
                                                                 
setosa     a  1.068204e+00 1.729358e-01 6.176882e+00 3.226364e-08
           b  2.322613e-01 3.189911e-02 7.281121e+00 5.122595e-10
versicolor a  1.412344e+00 2.302122e-01 6.134967e+00 1.314458e-07
           b  1.131989e-01 2.708865e-02 4.178831e+00 1.107640e-04
virginica  a  1.795394e+00 2.453925e-01 7.316418e+00 1.170751e-08
           b  7.643349e-02 2.046418e-02 3.734989e+00 9.965872e-04

也是S3 对象fm，summary(fm)即具有类变量的列表。列表的组成部分可用，可通过以下方式查看：

str(fm)
str(summary(fm))

fm具有类向量c("nlsList", "lmList")，因此您可以找到所有可用于操作的方法，fm如下所示：

methods(class = "nlsList")
methods(class = "lmList")

事实上fm，分组变量的每一级都有一个组件，每个组件都是一个nls对象，我们可以获得可以应用于nls对象的方法，如下所示：

methods(class = "nls")

同样summary(fm)有类向量c("summary.nlsList", "summary.lmList")，因此我们可以使用以下方法找到可以应用于它的方法：

methods(class = "summary.nlsList")
methods(class = "summary.lmList")

2) nls使用plain nls，我们可以通过将分组变量指定为每个参数的索引来创建单个模型。起始值是一个列表，其中包含每个参数的向量分量，每个参数都有一个分组向量的每个级别的条目。Species 有 3 个级别，因此我们有 3 个起始值a和 3个起始值b。

st <- list(a = rep(0.5, 3), b = rep(0.5, 3))
fm <- nls(Sepal.Width ~ a[Species] * exp(Sepal.Length*b[Species]), iris, start = st)
co <- coef(summary(fm))
co

给出这个矩阵：

     Estimate Std. Error  t value     Pr(>|t|)
a1 1.06820417 0.17293582 6.176882 6.329948e-09
a2 1.41233648 0.23021095 6.134967 7.804013e-09
a3 1.79539324 0.24539239 7.316418 1.638244e-11
b1 0.23226125 0.03189911 7.281121 1.983862e-11
b2 0.11319979 0.02708865 4.178864 5.058757e-05
b3 0.07643355 0.02046418 3.734992 2.697815e-04

您可以考虑将行名更改为：

rownames(co) <- c(t(outer(c("a", "b"), levels(iris$Species), paste, sep = ".")))
co

给予：

               Estimate Std. Error  t value     Pr(>|t|)
a.setosa     1.06820417 0.17293582 6.176882 6.329948e-09
a.versicolor 1.41233648 0.23021095 6.134967 7.804013e-09
a.virginica  1.79539324 0.24539239 7.316418 1.638244e-11
b.setosa     0.23226125 0.03189911 7.281121 1.983862e-11
b.versicolor 0.11319979 0.02708865 4.178864 5.058757e-05
b.virginica  0.07643355 0.02046418 3.734992 2.697815e-04

如（1）中，我们可以检查 fm using 的组件，str(fm)并可以找到可以作用于fmusing的方法methods(class = "nls")。

于 2021-01-27T22:36:17.943 回答
1

更一般地，如果您想从模型中提取参数，您可以采用与从数据框中选择列类似的方式进行操作。例如，如果nls您可以键入?nls以检查可以从模型返回的内容（查看值）。例如，nls$model获取模型名称，nls$weights获取权重等等。值得注意的是，相同的方法可用于相关性。您可以键入?cor.test（并在值下查看）以查看可以提取的参数。cor.test(x,y)$estimate

于 2021-01-28T01:48:10.720 回答