# Visualising a linear model with 6 predictors in R

I've recently been wondering how to visualise a linear model with >1 predictor. Today I came across a paper that managed to visualise a linear model with 6 predictors (and one response variable) in a scatterplot.

The model was similar to this:

weight ~ days since 1st Sept + site + age + year + PC1 + arrival date

Weight was a continuous response variable. Days since 1st Sept was a integer predictor, site was a factor predictor (2 levels), age was a factor predictor (2 levels), year was a factor predictor (2 levels), PC1 was a continuous predictor and arrival date was a integer predictor.

On the y axis of the scatterplot was values of weight predicted by the model. On the x axis of the scatterplot was days since 1st Sept. Eight lines of best fit showed predicted values at each level of the three factor predictors. PC1 and arrival date were 'held constant at their mean values'.

The plot looked like this (made using Photoshop): I'm hoping someone can recreate a plot similar to this using a dataset installed in R? I'm particularly interested in how to keep two predictors 'constant at their mean values', to enable visualisation of the other predictors. The mtcars dataset might be appropriate.

• Are the various lines meant to be non-parallel or is that just the inherent imperfection of photoshop? – gung - Reinstate Monica Aug 21 '13 at 18:28
• To keep things more simple, just assume lines are parallel! – luciano Aug 21 '13 at 18:30
• We can give you something more specific if you can specify the distributions (ranges, etc.) of the variables. But we can also provide a generic solution which you can adapt by changing parts to your liking. – gung - Reinstate Monica Aug 21 '13 at 18:41
• The more generic solution would be fine. I'm just after the R code to hold two predictors constant at their mean whilst visualising the rest of the model. – luciano Aug 21 '13 at 18:48

Here is some code that is hopefully self-explanatory:

set.seed(20987)     # for reproducability

N = 200

# variables
days_since   = rpois(N, lambda=60)
site         = factor(sample(c("site1", "site2"), N, replace=T), c("site1", "site2"))
year         = factor(sample(c("2012", "2013"),   N, replace=T), c("2012", "2013"))
PC1          = rnorm(N, mean=100, sd=25)
arrival_date = sample.int(365, N, replace=T)

# betas
B0  =  13
Bds =  74
Bs  = 114
Ba  = 160
By  = 191
Bpc =  59

# response variable
weight = B0 + Bds*days_since + Bs*(site=="site2") + Ba*(age=="adult") +
By*(year=="2013") + Bpc*PC1 + Bad*arrival_date + rnorm(N, mean=0, sd=10)

model = lm(weight~days_since+site+age+year+PC1+arrival_date)

# predicted values for plot
ds    = seq(min(days_since), max(days_since))
ds1j2 = predict(model, data.frame(days_since=ds, site="site1", age="juv",
year="2012", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds1j3 = predict(model, data.frame(days_since=ds, site="site1", age="juv",
year="2013", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds1a2 = predict(model, data.frame(days_since=ds, site="site1", age="adult",
year="2012", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds1a3 = predict(model, data.frame(days_since=ds, site="site1", age="adult",
year="2013", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds2j2 = predict(model, data.frame(days_since=ds, site="site2", age="juv",
year="2012", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds2j3 = predict(model, data.frame(days_since=ds, site="site2", age="juv",
year="2013", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds2a2 = predict(model, data.frame(days_since=ds, site="site2", age="adult",
year="2012", PC1=mean(PC1), arrival_date=mean(arrival_date)))
ds2a3 = predict(model, data.frame(days_since=ds, site="site2", age="adult",
year="2013", PC1=mean(PC1), arrival_date=mean(arrival_date)))

# plot
windows()
plot(x=ds, y=ds1j2, ylim=c(11000, 14500), type="l", lty=1,
ylab="predicted weight", xlab="days since 1st Sept")
points(range(ds), range(ds1j2), pch=5)
lines(x=ds, y=ds1j3, lty=1);  points(range(ds), range(ds1j3), pch=18)
lines(x=ds, y=ds1a2, lty=2);  points(range(ds), range(ds1a2), pch=5)
lines(x=ds, y=ds1a3, lty=2);  points(range(ds), range(ds1a3), pch=18)
lines(x=ds, y=ds2j2, lty=1);  points(range(ds), range(ds2j2), pch=1)
lines(x=ds, y=ds2j3, lty=1);  points(range(ds), range(ds2j3), pch=16)
lines(x=ds, y=ds2a2, lty=2);  points(range(ds), range(ds2a2), pch=1)
lines(x=ds, y=ds2a3, lty=2);  points(range(ds), range(ds2a3), pch=16)

legend("bottomright", lty=rep(1:2, 4), pch=c(5,18,5,18,1,16,1,16),
legend=c("site 1, juveniles, 2012", "site 1, juveniles, 2013", 