# Difficulties obtaining valid predictions when using interactions

I examine long term trends (2003 to 2014) for a continuous dependent variable. I want to predict the mean each year in relation to income category. Income is arranged in quintiles, from 1 (poorest) to 5 (richest). The richest should have the lowest predicted values, there should be a stepwise increase in predicted values with declining income.

When I fit the models separately for each year I can clearly see that pattern. However, it is desirable to fit all the data at once by using interactions. This would allow me to examine both trends and compare the income categories each year. However, when I try using interactions for this purpose, I obtain implausible estimates.

I will illustrate how I fit the analysis stratified by year

Variables:

dependent_var = the outcome variable

id = person identifier (used to account for repeated measurements)

biomarker = a blood biomarker

income = income quintiles; 1 = the poorest, 5 = the richest.

(remaining variables are obvious)

    # As is evident from this analysis, the predicted mean increases  with declining income, and we should see this in a stepwise fashion.
# We start by fitting separate regressions for each year, from 2003 to 2014

library(dplyr)
library(lme4)
library(lsmeans)
library(ggplot2)

y2003 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2003"))
y2004 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2004"))
y2005 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2005"))
y2006 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2006"))
y2007 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2007"))
y2008 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2008"))
y2009 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2009"))
y2010 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2010"))
y2011 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2011"))
y2012 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2012"))
y2013 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2013"))
y2014 <- lmer(dependent_var ~ income +  age + sex + biomarker + education + (1 | id), data=filter(data, year=="2014"))

# Predicting mean value each year by income
r2003 <- summary(lsmeans(y2003, "income")); r2003$category <- "2003" r2004 <- summary(lsmeans(y2004, "income")); r2004$category <- "2004"
r2005 <- summary(lsmeans(y2005, "income")); r2005$category <- "2005" r2006 <- summary(lsmeans(y2006, "income")); r2006$category <- "2006"
r2007 <- summary(lsmeans(y2007, "income")); r2007$category <- "2007" r2008 <- summary(lsmeans(y2008, "income")); r2008$category <- "2008"
r2009 <- summary(lsmeans(y2009, "income")); r2009$category <- "2009" r2010 <- summary(lsmeans(y2010, "income")); r2010$category <- "2010"
r2011 <- summary(lsmeans(y2011, "income")); r2011$category <- "2011" r2012 <- summary(lsmeans(y2012, "income")); r2012$category <- "2012"
r2013 <- summary(lsmeans(y2012, "income")); r2013$category <- "2013" r2014 <- summary(lsmeans(y2014, "income")); r2014$category <- "2014"

r2003$year_cat <- "Y 2003" r2004$year_cat <- "Y 2004"
r2005$year_cat <- "Y 2005" r2006$year_cat <- "Y 2006"
r2007$year_cat <- "Y 2007" r2008$year_cat <- "Y 2008"
r2009$year_cat <- "Y 2009" r2010$year_cat <- "Y 2010"
r2011$year_cat <- "Y 2011" r2012$year_cat <- "Y 2012"
r2013$year_cat <- "Y 2013" r2014$year_cat <- "Y 2014"

# Binding together
data2 <- rbind(r2003, r2004, r2005, r2006, r2007, r2008, r2009, r2010, r2011, r2012, r2013, r2014)

# Plot
ggplot() +
geom_smooth(data=data2, aes(x=factor(year_cat), y=lsmean, colour = income, group=income), se=F, size=2) +
geom_point(data=data2, aes(x=factor(year_cat), y=lsmean, colour = income, group=income), size=3)


I need to fit all the data simultaneously, which can presumably only be done by using interactions. This will allow me to compare both the trends over time and compare income categories annually.

I have tried many possible interactions, both simple terms and more complex.

Can anyone figure out how to fit the model I desire?

• I appreciate you supplying the data and code. However, it is not reproducible as-is: please indicate which libraries you are using for functions such as lmer, lsmeans, filter, and ggplot. Indicate how to read in the data file (which is large--it has 100,000 rows--and uses nonstandard delimiters). Couldn't you possibly strip the input to a much smaller example so that your readers don't have to wait to perform so much calculation just to see what's going on?
– whuber
Commented Jun 29, 2015 at 22:08
• @whuber Thanks. I've made the post more succinct now and replaced the data file with a conventional .RData file. I did not, however, figure out how to reduce the amount of code, but I believe it is reproducible as is, if the file is loaded to the environment. =) Commented Jun 29, 2015 at 22:41
• Isn't this the same question you had here? Commented Jul 2, 2015 at 22:25
• Two things. First, if your sample size varies between years and other random effects then the partial pooling phenomena of random effects model will cause their results to be different than data fit on subsets of the data, but usually this won't be too huge. I forget how lsmeans() works exactly but it may produce different results when fed models with different amounts of data. Commented Apr 20 at 23:01

It seems to me that you need to fit something like

comb <- lmer(dependent_var ~ (income +  age + sex + biomarker +
education)*year + (1 | id), data = data)


... and subsequently do things like

lsmeans(comb, ~ income | year)


which should give you about the same estimates as you had before, but will be based on combining the variance components for all the years.

The reason is that with year interacting with everything, you are basically saying you want to fit the rest of the fixed effects separately for each year. You can do

library(car)
Anova(comb)


and perhaps identify interaction terms that are not really needed.

• Thanks @rvl for the answer. It does improve the plausibility of the results but it is still not fully plausible; we should be seing a clear stepwise gradient toward higher values as income declines. this is evident from crude data, from subject matter knowledge and when stratifyng on calendar year. 😳 Commented Jul 2, 2015 at 15:06
• If I'm thinking straight, you could use (1 + year | id) as the error part and you should get the same fits as separate models using (1 | id). Fit that model and give it a different name, say comb2), and use anova(comb, comb2) to compare them. Commented Jul 2, 2015 at 20:48
• It returns an error implyimg that the number of observations are less than number of random effects. I'll see if I can figure out how to solve it. Commented Jul 3, 2015 at 10:18
• Oops, I think I had it wrong. Try (1|year:id) - though now I have no confidence at all ... Commented Jul 3, 2015 at 15:42