如何从 R 中列出的模型中提取参数？
How do I extract parameters from listed models in R?
我正在尝试提取为我的工作设计的模型的参数，但在这样做时遇到了问题。我将使用 R 提供的虹膜数据作为示例数据集。
我拉入数据：
library(dplyr)
df <- iris
我设计了一个采用萼片宽度和长度的模型，并确定给定萼片长度的最佳模型拟合是确定萼片宽度（这可能是一个糟糕的数据集示例，我没有付出太多注意起始参数，但我几乎可以肯定这不会影响回答这个问题的难度）。每个物种生产一个模型：
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()
这会生成一个数据框，其中一列是物种，另一列是模型。在模型列中，每个模型在每一行中都由一个大列表表示，我很难从中提取其基本数据。我想直接使用参数，作为使用这些模型表示和生成数据的更简洁和可重现的方式。如何从这些模型中提取 a 和 b 参数？是否有这样做的功能，或者我是否需要深入研究每个模型的代码？此外，该列表中是否内置任何数据将其与特定物种相关联，或者它是否仅与同一行中发现的物种直接相关？
一种选择是将 broom::tidy 与 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
这里有两种方法。
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 可用于以各种二维方式查看此数组。此 link 还提供了一个函数 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
fm 和 summary(fm) 也是 S3 对象，即带有 class 变量的列表。列表的组成部分可用，可通过以下方式查看：
str(fm)
str(summary(fm))
fm 有 class 向量 c("nlsList", "lmList") 所以你可以找到所有可用于作用于 fm 的方法，如下所示：
methods(class = "nlsList")
methods(class = "lmList")
事实上 fm 对分组变量的每个级别都有一个组件，每个组件都是一个 nls 对象，我们可以得到可以应用于 [=31= 的方法] 这样的对象：
methods(class = "nls")
类似地，summary(fm) 有 class 向量 c("summary.nlsList", "summary.lmList")，所以我们可以找到可以应用于它的方法：
methods(class = "summary.nlsList")
methods(class = "summary.lmList")
2) nls 使用普通 nls 我们可以通过将分组变量指定为每个参数的索引来创建单个模型。起始值是一个列表，其中包含每个参数的矢量分量，每个参数都有一个分组矢量的每个级别的条目。物种有 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) 一样，我们可以使用 str(fm) 检查 fm 的组件，并可以使用 methods(class = "nls") 找到可以作用于 fm 的方法。
更一般地说，如果您想从模型中提取参数，可以采用与从数据框中选择列类似的方式进行。例如，在 nls 的情况下，您可以键入 ?nls 来检查您可以从模型中 return 得到什么（查看值）。比如nls$model获取模型名称，nls$weights获取权重等。值得注意的是，相同的方法可用于相关性。您可以键入 ?cor.test（并查看值）以查看可以提取的参数。 cor.test(x,y)$estimate
我正在尝试提取为我的工作设计的模型的参数，但在这样做时遇到了问题。我将使用 R 提供的虹膜数据作为示例数据集。
我拉入数据：
library(dplyr)
df <- iris
我设计了一个采用萼片宽度和长度的模型，并确定给定萼片长度的最佳模型拟合是确定萼片宽度（这可能是一个糟糕的数据集示例，我没有付出太多注意起始参数，但我几乎可以肯定这不会影响回答这个问题的难度）。每个物种生产一个模型：
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()
这会生成一个数据框，其中一列是物种，另一列是模型。在模型列中，每个模型在每一行中都由一个大列表表示，我很难从中提取其基本数据。我想直接使用参数，作为使用这些模型表示和生成数据的更简洁和可重现的方式。如何从这些模型中提取 a 和 b 参数？是否有这样做的功能，或者我是否需要深入研究每个模型的代码？此外，该列表中是否内置任何数据将其与特定物种相关联，或者它是否仅与同一行中发现的物种直接相关？
一种选择是将 broom::tidy 与 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
这里有两种方法。
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 可用于以各种二维方式查看此数组。此 link 还提供了一个函数 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
fm 和 summary(fm) 也是 S3 对象，即带有 class 变量的列表。列表的组成部分可用，可通过以下方式查看：
str(fm)
str(summary(fm))
fm 有 class 向量 c("nlsList", "lmList") 所以你可以找到所有可用于作用于 fm 的方法，如下所示：
methods(class = "nlsList")
methods(class = "lmList")
事实上 fm 对分组变量的每个级别都有一个组件，每个组件都是一个 nls 对象，我们可以得到可以应用于 [=31= 的方法] 这样的对象：
methods(class = "nls")
类似地，summary(fm) 有 class 向量 c("summary.nlsList", "summary.lmList")，所以我们可以找到可以应用于它的方法：
methods(class = "summary.nlsList")
methods(class = "summary.lmList")
2) nls 使用普通 nls 我们可以通过将分组变量指定为每个参数的索引来创建单个模型。起始值是一个列表，其中包含每个参数的矢量分量，每个参数都有一个分组矢量的每个级别的条目。物种有 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) 一样，我们可以使用 str(fm) 检查 fm 的组件，并可以使用 methods(class = "nls") 找到可以作用于 fm 的方法。
更一般地说，如果您想从模型中提取参数，可以采用与从数据框中选择列类似的方式进行。例如，在 nls 的情况下，您可以键入 ?nls 来检查您可以从模型中 return 得到什么（查看值）。比如nls$model获取模型名称，nls$weights获取权重等。值得注意的是，相同的方法可用于相关性。您可以键入 ?cor.test（并查看值）以查看可以提取的参数。 cor.test(x,y)$estimate