Solved – R Using linear and quadratic term in regression model

quadratic formrregression

Say I have a model:

log(Y) ~ X1 + I(X1^2)

1) Do I still need to take the log of Y if I am using the quadratic of X. I thought we need to take log to linearise the relationship between Y and X and therefore we don't need to include the squared term of X1. In a way, which of the following three models is the right model:

log(Y) ~ X1 + I(X1^2) 
Y ~ X1 + I(X1^2)  
log(Y) ~ X1

2) Let's say my model is

Y ~ X1 + I(X1^2)

Both the X1 and I(X1^2) comes out to be significant. How do I interpret my regression slopes? How do I make a statement saying.." one unit change in X changes Y by the value given by the regression slope?" How is the quadratic term interpreted?

Thanks

Best Answer

For your first question, you should take log of Y if log of Y makes sense from a substantive point of view. One example where a log often makes sense is money because we often think of money in terms of proportions rather than absolute differences. The difference between a \$100,000 house and a \$150,000 house is large; the difference between a \$10,000,000 house and a \$10,050,000 house is rounding error. Decide whether to transform variables before you start modeling.

For your second question, you cannot make a statement of the type you give, because a 1 unit change in X will have a different effect on Y at different levels of X.

Related Solutions

Solved – Quadratic terms in multiple linear regression

Partially because your other variables 'take credit' of X4sq or X8sq when they are not present in your models yet. Considering the situation when number of shark attack is significantly increased by the increase in ice-cream sales at the beach. But when we add temperature variables to the equation, which in turn causing more people to go to the beach, the ice-cream sale effect is no more.

It would be interesting to check the correlation among variables by using pairs(data)

Further improvement: Only some variables are significantly predictors to age. So you can consider to remove some of the variables and check the R squares using Bayesian Information Criterion.

R – Interpreting Interactions in Linear Models vs Quadratic Models Using lme4-nlme

My favorite way to understand an interaction between two continuous predictors (e.g. heat and year) is to plot it with the full range of one predictor on the x-axis and a few different lines representing potentially interesting values of the other predictor --- I usually pick three, representing low, medium, and high levels of the predictor, but you can play around with it depending on your data.

For your specific case, here is some example code. First, since you didn't provide example data, I'll generate some that (very roughly) matches the parameter estimates you provide:

set.seed(24601) # to get exactly the same results as mine, use this set.seed()

# generate data
year <- 1996:2006
heat <- rnorm(n=20, mean=50, sd=15)
location <- letters[1:10]

library(tidyverse)
my_data <- base::expand.grid(year=year, heat=heat, location=location) %>% 
  mutate(yield = -75411.13 + 38.89*year - 82.56*heat + 0.04*year*heat + rnorm(nrow(.), sd = 50))

Basically, I just used the fixed effect estimates from your model and added a little normal noise. For your own purposes, of course, you will use your real data and real fitted models.

Model 1 (interaction, linear effects only)

# run model
library(lme4)
mdl <- lmer(yield ~ year * heat  + (1|location), data = my_data)

From your question, it seems like you're interested primarily in the overall fixed effects (not the estimates for individual locations, etc.), so I'll focus on that. If you want to get into visualizing the results for individual locations, I recommend this awesome blog post: https://tjmahr.github.io/plotting-partial-pooling-in-mixed-effects-models/

Here are the fixed effects as I estimated them (they're a bit different from your estimates, so your plot will look a bit different, too, but I think it's close enough to be useful):

> (fe <- summary(mdl)$coefficients[1:4,1])
  (Intercept)          year          heat     year:heat 
-74219.711477     38.293178   -117.339456      0.057416

Now let's plot it to take a look at that interaction. I'll choose to put heat on the x-axis and have three representative lines for low, medium, and high values of year, but you could switch and do it the other way if you prefer. I'll use those values (the full range of heat and the three example values of year) to generate predicted yield scores for each combination based on the fixed effects estimates from the model.

