다중 회귀를 통해 다른 변수를 제어하면서 변수 수준으로 분할 된 막대 그래프를 그리는 방법은 무엇입니까?

Aug 20 2020

회귀를 통해 다른 변수를 제어하면서 평균에 대한 막대 그래프를 어떻게 분할 막대 단위 방식으로 그릴 수 있습니까?

내 일반적인 문제

망고, 바나나, 사과 중 어떤 과일이 더 좋아하는지 알아보기 위해 조사를합니다. 이를 위해 100 명을 무작위로 샘플링합니다. 나는 그들에게 각각의 과일을 좋아하는 정도를 1-5의 척도로 평가하도록 요청한다. 또한 성별, 연령, 교육 수준, 색맹 여부와 같은 인구 통계 학적 정보도 수집합니다. 색각이 결과를 바꿀 수 있다고 생각하기 때문입니다. 그러나 내 문제는 데이터 수집 후 샘플이 일반 인구를 잘 나타내지 않을 수 있다는 것을 깨닫습니다. 나는 80 %의 남성이있는 반면, 인구의 성별은 더 고르게 나뉩니다. 내 샘플의 교육 수준은 매우 균일하지만 인구에서는 박사 학위를 소지하는 것보다 고등학교 졸업장 만 보유하는 것이 더 일반적입니다. 나이도 대표적이지 않습니다.

따라서 내 표본을 기반으로 과일 선호도를 계산하는 수단은 인구 수준에 대한 결론을 일반화하는 측면에서 제한 될 가능성이 있습니다. 이 문제를 처리하는 한 가지 방법은 편향된 인구 통계 데이터를 제어하기 위해 다중 회귀를 실행하는 것입니다.

색각 수준 (색맹 여부)에 따라 막대를 나란히 분할하는 막대 그래프에서 회귀 결과를 플로팅하고 싶습니다.

내 데이터

library(tidyverse)

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), 100, replace = TRUE),
    i_love_mango = sample(c(1:5), 100, replace = TRUE),
    age = sample(c(20:70), 100, replace = TRUE),
    is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
    education_level = sample(c(1:4), 100, replace = TRUE),
    is_colorblinded = sample(c(0, 1), 100, 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


각 과일 좋아하는 수준에 대한 평균 값을 얻으려면

fruit_liking_df_for_barplot <-
  fruit_liking_df %>%
  pivot_longer(.,
    cols = c(i_love_apple, i_love_banana, i_love_mango),
    names_to = "fruit",
    values_to = "rating") %>%
  select(id, fruit, rating, everything())

ggplot(fruit_liking_df_for_barplot, aes(fruit, rating, fill = as_factor(is_colorblinded))) +
  stat_summary(fun = mean,
               geom = "bar",
               position = "dodge") +
  ## errorbars
  stat_summary(fun.data = mean_se,
               geom = "errorbar",
               position = "dodge") +
  ## bar labels
  stat_summary(
    aes(label = round(..y.., 2)),
    fun = mean,
    geom = "text",
    position = position_dodge(width = 1),
    vjust = 2,
    color = "white") +
  scale_fill_discrete(name = "is colorblind?",
                      labels = c("not colorblind", "colorblind")) +
  ggtitle("liking fruits, without correcting for demographics")

그러나 인구를 더 잘 나타 내기 위해 이러한 의미를 수정하려면 어떻게해야합니까?

다중 회귀를 사용할 수 있습니다.

  • 45 세인 인구의 평균 연령을 수정하겠습니다.

  • 섹스에 대한 올바른 50-50 분할을 수정하겠습니다.

  • 고등학교 ( 2내 데이터에 코딩 됨) 인 일반 교육 수준으로 수정합니다.

  • 나이가 비선형적인 방식으로 과일의 취향에 영향을 미친다고 믿을만한 이유가 있기 때문에 이에 대해서도 설명하겠습니다.

lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)

동일한 모델을 통해 세 가지 과일 데이터 (사과, 바나나, 망고)를 실행하고 절편을 추출하여 인구 통계 데이터를 제어 한 후이를 수정 된 평균으로 간주합니다.

