Building linear regression model from map_*() in R

user_14512983 0
Viewed 155

I am using the attrition data for building a linear model for the relationship between Age & MonthlyIncome for each department in the attrition data. Further, I want to use a map function to add a column r_squared by extracting the real-valued r.squared element from lm_summary which quantifies the goodness of the fit.

My code:

attrition %>%
  group_by(Department) %>%
  mutate(lm_summary = list(summary(lm(MonthlyIncome ~ Age)))) %>%
  mutate(r_squared = map_dfc(lm_summary, "r.squared"))

I am getting the following error upon execution:

essentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableNew names:
* NA -> ...1
New names:
* NA -> ...1
New names:
* NA -> ...1

And this goes into an infinite loop. Not sure if I'm wrongly using the functions or should I interchange them.

> dput(head(attrition))
structure(list(Age = c(41L, 49L, 37L, 33L, 27L, 32L), Attrition = structure(c(2L, 
1L, 2L, 1L, 1L, 1L), .Label = c("No", "Yes"), class = "factor"), 
    BusinessTravel = structure(c(3L, 2L, 3L, 2L, 3L, 2L), .Label = c("Non-Travel", 
    "Travel_Frequently", "Travel_Rarely"), class = "factor"), 
    DailyRate = c(1102L, 279L, 1373L, 1392L, 591L, 1005L), Department = structure(c(3L, 
    2L, 2L, 2L, 2L, 2L), .Label = c("Human_Resources", "Research_Development", 
    "Sales"), class = "factor"), DistanceFromHome = c(1L, 8L, 
    2L, 3L, 2L, 2L), Education = structure(c(2L, 1L, 2L, 4L, 
    1L, 2L), .Label = c("Below_College", "College", "Bachelor", 
    "Master", "Doctor"), class = c("ordered", "factor")), EducationField = structure(c(2L, 
    2L, 5L, 2L, 4L, 2L), .Label = c("Human_Resources", "Life_Sciences", 
    "Marketing", "Medical", "Other", "Technical_Degree"), class = "factor"), 
    EnvironmentSatisfaction = structure(c(2L, 3L, 4L, 4L, 1L, 
    4L), .Label = c("Low", "Medium", "High", "Very_High"), class = c("ordered", 
    "factor")), Gender = structure(c(1L, 2L, 2L, 1L, 2L, 2L), .Label = c("Female", 
    "Male"), class = "factor"), HourlyRate = c(94L, 61L, 92L, 
    56L, 40L, 79L), JobInvolvement = structure(c(3L, 2L, 2L, 
    3L, 3L, 3L), .Label = c("Low", "Medium", "High", "Very_High"
    ), class = c("ordered", "factor")), JobLevel = c(2L, 2L, 
    1L, 1L, 1L, 1L), JobRole = structure(c(8L, 7L, 3L, 7L, 3L, 
    3L), .Label = c("Healthcare_Representative", "Human_Resources", 
    "Laboratory_Technician", "Manager", "Manufacturing_Director", 
    "Research_Director", "Research_Scientist", "Sales_Executive", 
    "Sales_Representative"), class = "factor"), JobSatisfaction = structure(c(4L, 
    2L, 3L, 3L, 2L, 4L), .Label = c("Low", "Medium", "High", 
    "Very_High"), class = c("ordered", "factor")), MaritalStatus = structure(c(3L, 
    2L, 3L, 2L, 2L, 3L), .Label = c("Divorced", "Married", "Single"
    ), class = "factor"), MonthlyIncome = c(5993L, 5130L, 2090L, 
    2909L, 3468L, 3068L), MonthlyRate = c(19479L, 24907L, 2396L, 
    23159L, 16632L, 11864L), NumCompaniesWorked = c(8L, 1L, 6L, 
    1L, 9L, 0L), OverTime = structure(c(2L, 1L, 2L, 2L, 1L, 1L
    ), .Label = c("No", "Yes"), class = "factor"), PercentSalaryHike = c(11L, 
    23L, 15L, 11L, 12L, 13L), PerformanceRating = structure(c(3L, 
    4L, 3L, 3L, 3L, 3L), .Label = c("Low", "Good", "Excellent", 
    "Outstanding"), class = c("ordered", "factor")), RelationshipSatisfaction = structure(c(1L, 
    4L, 2L, 3L, 4L, 3L), .Label = c("Low", "Medium", "High", 
    "Very_High"), class = c("ordered", "factor")), StockOptionLevel = c(0L, 
    1L, 0L, 0L, 1L, 0L), TotalWorkingYears = c(8L, 10L, 7L, 8L, 
    6L, 8L), TrainingTimesLastYear = c(0L, 3L, 3L, 3L, 3L, 2L
    ), WorkLifeBalance = structure(c(1L, 3L, 3L, 3L, 3L, 2L), .Label = c("Bad", 
    "Good", "Better", "Best"), class = c("ordered", "factor")), 
    YearsAtCompany = c(6L, 10L, 0L, 8L, 2L, 7L), YearsInCurrentRole = c(4L, 
    7L, 0L, 7L, 2L, 7L), YearsSinceLastPromotion = c(0L, 1L, 
    0L, 3L, 2L, 3L), YearsWithCurrManager = c(5L, 7L, 0L, 0L, 
    2L, 6L)), row.names = c(NA, -6L), class = c("tbl_df", "tbl", 
"data.frame"))
purrr r
4 Answers

