These data support Harrison et al (2017) Best Practice for Mixed Effects Models in Ecology and Evolution. This document contains all the code necessary to reproduce the simulations and Figure 2
Set number of simulations, and make two matrices to store the results of the simulations
nsims<-10000
#Data Storage
colmat.weak<-colmat.strong<-matrix(NA,ncol=2,nrow=nsims)
r.weak<-0.5
r.strong<-0.9
This section starts the loop that stores coefficient values for models where there is strong or weak collinearity between two predictors. The basic premise for these simulations, including how to specify correlations between predictors, are adapted from Freckleton (2011). This paper is definitely worth a read.
Freckleton RP (2011). Dealing with collinearity in behavioural and ecological data: model averaging and the problems of measurement error. Behavioral Ecology and Sociobiology, 65, 91-101.
for (k in 1:nsims){
##### Data Generation
#Collinear Models
x1=rnorm(100,0,2)
beta.x1<-1
alpha<-0
y=rnorm(100,alpha + beta.x1*x1,1)
#Weak
x2= r.weak * x1 +rnorm(100,0,(1-r.weak))
m1<-lm(y~x1+x2)
colmat.weak[k,1]<-m1$coefficients[2]
colmat.weak[k,2]<-m1$coefficients[3]
#Strong
x2.strong= r.strong * x1 +rnorm(100,0,(1-r.strong))
m1.strong<-lm(y~x1+x2.strong)
colmat.strong[k,1]<-m1.strong$coefficients[2]
colmat.strong[k,2]<-m1.strong$coefficients[3]
}
#Store All The Data
weak<-data.frame(vals=c(colmat.weak[,1],colmat.weak[,2],colmat.strong[,1],colmat.strong[,2]),group=rep(c("r = 0.5","r = 0.9"),each=2*nsims))
#Add Predictor IDs
weak$param<-rep(c("x1","x2"),each=nsims,2)
You will need these packages to make the plot: ggplot2 for plotting, RColorBrewer for its nice colour palettes, and cowplot because it has a version of the ‘ggsave’ function that is more flexible than the ggplot2 version.
library(ggplot2)
library(RColorBrewer)
library(cowplot)
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
##
## ggsave
Plot Code:
#Choose Palette Colours
colcols<-brewer.pal(5,"Pastel1")[1:2]
#Generate Plot
h1<-ggplot(weak,aes(x=vals,group=param,fill=param,colour=param)) + geom_density(alpha=0.5) + facet_grid(.~group)
h1 <- h1+ theme_bw() + scale_fill_manual(values=colcols,name="Predictor") + scale_colour_manual(values=colcols,guide=FALSE)
h1 <- h1 + ggtitle("Predictor Correlation") + xlab("Coefficient Values") + geom_segment(x=1,y=0,xend=1,yend=3.5,colour="black",size=1,linetype=2)
h1 <- h1 + theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),strip.text.x = element_text(size=20),legend.text = element_text(size=20),legend.title = element_text(size=20),axis.text=element_text(size=15),axis.title = element_text(size=15),plot.title = element_text(hjust = 0.5))
h1
#Save Plots as TIFF and PDF
#ggsave("Collinearity Plot v2.pdf",h1,width=10,height=5)
#ggsave("Collinearity Plot v2.tiff",h1,width=10,height=5)