2

I want to make this code into a function with a for loop so that I can plug a variable into the function (where "age" is in this code) and get an ANOVA table out that compares the fit of each polynomial. I would like the fits examined to go up to fit_10.

fit_1 = lm(wage~age, data = Wage)
fit_2 = lm(wage~poly(age,2), data = Wage)
fit_3 = lm(wage~poly(age,3), data = Wage)
fit_4 = lm(wage~poly(age,4), data = Wage)
fit_5 = lm(wage~poly(age,5), data = Wage)
print(anova(fit_1,fit_2,fit_3,fit_4,fit_5))

Any help appreciated on this basic question!

Rachel
  • 21
  • 2

1 Answers1

2

You may use update.formula in an lapply. To update the formula, we loop over 1:5 for desired polynomials and paste/sprintf a formula together in each iteration.

res <- lapply(1:5, \(x) lm(update(mpg ~ 1, paste('. ~', sprintf('poly(hp, %s)', x))), mtcars))

(Note, that if we just do lapply(1:5, \(x) lm(mpg ~ poly(hp, x), mtcars)), only the poly(hp, x) is stored in the call and not e.g. poly(hp, 2).)

Then we do.call anova with the res-list as argument.

do.call('anova', res)
# Analysis of Variance Table
# 
# Model 1: mpg ~ poly(hp, 1)
# Model 2: mpg ~ poly(hp, 2)
# Model 3: mpg ~ poly(hp, 3)
# Model 4: mpg ~ poly(hp, 4)
# Model 5: mpg ~ poly(hp, 5)
#   Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
# 1     30 447.67                                   
# 2     29 274.63  1   173.043 16.8050 0.0003606 ***
# 3     28 269.61  1     5.026  0.4881 0.4909752    
# 4     27 269.52  1     0.087  0.0085 0.9272928    
# 5     26 267.72  1     1.794  0.1742 0.6798122    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
jay.sf
  • 60,139
  • 8
  • 53
  • 110