# select values for plotting
# full range of heat, and selected low, med and high values for year
plot_df <- expand.grid(heat=heat, year=c(min(year), mean(year), max(year))) %>% 
  mutate(yield = fe[1] + fe[2]*year + fe[3]*heat + fe[4]*year*heat)

Now I'll plot those predicted yield values as lines.

plot_interaction <- ggplot(plot_df, aes(y=yield, x = heat, color = as.factor(year))) + 
  geom_line(size=2) + 
  labs(color = "year") + 
  # tweak plot appearance
  scale_color_brewer(palette = "Dark2") + 
  theme(legend.position = "top") + 
  theme_classic()

plot_interaction

So what do we see here? First of all, the interaction is very small relative to the individual effects of year and heat --- the lines look almost parallel. But you know from the parameter estimates that it is significantly different from zero, so it's there even if it's small. heat has a negative effect when year = 0 (ideally, you will have centered both predictors before estimating the model, so this would be "at the mean of year"), and the positive interaction term indicates that as year increases, the effect of heat gets less negative, so it weakens. What you should look for in the plot is the slope to be a bit shallower for higher years than for low years.

Model 2 (linear and quadratic effects, with interactions)

The second model seems much more complicated, but it's actually just as easy to plot, almost! We'll do the same procedure as before, first getting the fixed effects from the model then generating some plotting data with the full range of heat, selected values of year, and predicted values for yield based on the fixed effect estimates.

First, since I don't have your real data, I'm re-generating another dataset that will (roughly) match the parameter estimates of the second model, so we can get a more or less accurate plot.

my_data <- base::expand.grid(year=year, heat=heat, location=location) %>% 
  mutate(yield = -97816.61 + 50.07*year + 20499.87*heat - 632.2*heat*heat - 10.22*year*heat +.31*year*heat*heat + rnorm(nrow(.), sd = 50))

Run the second model:

mdl <- lmer(yield ~ year*(heat + I(heat^2)) + (1|location), data = my_data)

Get the fixed effects:

fe <- summary(mdl)$coefficients[1:6,1]

> round(fe, 3)
   (Intercept)           year           heat      I(heat^2)      year:heat year:I(heat^2) 
    -95838.400         49.083      20410.924       -631.207        -10.176          0.310

Close enough. :) Now generate the plotting data.

# full range of heat, and selected low, med and high values for year
plot_df <- expand.grid(heat=heat, year=c(min(year), mean(year), max(year))) %>% 
  mutate(yield = fe[1] + fe[2]*year + fe[3]*heat + fe[4]*heat*heat + fe[5]*year*heat + fe[6]*year*heat*heat)

We can use exactly the same plotting code as before, just feed it the updated dataframe:

plot_interaction %+% plot_df

Now we can see that the general effect of heat on yield is negative, with an accelerating effect, so that there's a steeper negative slope at higher levels of heat. Moreover, we can see that the effect of heat depends on year (i.e. the interaction), such that heat drops off more sharply in earlier years than later years.

You should be able to generate plots for your own data using this code. I always find having a visual is extremely helpful when interpreting interactions, especially in models with more than a few coefficients (such as your model 2).

(Note that if the ranges for heat and year are pretty different in your data to what I came up with generating these data, then your plots might look quite different even though the parameter estimates are similar.)

A final note on plotting and interpreting interactions

My guess from your parameter estimates is that you didn't center year and heat before estimating the model, which is why I didn't either. But you probably should. Not only does it ease some model estimation issues like multicollinearity, it makes interpreting your interactions easier since 0 values are more meaningful. For example, the parameter estimates you're getting for heat are the predicted change in yield for each additional unit of heat when year = 0. Unless you've centered year (or you're studying heat and yields in antiquity), you're probably talking about values far outside of the range of your data. Your estimate for the effect of heat will be much more meaningful if you center first.