You can use pluck to extract the "r.squared" from each lm model.

library(dplyr)
library(purrr)

attrition %>%
  group_by(Department) %>%
  summarise(lm_summary = list(summary(lm(MonthlyIncome ~ Age))), 
            r_squared = map_dbl(lm_summary, pluck, "r.squared"))

#  Department           lm_summary r_squared
#  <fct>                <list>         <dbl>
#1 Research_Development <smmry.lm>     0.389
#2 Sales                <smmry.lm>     0

You can use mutate instead of summarise if you want to maintain all the rows in the data.

user_39629140

If you are just interested in the R2 I don't think you need the map function: You could just group by department and then extract the R2 directly:

attrition %>%
    group_by(Department) %>%
    mutate(r_squared = summary(lm(MonthlyIncome ~ Age))[['r.squared']])

If you insist on using a map function, you must make sure that you really supply a function:

attrition %>%
    group_by(Department) %>%
    mutate(lm_summary = list(summary(lm(MonthlyIncome ~ Age)))) %>%
    mutate(r_squared = purrr::map_dbl(lm_summary, function(x) x[["r.squared"]]))
user_100308450

Using base R

do.call(rbind, lapply(split(attrition, attrition$Department, drop = TRUE), 
  function(x) {
            smry <- summary(lm(MonthlyIncome ~ Age, data = x))
            r_squared <- smry$r.squared
             data.frame(lm_summary = I(list(smry)), r_squared)

        }))

-output

#                     lm_summary r_squared
#Research_Development lm(formu.... 0.3890383
#Sales                lm(formu.... 0.0000000

Also, we could just use dplyr to do this. Advantage is that we can still get the model summary as a list along with the r.squared without using another package

library(dplyr)
attrition %>%
   nest_by(Department) %>% 
   mutate(lm_summary = list(summary(lm(MonthlyIncome ~ Age, data = data))),
         r_squared = lm_summary$r.squared) %>%
   ungroup
# A tibble: 2 x 4
#  Department                     data lm_summary r_squared
#  <fct>                <list<tibble>> <list>         <dbl>
#1 Research_Development       [5 × 30] <smmry.lm>     0.389
#2 Sales                      [1 × 30] <smmry.lm>     0
user_37322710

You can also use the following solution. Although very similar to some aspects to already posted solutions, I thought it would be nice to post it.

library(dplyr)
library(purrr)
library(tidyr)
library(broom)

attrition %>%
  group_nest(Department) %>% 
  mutate(model = map(data, ~ lm(MonthlyIncome ~ Age, data = .)), 
         summary = map(model, ~ glance(.x))) %>%
  unnest(summary) %>%
  select(Department, r.squared)

# A tibble: 2 x 2
  Department           r.squared
  <fct>                    <dbl>
1 Research_Development     0.389
2 Sales                    0
user_143145200
Related