# Why does lm() estimate different coefficients for polynomial or quadratic terms?

I have a survival dataset H.month where I look at hazard per month.

> monthly
# A tibble: 144 x 11
# Groups:   company 
time n.risk n.event n.censor      surv     std.err     upper     lower        company  hazardrate
<dbl>  <dbl>   <dbl>    <dbl>     <dbl>       <dbl>     <dbl>     <dbl>          <chr>       <dbl>
1     0   2538       0        0 1.0000000          NA        NA        NA Consumer Goods          NA
2     1   2538       9        9 0.9964539 0.001179935 0.9987692 0.9941440 Consumer Goods 0.003546099
3     2   2520       5        7 0.9944768 0.001472058 0.9973662 0.9915958 Consumer Goods 0.001984127
4     3   2508       9       14 0.9909081 0.001887172 0.9946138 0.9872162 Consumer Goods 0.003588517
5     4   2485       9        4 0.9873193 0.002227452 0.9916947 0.9829632 Consumer Goods 0.003621730
6     5   2472       8        6 0.9841241 0.002490285 0.9890171 0.9792553 Consumer Goods 0.003236246
7     6   2458       8       16 0.9809211 0.002727536 0.9862816 0.9755898 Consumer Goods 0.003254679
8     7   2434       7      109 0.9781000 0.002920678 0.9838412 0.9723923 Consumer Goods 0.002875924
9     8   2318      14       37 0.9721926 0.003302311 0.9786866 0.9657417 Consumer Goods 0.006039689
10     9   2267       5      157 0.9700484 0.003431432 0.9767972 0.9633462 Consumer Goods 0.002205558
# ... with 134 more rows, and 1 more variables: time2 <dbl>


I want to examine whether there is a curvilinear relationship between the time (in months) and the hazard.

When I estimate the model with a quadratic or a polynomial term, I get similar model statistics but different coefficients. The also happens when I estimate an unweighted model

> monthly$time2 <- monthly$time*monthly\$time
> lm(hazardrate ~ time + time2, monthly) %>% summary()

Call:
lm(formula = hazardrate ~ time + time2, data = monthly)

Residuals:
Min         1Q     Median         3Q        Max
-0.0061292 -0.0036970 -0.0009177  0.0025068  0.0263357

Coefficients:
Estimate    Std. Error t value  Pr(>|t|)
(Intercept)  0.0006151737  0.0012128566   0.507     0.613
time         0.0003116329  0.0000685731   4.545 0.0000119 ***
time2       -0.0000033660  0.0000008221  -4.095 0.0000714 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.005017 on 139 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared:  0.135, Adjusted R-squared:  0.1226
F-statistic: 10.85 on 2 and 139 DF,  p-value: 0.00004187

> lm(hazardrate ~ poly(time, 2), monthly) %>% summary()

Call:
lm(formula = hazardrate ~ poly(time, 2), data = monthly)

Residuals:
Min         1Q     Median         3Q        Max
-0.0061292 -0.0036970 -0.0009177  0.0025068  0.0263357

Coefficients:
Estimate Std. Error t value             Pr(>|t|)
(Intercept)     0.0057716  0.0004213  13.699 < 0.0000000000000002 ***
poly(time, 2)1  0.0123404  0.0051170   2.412               0.0172 *
poly(time, 2)2 -0.0211600  0.0051678  -4.095            0.0000714 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.005017 on 139 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared:  0.135, Adjusted R-squared:  0.1226
F-statistic: 10.85 on 2 and 139 DF,  p-value: 0.00004187


Can someone explain why this is happening, how to interpret, and which to use?

• You want to explain hazard by month2 and on the other hand by month and month2. Without going into detail (since this question should be posted over at crossvalidated), isnt it obvious that the explanatory value of month2 changes if you use another explanatory variable as well.... – Martin Schmelzer Aug 4 '17 at 12:23
• google for 'omitted variable bias' , just an example: en.wikipedia.org/wiki/Omitted-variable_bias – user83346 Aug 4 '17 at 13:45

If we append another explanatory variable then it changes the coefficients of the original explanatory variables unless the new variable is orthogonal (think independent) to the original explanatory variables. If the new variable is not orthogonal then both the original variables and then new variable redundantly explain a certain portion of the dependent variable and so it effectively allocates a bit to each. If the new variable is orthongonal then it explains entirely different part of the dependent variable so the coefficients will not change when the new one is added. Thus in most cases the coefficient of an explanatory variable is only meaningful within the context of what other explanatory variables were used.

Regarding which to use, i.e. whehter to use an explanatory variable or not, you can use the anova command to perform a statistical test.

For example, we regress mpg on an intercept and cyl in (1) and then on the intercept, cyl and also disp in (2) and we note the coefficients have changed.

In the example (3) we take the projection of disp orthogonal to the plane of the intercept and cyl and and then double check that the resulting disp_orth is orthongal to cyl and orthogonal to the intercept. The projection of disp orthognonal to the intercept and cyl is that part of disp that is not already explained by the intercept and cyl. Two vectors v1 and v2 are orthogonal if crossprod(v1, v2) = 0 and in the case of the intercept crossprod(intercept, v) = sum(v) = 0. Now, using this new orthogonal column when we run the regression of mpg against cyl and disp_orth we see that the coefficients of the intercept and cyl are the same as in (1). That is the intercept is 37.88457649 in both (1) and (3) and the coefficient of cyl is -2.87579014 in both (1) and (3).

Note that (2) and (3) are actually the same regressions except for paramterization. For example, they both give the same fitted values as predictions.

# 1
coef(fm1 <- lm(mpg ~ cyl, mtcars))
## (Intercept)         cyl
##    37.88458    -2.87579

# 2
coef(fm2 <- lm(mpg ~ cyl + disp, mtcars))
## (Intercept)         cyl        disp
## 34.66099474 -1.58727681 -0.02058363

# 3
disp_orth <- resid(lm(disp ~ cyl, mtcars))

# check orthogonality of disp_orth and cyl
with(mtcars, crossprod(cyl, disp_orth))
##              [,1]
## [1,] 2.557954e-13

# check orthogonality of disp_orth and intercept
sum(disp_orth)
##  2.007283e-13

coef(fm3 <- lm(mpg ~ cyl + disp_orth, mtcars))
## (Intercept)         cyl   disp_orth
## 37.88457649 -2.87579014 -0.02058363

all.equal(fitted(fm2), fitted(fm3))
##  TRUE


Also note that the anova values are the same for comparing fm1 vs. fm2 and fm1 vs. fm3. In both case the p-value is slighly greater than 5% so on that bassis we would use fm1 as our model rather than fm2 or fm3.

> anova(fm1, fm2)
Analysis of Variance Table

Model 1: mpg ~ cyl
Model 2: mpg ~ cyl + disp
Res.Df    RSS Df Sum of Sq      F  Pr(>F)
1     30 308.33
2     29 270.74  1    37.594 4.0268 0.05419 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(fm1, fm3)
Analysis of Variance Table

Model 1: mpg ~ cyl
Model 2: mpg ~ cyl + disp_orth
Res.Df    RSS Df Sum of Sq      F  Pr(>F)
1     30 308.33
2     29 270.74  1    37.594 4.0268 0.05419 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1