1

データ1~3に対して複数の回帰分析を行い、その結果をこのサイトにあるような図にまとめたいと考えています。

回帰分析の結果は以下のようにまとまっています。

> lmRes_set
Source: local data frame [5 x 3]
Groups: <by row>

# A tibble: 5 × 3
   iter            lmResRand            lmResBias
* <dbl>               <list>               <list>
1     1 <data.frame [7 × 5]> <data.frame [7 × 5]>
2     2 <data.frame [6 × 5]> <data.frame [6 × 5]>
3     3 <data.frame [6 × 5]> <data.frame [6 × 5]>

一つ一つの結果は、

> lmRes_set[[2]]
[[1]]
              term    estimate  std.error   statistic    p.value
1      (Intercept) 44.08651614 12.7941054  3.44584594 0.04106068
2   treatment_rand -1.00512172 10.4939360 -0.09578119 0.92973389
3              age  0.06060037  0.1618392  0.37444793 0.73297251
4 factor(partyID)4  3.32417091 10.8821091  0.30547120 0.77997631
5 factor(partyID)7 -1.10496660 11.2153538 -0.09852267 0.92773124
6 factor(partyID)8 23.71469512 10.9822559  2.15936465 0.11964723
7           gender  9.75602334  7.2412885  1.34727726 0.27061498

のようになっていて、このうちtreatment_randだけを使い、先ほどのページの図の(Intercept)のように、一つの変数に対して3つのデータセットを結果を色を変えつつ表示したいと考えています。

この場合、結果の入ったlmRes_Setをどのような形でggolot2に流すべきなのでしょうか。いくつかためしましたが、上手くいきませんでした。

以下私が使いたいデータに近い擬似データを作るコードと、試したコードになっています

library(tidyverse);library(broom)
# Data Creation
population_num <- 50
population <- tibble(
    gender = as.numeric(rbinom(population_num, 1, 0.5)),
    age=rnorm(population_num, mean=50, sd=20),
    score=rnorm(population_num, mean=80, sd=30),
    setid=sample(c(1,2,3), size=population_num, replace=T)
    )
# Analysis
lmRes_set <- population %>%
    group_by(setid) %>%
    do(model1=tidy(lm(score~age, data=.)),
         model2=tidy(lm(score~age+gender, data=.)))
# Make Plot for model 1 # Error!!!
temp <- lmRes_set %>%
    group_by(setid) %>%
    select(model1)

よろしくお願い致します。
(こちらの質問は英語版での質問に反応が悪かったため、日本語で再質問したものになります。どちらかで回答頂けた場合は質問者が責任を持ってもう片方に翻訳致します。)

2 件の回答 2

1

lmRes_setはロング形式ではない、処理し難いデータとなっています。まずはこれをロング形式化し、以降の処理を行っていく形となります。

temp <- population %>%
  group_by(setid) %>%
  do(model1 = tidy(lm(score ~ age, data = .)),
     model2 = tidy(lm(score ~ age + gender, data = .))) %>%   ## ここまで質問者様コード
  gather(model_name, model, -setid) %>%                        ## ロング形式化
  unnest() %>%                                                 ## 展開
  filter(term == "age")                                        ## 必用な項のみ選択


  ## 以降は質問者様リンク先サイトのコードとほぼ同様です (一部割愛)
interval1 <- -qnorm((1-0.9)/2)

ggplot(temp, aes(colour = as.factor(setid))) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = model_name, ymin = estimate - std.error*interval1,
                     ymax = estimate + std.error*interval1),
                 lwd = 1, position = position_dodge(width = 1/2)) +
  coord_flip()

a

1
  • 回答ありがとうございました!無事に解決致しました。質問の末尾に記しておりました通り、頂いたご回答を英語版のstackoverflowに投稿しております。もし不都合がございましたら、お知らせください。訂正・削除致します。
    – user51966
    Commented 2017年2月3日 9:10
0

map()で要素名、flatten_*()で型をそれぞれ指定し、tibble::tibbleで再びデータフレームとして格納します。
あとは可視化対象の係数名をdplyr::filterなどを使って抽出します。

可視化はgistのリンクと同様ですので省略します。

allModelFrame <- tibble::tibble(
     Variable = lmRes_set$model1 %>% purrr::map("term") %>% purrr::flatten_chr(),
     Coefficient = lmRes_set$model1 %>% purrr::map("estimate") %>% purrr::flatten_dbl(),
     SE  = lmRes_set$model1 %>% purrr::map("std.error") %>% purrr::flatten_dbl(),
     modelName = rep(lmRes_set$setid, each = 2) %>% as.character()
 ) %>%  filter(Variable == "age")

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。