다중 회귀를 통해 다른 변수를 제어하면서 변수 수준으로 분할 된 막대 그래프를 그리는 방법은 무엇입니까?
회귀를 통해 다른 변수를 제어하면서 평균에 대한 막대 그래프를 어떻게 분할 막대 단위 방식으로 그릴 수 있습니까?
내 일반적인 문제
망고, 바나나, 사과 중 어떤 과일이 더 좋아하는지 알아보기 위해 조사를합니다. 이를 위해 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그래픽 에 대한 질문 입니다. 그러나 알 수 있듯이 내 방법은 길고 (즉, 간결하지 않음) 반복적입니다. 특히 처음에 설명한 것처럼 수정되지 않은 수단에 대한 막대 그래프를 얻는 단순성과 관련이 있습니다. 누군가가 코드를 더 짧고 간단하게 만드는 방법에 대한 아이디어를 가지고 있다면 매우 기쁠 것입니다.
답변
데이터 하위 집합에 모델을 맞출 때 원하는 통계적 수량을 얻고 있다고 확신하지 않습니다. 질문하고 싶은 질문을하는 더 좋은 방법은보다 완전한 모델 (모델에 맹인 포함)을 사용한 다음 각 그룹 간의 평균 점수 차이에 대한 모델 대비 를 계산하는 것 입니다.
즉, 여기에 원하는 작업을 수행하는 코드가 있습니다.
- 먼저
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)대신 name및 is_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에 가까워지는 것을 알 수 있습니다. 이것은 부분 군 모델보다 약간 더 잘하고 있다는 신호로 볼 수 있습니다. , 샘플링 된 분포와 관련하여 지상 진실을 알고 있기 때문입니다.