Уровни факторов и моделирование в R

Question

Следующий код работает очень просто lm() и пытается суммировать результаты (уровень фактора, коэффициент) в небольшом кадре данных:

df <- data.frame(star_sign = c("Aries", "Taurus", "Gemini", "Cancer", "Leo", "Virgo", "Libra", "Scorpio", "Sagittarius", "Capricorn", "Aquarius", "Pisces"),
                 y = c(1.1, 1.2, 1.4, 1.3, 1.8, 1.6, 1.4, 1.3, 1.2, 1.1, 1.5, 1.3))

levels(df$star_sign) #alphabetical order

# fit a simple linear model

my_lm <- lm(y ~ 1 + star_sign, data = df)
summary(my_lm) # intercept is based on first level of factor, aquarius

# I want the levels to work properly 1..12 = Aries, Taurus...Pisces so I'm going to redefine the factor levels

df$my_levels <- c("Aries", "Taurus", "Gemini", "Cancer", "Leo", "Virgo", "Libra", "Scorpio", "Sagittarius", "Capricorn", "Aquarius", "Pisces")

df$star_sign <- factor(df$star_sign, levels = df$my_levels)

my_lm <- lm(y ~ 1 + star_sign_, data = df)
summary(my_lm) # intercept is based on first level of factor which is now Aries

# but for my model fit I want the reference level to be Virgo (because reasons)

df$star_sign_2 <- relevel(df$star_sign, ref = "Virgo")

my_lm <- lm(y ~ 1 + star_sign_2, data = df)
summary(my_lm)

df_results <- data.frame(factor_level = names(my_lm$coefficients), coeff = my_lm$coefficients )

# tidy up
rownames(df_results) <- 1:12
df_results$factor_level <- as.factor(gsub("star_sign_2", "", df_results$factor_level))

# change label of "(Intercept)" to "Virgo"
df_results$factor_level <- plyr::revalue(df_results$factor_level, c("(Intercept)" = "Virgo"))

levels(df_results$factor_level) # the levels are alphabetical + Virgo at the front (not same as display order from lm)

Уровни фактора не в правильном порядке : Я хочу отсортировать df_results так, чтобы знаки зодиака появлялись в том же порядке, в котором они изначально (Овен, Телец ... Рыбы), как показано в столбце df$my_levels. Я не думаю, что у меня есть хорошее понимание манипулирования факторами и их ярлыками / уровнями и т. Д. c. поэтому я изо всех сил пытаюсь узнать, как это сделать.

Кроме того, это довольно скучный и неуклюжий фрагмент кода. Есть ли более краткие способы сделать что-то подобное?

Спасибо.

(ps математически модель явно тривиальна, но это нормально для этих целей - меня просто интересует, как манипулировать выходами)

Answer 1

Вот как я бы подошел к извлечению коэффициента модели, используя пакет broom (и dplyr):

library(broom)
library(dplyr)
broom::tidy(my_lm) %>%
  mutate(term = sub("star_sign_2", "", term),
         term = ifelse(term == "(Intercept)", "Virgo", term),
         term = factor(term, levels = unique(term)))
# A tibble: 12 x 5
   term        estimate std.error statistic p.value
   <fct>          <dbl>     <dbl>     <dbl>   <dbl>
 1 Virgo          1.6         NaN       NaN     NaN
 2 Aries         -0.500       NaN       NaN     NaN
 3 Taurus        -0.4         NaN       NaN     NaN
 4 Gemini        -0.2         NaN       NaN     NaN
 5 Cancer        -0.300       NaN       NaN     NaN
 6 Leo            0.20        NaN       NaN     NaN
 7 Libra         -0.2         NaN       NaN     NaN
 8 Scorpio       -0.3         NaN       NaN     NaN
 9 Sagittarius   -0.4         NaN       NaN     NaN
10 Capricorn     -0.500       NaN       NaN     NaN
11 Aquarius      -0.1         NaN       NaN     NaN
12 Pisces        -0.300       NaN       NaN     NaN

Установка levels = unique(term) - это хороший трюк для помещения уровней в порядок, в котором они встречаются.

Другой совет, который я имею, - держать вектор уровней в порядке, который вы хотите , а не в кадре данных, и затем обращаться к нему всякий раз, когда вам нужно установить sh порядок. Например,

astro_order = c("Aries", "Taurus", "Gemini", "Cancer", "Leo", "Virgo", "Libra", "Scorpio", "Sagittarius", "Capricorn", "Aquarius", "Pisces")

# messy but effective:
astro_order_virgo1 = factor(astro_order, levels = astro_order) %>% 
  relevel("Virgo") %>%
  levels()

Таким образом, вы можете заменить последний шаг, описанный выше, на term = factor(term, levels = astro_order_virgo1).

. Этот метод разделения порядка уровней хорош, потому что (a) он не будет изменить, если вы переупорядочиваете свой фрейм данных, и (b) он работает так же хорошо, если ваш фрейм данных длинный и у вас есть повторяющиеся записи ваших уровней факторов.

Answer 2

Если я понял, что вам нужно сделать, это довольно просто. Просто добавьте следующий код в конец вашего скрипта. Я также призываю вас погрузиться в dplyr или tidyverse. Дайте мне знать, если у вас есть какие-либо вопросы:)

## ADDED: 

#WE CREATE AN ID to maintain order in df_results 
df$id <- 1:nrow(df)


library(dplyr)
#Perform left _ join (you could also do inner or right, you'll get the same result in this case )
df_results = left_join(df_results,df, by=c('factor_level'='star_sign_2'))
df_results = df_results %>% arrange(id)

# select desired columns (optionally) 
df_results = df_results %>% select(factor_level,coeff) 


head(df_results)

 factor_level coeff
1        Aries  -0.5
2       Taurus  -0.4
3       Gemini  -0.2
4       Cancer  -0.3
5          Leo   0.2
6        Virgo   1.6