Skip to content

Instantly share code, notes, and snippets.

@ibartomeus
Created January 21, 2015 14:31
Show Gist options
  • Save ibartomeus/f493bf142073c704391b to your computer and use it in GitHub Desktop.
Save ibartomeus/f493bf142073c704391b to your computer and use it in GitHub Desktop.
This shows how to get the random slopes and CI's for each level in a hierarchical model
#This shows how to get the random slopes and CI's for each level in a hierarchical model
#dataset used
head(iris)
#what we want to investigate
#Is there a general relationship? and how it differs by species
plot(iris$Sepal.Width ~ iris$Petal.Width, col = iris$Species, las =1)
#Our model with random slope and intercept
library(lmer)
m2 <- lmer(data = iris, Sepal.Width ~ Petal.Width + (1 + Petal.Width|Species))
summary(m2)
#extract fixed effects
a=fixef(m2)
a
#extract random effects
b=ranef(m2, condVar=TRUE)
b
# Extract the variances of the random effects
qq <- attr(b[[1]], "postVar")
qq
e=(sqrt(qq))
e=e[2,2,] #here we want to access the Petal.Weigth, which is stored in column 2 in b[[1]].
e
#calculate CI's
liminf=(b[[1]][2]+a[2])-(e*2)
liminf
mean_=(b[[1]][2]+a[2])
mean_
limsup=(b[[1]][2]+a[2])+(e*2)
limsup
#Plot betas and its errors
dotchart(mean_$Petal.Width,
labels = rownames(mean_), cex = 0.5,
xlim = c(0.4,1.4),
xlab = "betas")
#add CI's...
for (i in 1:nrow(mean_)){
lines(x = c(liminf[i,1],
limsup[i,1]), y = c(i,i))
}
#make final plot
plot(iris$Sepal.Width ~ iris$Petal.Width, col = iris$Species, las = 1)
#and plot each random slope
abline(a = b[[1]][1,1]+a[1], b= mean_$Petal.Width[1], col = "black")
abline(a = b[[1]][2,1]+a[1], b= mean_$Petal.Width[2], col = "red")
abline(a = b[[1]][3,1]+a[1], b= mean_$Petal.Width[3], col = "green")
#and general response
abline(a, lty = 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment