为了简化再现性,我使用goats
数据集来自ResourceSelection数据集,该数据集具有用于山羊的二手(STATUS == 1
)和'可用'(STATUS == 0
)GPS位置的空间数据。 ID
用于个体(n = 10)和ELEVATION
,...,TASP
是这些点的属性。
library(ResourceSelection)
head(goats)
STATUS ID ELEVATION SLOPE ET ASPECT HLI TASP
1 1 1 651 38.5216 35.3553 243.1131 0.9175926 0.9468804
2 1 1 660 39.6927 70.7107 270.0000 0.8840338 0.6986293
3 1 1 316 20.5477 50.0000 279.2110 0.7131423 0.5749115
4 1 1 334 34.0783 35.3553 266.1859 0.8643775 0.7447368
5 1 1 454 41.6187 25.0000 258.3106 0.9349181 0.8292587
6 1 1 343 28.4694 103.0776 237.0426 0.8254866 0.9756112
我想为glm
(下面创建)分组的每个人拟合Season
,但模型结构根据季节而有所不同。我一直在使用一些有用的SO posts和其他resources但他们都适合每个组的单一模型,因为我想为每个Season
分组适合不同的模型。
#Add a new `Season` field
library(tidyverse)
goats <- goats %>%
mutate(Season = if_else(ID %in% 1:3, "Summer",
if_else(ID %in% 4:7, "Winter", "Fall")))
下面我创建一个模型构建函数,并使用if else
指定每个季节特定的模型。
SeasonalMods <- function(df) {
#Models for Archery
if(Season == "Summer") {
glm(STATUS ~ SLOPE + I(SLOPE^2), data = df)
#Models for Winter
} else if (Season == "Winter") {
glm(STATUS ~ SLOPE + ASPECT + TASP, data = df)
#Models for Fall
} else if (Season == "Fall") {
glm(STATUS ~ ELEVATION + SLOPE + ET + ASPECT + HLI + TASP, data = df)}
}
然后我尝试将函数映射到分组数据并创建新的列表列,如下所示。
ModelFits <- goats %>%
group_by(Season, ID) %>%
nest() %>%
mutate(fits = map(data, SeasonalMods),
tidied = map(fits, tidy),
glanced = map(fits, glance),
augmented = map(fits, augment))
这会产生以下错误:
Error in mutate_impl(.data, dots) :
Evaluation error: object 'Season' not found
我不确定如何在Season
函数中正确指定SeasonalMods
,以便map()
可以解释它。
我试图在df$
和Season
语句中在if
前添加else if
,但这也会产生错误。
在我建议使用modelr
后,我没有它就继续前进。就像我上面说的那样,你的函数SeasonalMods
不知道Season
是数据框中的一个列,它作为一个参数,所以你会得到它未定义的错误。换句话法的一种方法是在函数中添加第二个参数以获取季节。由于您正在嵌套数据,现在很容易将数据和季节传递给建模函数。我正在使用map2
,因为列data
和Season
长度相同。
library(tidyverse)
library(broom)
这个函数的所有内部都是相同的 - 我只添加了第二个参数。
SeasonalMods <- function(df, Season) {
...
}
为了说明您使用的broom
函数,我添加了tidied
列,并保存此数据框:
models <- goats %>%
group_by(Season, ID) %>%
nest() %>%
mutate(fits = map2(data, Season, ~SeasonalMods(.x, .y))) %>%
mutate(tidied = map(fits, tidy))
head(models)
#> # A tibble: 6 x 5
#> Season ID data fits tidied
#> <chr> <int> <list> <list> <list>
#> 1 Summer 1 <tibble [2,106 × 7]> <S3: glm> <tibble [3 × 5]>
#> 2 Summer 2 <tibble [1,668 × 7]> <S3: glm> <tibble [3 × 5]>
#> 3 Summer 3 <tibble [1,539 × 7]> <S3: glm> <tibble [3 × 5]>
#> 4 Winter 4 <tibble [951 × 7]> <S3: glm> <tibble [4 × 5]>
#> 5 Winter 5 <tibble [1,908 × 7]> <S3: glm> <tibble [4 × 5]>
#> 6 Winter 6 <tibble [2,184 × 7]> <S3: glm> <tibble [4 × 5]>
只是为了检查模型是否有不同季节的不同公式:
models$fits[[1]]
#>
#> Call: glm(formula = STATUS ~ SLOPE + I(SLOPE^2), data = df)
#>
#> Coefficients:
#> (Intercept) SLOPE I(SLOPE^2)
#> -0.042618 -0.000989 0.000375
#>
#> Degrees of Freedom: 2105 Total (i.e. Null); 2103 Residual
#> Null Deviance: 468
#> Residual Deviance: 337.2 AIC: 2127
models$fits[[6]]
#>
#> Call: glm(formula = STATUS ~ SLOPE + ASPECT + TASP, data = df)
#>
#> Coefficients:
#> (Intercept) SLOPE ASPECT TASP
#> 0.024625 0.017838 -0.001768 0.215217
#>
#> Degrees of Freedom: 2183 Total (i.e. Null); 2180 Residual
#> Null Deviance: 485.3
#> Residual Deviance: 385.7 AIC: 2421