먼저 색맹 인 데이터에 대해서만 회귀 분석을 실행하겠습니다.

library(broom)

dep_vars <- c("i_love_apple",
              "i_love_banana",
              "i_love_mango")

regresults_only_colorblind <-
  lapply(dep_vars, function(dv) {
    tmplm <-
      lm(
        get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2), 
        data = filter(fruit_liking_df, is_colorblinded == 1)
      )
    
    broom::tidy(tmplm) %>%
      slice(1) %>%
      select(estimate, std.error)
  })

data_for_corrected_barplot_only_colorblind <-
  regresults_only_colorblind %>%
  bind_rows %>%
  rename(intercept = estimate) %>%
  add_column(dep_vars, .before = c("intercept", "std.error")) 

## # A tibble: 3 x 3
##   dep_vars      intercept std.error
##   <chr>             <dbl>     <dbl>
## 1 i_love_apple       3.07     0.411
## 2 i_love_banana      2.97     0.533
## 3 i_love_mango       3.30     0.423

그런 다음 색맹에 대해서만 수정 된 막대 그래프를 플로팅합니다.

ggplot(data_for_corrected_barplot_only_colorblind, 
       aes(x = dep_vars, y = intercept)) +
  geom_bar(stat = "identity", width = 0.7, fill = "firebrick3") +
  geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
                width = 0.2) +
  geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
  ggtitle("liking fruits after correction for demogrpahics \n colorblind subset only")

둘째, 색각 만있는 데이터에 대해 동일한 회귀 프로세스를 반복합니다.

dep_vars <- c("i_love_apple",
              "i_love_banana",
              "i_love_mango")

regresults_only_colorvision <-
  lapply(dep_vars, function(dv) {
    tmplm <-
      lm(
        get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2), 
        data = filter(fruit_liking_df, is_colorblinded == 0) ## <- this is the important change here
      )
    
    broom::tidy(tmplm) %>%
      slice(1) %>%
      select(estimate, std.error)
  })


data_for_corrected_barplot_only_colorvision <-
  regresults_only_colorvision %>%
  bind_rows %>%
  rename(intercept = estimate) %>%
  add_column(dep_vars, .before = c("intercept", "std.error")) 

ggplot(data_for_corrected_barplot_only_colorvision, 
       aes(x = dep_vars, y = intercept)) +
  geom_bar(stat = "identity", width = 0.7, fill = "orchid3") +
  geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
                width = 0.2) +
  geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
  ggtitle("liking fruits after correction for demogrpahics \n colorvision subset only")



내가 궁극적으로 찾고있는 것은 수정 된 플롯을 결합하는 것입니다.


마지막 메모

이것은 주로 ggplot그래픽 에 대한 질문 입니다. 그러나 알 수 있듯이 내 방법은 길고 (즉, 간결하지 않음) 반복적입니다. 특히 처음에 설명한 것처럼 수정되지 않은 수단에 대한 막대 그래프를 얻는 단순성과 관련이 있습니다. 누군가가 코드를 더 짧고 간단하게 만드는 방법에 대한 아이디어를 가지고 있다면 매우 기쁠 것입니다.

답변

1 BrianLang Aug 20 2020 at 16:37

데이터 하위 집합에 모델을 맞출 때 원하는 통계적 수량을 얻고 있다고 확신하지 않습니다. 질문하고 싶은 질문을하는 더 좋은 방법은보다 완전한 모델 (모델에 맹인 포함)을 사용한 다음 각 그룹 간의 평균 점수 차이에 대한 모델 대비 를 계산하는 것 입니다.

