问题描述
@H_404_0@我正在计算由
我想以百分比形式测试我的数据，因此首先我将其减去1，然后除以4，以进行转换：
@H_404_0@取消嵌套
@H_404_0@ 演示如何获取没有管道或数据帧的信息 @H_404_0@以下方法依赖于为过程分配多个步骤的变量。为了演示起见，它显示了如何运行模型并仅获取一种水果的
glm生成的模型的输出。模型输出存储在嵌套的小节中。我想通过从type =“ link”到反向链接（使用$family$linkinv）的转换来计算置信区间。但是，我无法使其与dplyr::mutate一起嵌套嵌套使用，因为从$family$linkinv的模型对象中拉出model$family$linkinv(x)的方式似乎不符合预期嵌套格式。
背景
@H_404_0@当前问题基于我发布的previous question（并选择了答案），有关使用线性模型通过不同预测变量测试喜欢水果的水平。我进行了一项研究，以确定哪种水果更受欢迎：芒果，香蕉或苹果。为此，我继续进行随机抽样100人。我要求他们以1-5的等级来评价每种水果的喜好程度。 @H_404_0@尽管上一个问题与lm有关，但在这里我尝试使用准二项式glm。问题是我想获得置信区间，但是我的方法（glm %>% predict）在“链接空间”中输出SE，因此我必须经过转换过程（detailed in this SO answer）才能获得所需的结果。
数据
library(tidyverse)
library(magrittr)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,i_love_apple = sample(c(1:5),100,replace = TRUE),i_love_banana = sample(c(1:5),i_love_mango = sample(c(1:5),age = sample(c(20:70),is_male = sample(c(0,1),prob = c(0.2,0.8),education_level = sample(c(1:4),is_colorblinded = sample(c(0,replace = TRUE)
)
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <int> <int> <int> <int> <dbl> <int> <dbl>
## 1 1 3 5 2 50 1 2 0
## 2 2 3 3 1 49 1 1 0
## 3 3 2 1 5 70 1 1 1
## 4 4 2 2 5 41 1 3 1
## 5 5 3 1 1 49 1 4 0
## 6 6 5 2 1 29 0 1 0
## 7 7 4 5 5 35 1 3 0
## 8 8 1 3 5 24 0 3 0
## 9 9 2 4 2 55 1 2 0
## 10 10 3 4 2 69 1 4 0
## # ... with 90 more rows
我想以百分比形式测试我的数据，因此首先我将其减去1，然后除以4，以进行转换：
fruit_liking_df %<>%
mutate_at(vars(starts_with("i_love_")),~ subtract(.,1) %>% divide_by(.,4))
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <dbl> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 1 0.5 1 0.25 50 1 2 0
## 2 2 0.5 0.5 0 49 1 1 0
## 3 3 0.25 0 1 70 1 1 1
## 4 4 0.25 0.25 1 41 1 3 1
## 5 5 0.5 0 0 49 1 4 0
## 6 6 1 0.25 0 29 0 1 0
## 7 7 0.75 1 1 35 1 3 0
## 8 8 0 0.5 1 24 0 3 0
## 9 9 0.25 0.75 0.25 55 1 2 0
## 10 10 0.5 0.75 0.25 69 1 4 0
## # ... with 90 more rows
现在，我使用管道为每个水果运行glm模型，在链接空间中获取SE，并将SE转换为CI
## will be needed later
my_new_data_for_pred <- expand_grid(
age = 45,is_male = .5,education_level = 2.5,is_colorblinded = 0.5
)
## will be needed later
critval <- 1.96
model_fits_grouped <-
fruit_liking_df %>%
pivot_longer(starts_with("i_love"),values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,~ glm(
data = .x,fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,~ bind_cols(my_new_data_for_pred,as.data.frame(
predict(
newdata = my_new_data_for_pred,.x,type = "link",interval = "confidence",level = 0.95,se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,lower_ci_link = fit - critval * se.fit,upper_ci_link = fit + critval * se.fit
)
))
> model_fits_grouped
## # A tibble: 3 x 4
## # Groups: name [3]
## name data model_fit predicted_values
## <chr> <list> <list> <list>
## 1 i_love_apple <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 2 i_love_banana <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 3 i_love_mango <tibble [100 x 6]> <glm> <tibble [1 x 10]>
@H_404_0@取消嵌套fruit_liking_df %<>%
mutate_at(vars(starts_with("i_love_")),~ subtract(.,1) %>% divide_by(.,4))
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <dbl> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 1 0.5 1 0.25 50 1 2 0
## 2 2 0.5 0.5 0 49 1 1 0
## 3 3 0.25 0 1 70 1 1 1
## 4 4 0.25 0.25 1 41 1 3 1
## 5 5 0.5 0 0 49 1 4 0
## 6 6 1 0.25 0 29 0 1 0
## 7 7 0.75 1 1 35 1 3 0
## 8 8 0 0.5 1 24 0 3 0
## 9 9 0.25 0.75 0.25 55 1 2 0
## 10 10 0.5 0.75 0.25 69 1 4 0
## # ... with 90 more rows
## will be needed later
my_new_data_for_pred <- expand_grid(
age = 45,is_male = .5,education_level = 2.5,is_colorblinded = 0.5
)
## will be needed later
critval <- 1.96
model_fits_grouped <-
fruit_liking_df %>%
pivot_longer(starts_with("i_love"),values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,~ glm(
data = .x,fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,~ bind_cols(my_new_data_for_pred,as.data.frame(
predict(
newdata = my_new_data_for_pred,.x,type = "link",interval = "confidence",level = 0.95,se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,lower_ci_link = fit - critval * se.fit,upper_ci_link = fit + critval * se.fit
)
))
> model_fits_grouped
## # A tibble: 3 x 4
## # Groups: name [3]
## name data model_fit predicted_values
## <chr> <list> <list> <list>
## 1 i_love_apple <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 2 i_love_banana <tibble [100 x 6]> <glm> <tibble [1 x 10]>
## 3 i_love_mango <tibble [100 x 6]> <glm> <tibble [1 x 10]>
predicted_values会得到：
> model_fits_grouped %>% unnest(predicted_values)
## # A tibble: 3 x 13
## # Groups: name [3]
## name data model_fit age is_male education_level is_colorblinded fit se.fit residual.scale estimate lower_ci_link upper_ci_link
## <chr> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 i_love_apple <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 0.0843 0.261 0.709 0.0843 -0.427 0.595
## 2 i_love_banana <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 -0.0718 0.286 0.781 -0.0718 -0.633 0.489
## 3 i_love_mango <tibble [100 x 6~ <glm> 45 0.5 2.5 0.5 -0.140 0.279 0.762 -0.140 -0.687 0.407
predicted_values内的{em> {1>}和lower_ci_link的反向链接转换中再增加两列，但这失败了
upper_ci_link@H_404_0@错误：@H_404_0@我认为问题是我正在尝试对model_fits_grouped <- fruit_liking_df %>% pivot_longer(starts_with("i_love"),upper_ci_link = fit + critval * se.fit ) %>% ######################### this addition fails ########################### mutate( lower_ci_inverse_link = model_fit$family$linkinv(lower_ci_link),upper_ci_inverse_link = model_fit$family$linkinv(upper_ci_link) ) ######################################################################### ))输入mutate()出现问题。 x问题 使用predicted_values输入mutate()。 x尝试申请 非功能我输入lower_ci_inverse_link是lower_ci_inverse_link。我的错误发生在行中
- i输入
model_fit$family$linkinv(lower_ci_link)是predicted_values。我的错误发生在第1行。
map(...)中的新列进行突变，但是使用predicted_values是指model_fit$family$linkinv(lower_ci_link)，它在嵌套小标题中处于较高级别。
底线问题
@H_404_0@如何使用model_fit和predicted_values来对 model_fit$family$linkinv(lower_ci_link)中的内的反向链接列进行突变以最终获得（一直滚动到最右边的两个列）：
model_fit$family$linkinv(upper_ci_link)附录
@H_404_0@ 演示如何获取没有管道或数据帧的信息 @H_404_0@以下方法依赖于为过程分配多个步骤的变量。为了演示起见，它显示了如何运行模型并仅获取一种水果的
> model_fits_grouped %>% unnest(predicted_values)
## # A tibble: 3 x 15
## # Groups: name [3]
## name data model_fit age is_male education_level is_colorblinded fit se.fit residual.scale estimate lower_ci_link upper_ci_link lower_ci_inverse_link_*DEMO* upper_ci_inverse_link_*DEMO*
## <chr> <list> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 i_love_apple <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.521 0.0632 0.349 0.521 0.397 0.645 0.111 0.111
## 2 i_love_banana <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.482 0.0701 0.387 0.482 0.345 0.620 0.222 0.222
## 3 i_love_mango <tibble [100 x 6]> <glm> 45 0.5 2.5 0.5 0.465 0.0683 0.377 0.465 0.331 0.599 0.333 0.333
。
数据
@H_404_0@像以前一样，它是经过$family$linkinv的算术转换为小数之后，因此：
fruit_liking_df型号
@H_404_0@我将仅关注> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <dbl> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 1 0.5 1 0.25 50 1 2 0
## 2 2 0.5 0.5 0 49 1 1 0
## 3 3 0.25 0 1 70 1 1 1
## 4 4 0.25 0.25 1 41 1 3 1
## 5 5 0.5 0 0 49 1 4 0
## 6 6 1 0.25 0 29 0 1 0
## 7 7 0.75 1 1 35 1 3 0
## 8 8 0 0.5 1 24 0 3 0
## 9 9 0.25 0.75 0.25 55 1 2 0
## 10 10 0.5 0.75 0.25 69 1 4 0
## # ... with 90 more rows
列数据，并对其运行i_love_apple。
glm预测
@H_404_0@现在，我使用my_model <-
glm(
i_love_apple ~
I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - 0.5) +
I(education_level - 2) +
I(is_colorblinded - 0.5),family = quasibinomial,data = fruit_liking_df
)
的预测数据在predict()上运行my_model：
my_new_data_for_predprediction_link_type <-
predict(object = my_model,newdata = my_new_data_for_pred,## <------------ type = "link" is crucial to note
interval = "confidence",se.fit = TRUE)
> prediction_link_type
## $fit
## 1
## 0.08427577
## $se.fit
## [1] 0.2606326
## $residual.scale
## [1] 0.7090294
（已分配给prediction_link_type），将我从critval获得的SE度量转换为置信区间（CI）。我分配了两个单独的向量：一个具有上限CI，另一个具有下限CI：
1.96lower_ci_link <- prediction_link_type$fit - (critval * prediction_link_type$se.fit)
upper_ci_link <- prediction_link_type$fit + (critval * prediction_link_type$se.fit)
使用了predict()）。要将CI值从“链接”转换回去，我使用了反向链接功能：
type = "link"解决方法
要引用在map中传递的数据，您需要使用.x。尝试以下答案。
library(tidyverse)
result <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"),values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(
data,~ glm(
data = .x,fruit ~ I(age - 45) +
I((age - 45) ^ 2) +
I(is_male - .5) +
I(education_level - 2) +
is_colorblinded,family = quasibinomial
)
)) %>%
mutate(predicted_values = map(
model_fit,~ bind_cols(my_new_data_for_pred,as.data.frame(
predict(
newdata = my_new_data_for_pred,.x,type = "link",interval = "confidence",level = 0.95,se.fit = T
)
)) %>%
rowwise() %>%
mutate(
estimate = fit,lower_ci_link = fit - critval * se.fit,upper_ci_link = fit + critval * se.fit,lower_ci_inverse_link = .x$family$linkinv(lower_ci_link),upper_ci_inverse_link = .x$family$linkinv(upper_ci_link)
)))
result看起来像：
result
# name data model_fit predicted_values
# <chr> <list> <list> <list>
#1 i_love_apple <tibble [100 × 6]> <glm> <tibble [1 × 12]>
#2 i_love_banana <tibble [100 × 6]> <glm> <tibble [1 × 12]>
#3 i_love_mango <tibble [100 × 6]> <glm> <tibble [1 × 12]>
要获取所有值作为单独的列，可以使用unnest_wider：
result %>% unnest_wider(predicted_values)
# name data model_fit age is_male education_level is_colorblinded fit se.fit
# <chr> <lis> <list> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#1 i_lo… <tib… <glm> 45 0.5 2.5 0.5 0.0843 0.261
#2 i_lo… <tib… <glm> 45 0.5 2.5 0.5 -0.0718 0.286
#3 i_lo… <tib… <glm> 45 0.5 2.5 0.5 -0.140 0.279
# … with 6 more variables: residual.scale <dbl>,estimate <dbl>,lower_ci_link <dbl>,# upper_ci_link <dbl>,lower_ci_inverse_link <dbl>,upper_ci_inverse_link <dbl>