즉, 여기에 원하는 작업을 수행하는 코드가 있습니다.

  • 먼저 pivot_longer데이터가 긴 형식이되도록 과일 열입니다.
  • 그런 다음 group_by과일 유형과 실명 변수 및 nest각 과일 유형 및 실명 범주에 대한 별도의 데이터 세트를 제공하는 호출 입니다.
  • 그런 다음 purrr::map각 데이터 세트에 모델을 맞추는 데 사용 합니다.
  • broom::tidy그리고 broom::confint_tidy우리에게 우리가 모델에 대해 원하는 통계를 제공합니다.
  • 그런 다음 모델 요약의 중첩을 해제하고 절편에 해당하는 행으로 구체적으로 필터링해야합니다.
  • 이제 그림을 만드는 데 필요한 데이터가 있습니다. 나머지는 여러분에게 맡기겠습니다.
library(tidyverse)

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), 100, replace = TRUE),
    i_love_mango = sample(c(1:5), 100, replace = TRUE),
    age = sample(c(20:70), 100, replace = TRUE),
    is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
    education_level = sample(c(1:4), 100, replace = TRUE),
    is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
  )

model_fits <- fruit_liking_df %>%
  pivot_longer(starts_with("i_love"), values_to = "fruit") %>% 
  group_by(name, is_colorblinded) %>%
  nest() %>% 
  mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
                                      I((age - 45)^2) +
                                      I(is_male - 0.5) + 
                                      I(education_level - 2))),
         model_summary = map(model_fit, ~ bind_cols(broom::tidy(.x), broom::confint_tidy(.x)))) 

model_fits %>%
  unnest(model_summary) %>%
  filter(term == "(Intercept)") %>% 
  ggplot(aes(x = name, y = estimate, group = is_colorblinded,
             fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
  geom_bar(stat = "identity", position = position_dodge(width = .95)) +
  geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
                colour = "black", width = .15, position = position_dodge(width = .95))

편집하다


차라리 단일 모델을 적용하려는 경우 (따라서 표본 크기를 늘리고 추정값을 줄임). is_colorblind를 모델에 factor.

lm(data = .x, fruit ~ I(age - 45) +
 I((age - 45)^2) + I(is_male - 0.5) + 
 I(education_level - 2) + 
 as.factor(is_colorblind))

그런 다음 "색맹이 아닌 평균적인 사람"과 "색맹이 아닌 평균적인 사람"이라는 두 가지 관찰에 대한 예측을 얻고 싶을 것입니다.

new_data <- expand_grid(age = 45, is_male = .5, 
                        education_level = 2.5, is_colorblinded = c(0,1))

그런 다음 이전과 같이 새 모델을 일부 함수형 프로그래밍으로 피팅 할 수 있지만 group_by(name)대신 nameis_colorblind.

model_fits_ungrouped <- fruit_liking_df %>%
  pivot_longer(starts_with("i_love"), values_to = "fruit") %>% 
  group_by(name) %>%
  tidyr::nest() %>% 
  mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
                                      I((age - 45)^2) +
                                      I(is_male - .5) + 
                                      I(education_level - 2) +
                                      as.factor(is_colorblinded))),
         predicted_values = map(model_fit, ~ bind_cols(new_data, 
                                                       as.data.frame(predict(newdata = new_data, .x, 
                                                                             type = "response", se.fit = T))) %>%
                                  rowwise() %>%
                                  mutate(estimate =  fit, 
                                         conf.low =  fit - qt(.975, df) * se.fit, 
                                         conf.high = fit + qt(.975, df) * se.fit)))

이를 통해 이전 플로팅 코드를 약간 변경합니다.

model_fits_ungrouped %>%
  unnest(predicted_values) %>%
  ggplot(aes(x = name, y = estimate, group = is_colorblinded,
             fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
geom_bar(stat = "identity", position = position_dodge(width = .95)) +
 geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
                colour = "black", width = .15, position = position_dodge(width = .95))

그룹화되고 부분 군화 된 두 그림을 비교하면 신뢰 구간이 줄어들고 평균에 대한 추정치가 대부분 3에 가까워지는 것을 알 수 있습니다. 이것은 부분 군 모델보다 약간 더 잘하고 있다는 신호로 볼 수 있습니다. , 샘플링 된 분포와 관련하여 지상 진실을 알고 있기 때문입니다.