Plot Random Effects from Lmer (Lme4 Package) Using Qqmath or Dotplot: How to Make It Look Fancy

Plot random effects from lmer (lme4 package) using qqmath or dotplot: How to make it look fancy?

One possibility is to use library ggplot2 to draw similar graph and then you can adjust appearance of your plot.

First, ranef object is saved as randoms. Then variances of intercepts are saved in object qq.

randoms<-ranef(fit1, postVar = TRUE)
qq <- attr(ranef(fit1, postVar = TRUE)[[1]], "postVar")

Object rand.interc contains just random intercepts with level names.

rand.interc<-randoms$Batch

All objects put in one data frame. For error intervals sd.interc is calculated as 2 times square root of variance.

df<-data.frame(Intercepts=randoms$Batch[,1],
sd.interc=2*sqrt(qq[,,1:length(qq)]),
lev.names=rownames(rand.interc))

If you need that intercepts are ordered in plot according to value then lev.names should be reordered. This line can be skipped if intercepts should be ordered by level names.

df$lev.names<-factor(df$lev.names,levels=df$lev.names[order(df$Intercepts)])

This code produces plot. Now points will differ by shape according to factor levels.

library(ggplot2)
p <- ggplot(df,aes(lev.names,Intercepts,shape=lev.names))

#Added horizontal line at y=0, error bars to points and points with size two
p <- p + geom_hline(yintercept=0) +geom_errorbar(aes(ymin=Intercepts-sd.interc, ymax=Intercepts+sd.interc), width=0,color="black") + geom_point(aes(size=2))

#Removed legends and with scale_shape_manual point shapes set to 1 and 16
p <- p + guides(size=FALSE,shape=FALSE) + scale_shape_manual(values=c(1,1,1,16,16,16))

#Changed appearance of plot (black and white theme) and x and y axis labels
p <- p + theme_bw() + xlab("Levels") + ylab("")

#Final adjustments of plot
p <- p + theme(axis.text.x=element_text(size=rel(1.2)),
axis.title.x=element_text(size=rel(1.3)),
axis.text.y=element_text(size=rel(1.2)),
panel.grid.minor=element_blank(),
panel.grid.major.x=element_blank())

#To put levels on y axis you just need to use coord_flip()
p <- p+ coord_flip()
print(p)

Sample Image

custom plot of lmer random intercepts and slopes

I think if you want the standard deviations of the conditional modes for the intercepts you should use sqrt(qq[1,1,]); this extracts the (1,1) element of each "slice" of the conditional-variance array.

df <- data.frame(Intercepts=randoms$Batch[,1],
sd.interc=2*sqrt(qq[1,1,]),
lev.names=rownames(rand.interc))

(I'm not sure why you doubled these values, I guess as a preliminary to using them as half-widths of a confidence bar?)

However, there are easier tools to use now.

broom.mixed::tidy(fit1, "ran_vals", conf.int = TRUE)

gives

# A tibble: 12 × 8
effect group level term estimate std.error conf.low conf.high
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 ran_vals Batch A (Intercept) -12.2 14.8 -41.2 16.8
2 ran_vals Batch B (Intercept) -1.93 14.8 -30.9 27.0
3 ran_vals Batch C (Intercept) 14.1 14.8 -14.9 43.1
...

or as.data.frame(ranef(fit1)) (this is even built into lme4)

Dotplot of two random effects in one graph

The problem you have is that ranef() stores the result of your analysis as a list, not data frame. You are trying to access this list as you would access data frame, that's why it doesn't work. The trick is to use double square bracket with dotplot to access the list. Then you can use grid.arrange to quickly combine your plots (or you can use @juba solution).

library(lme4)
fm3 <- lmer(strength ~ 1 + (1 | batch) + (1 | sample), Pastes)
d1 <- ranef(fm3, postVar=TRUE)
#double square bracket access the lists in d1
a <- dotplot(d1)[["batch"]]
b <- dotplot(d1)[["sample"]]
#grid.arrange combines your plots
library(gridExtra)
grid.arrange(a,b, nrow=1)

Sample Image

How can I sort random effects by value of the random effect (not the intercept) in dotplot or ggplot2

In the part of the function where ord and pDf are created the rows are reordered by size of the intercept. We can make that part optional like in the updated function below:

require(ggplot2)
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE, reorder=TRUE) {
require(ggplot2)
f <- function(x) {
pv <- attr(x, "postVar")
cols <- 1:(dim(pv)[1])
se <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
if (reorder) {
ord <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x))
pDf <- data.frame(y=unlist(x)[ord],
ci=1.96*se[ord],
nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
ID=factor(rep(rownames(x), ncol(x))[ord], levels=rownames(x)[ord]),
ind=gl(ncol(x), nrow(x), labels=names(x)))
} else {
pDf <- data.frame(y=unlist(x),
ci=1.96*se,
nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
ID=factor(rep(rownames(x), ncol(x)), levels=rownames(x)),
ind=gl(ncol(x), nrow(x), labels=names(x)))
}

if(QQ) { ## normal QQ-plot
p <- ggplot(pDf, aes(nQQ, y))
p <- p + facet_wrap(~ ind, scales="free")
p <- p + xlab("Standard normal quantiles") + ylab("Random effect quantiles")
} else { ## caterpillar dotplot
p <- ggplot(pDf, aes(ID, y)) + coord_flip()
if(likeDotplot) { ## imitate dotplot() -> same scales for random effects
p <- p + facet_wrap(~ ind)
} else { ## different scales for random effects
p <- p + facet_grid(ind ~ ., scales="free_y")
}
p <- p + xlab("Levels") + ylab("Random effects")
}

p <- p + theme(legend.position="none")
p <- p + geom_hline(yintercept=0)
p <- p + geom_errorbar(aes(ymin=y-ci, ymax=y+ci), width=0, colour="black")
p <- p + geom_point(aes(size=1.2), colour="blue")
return(p)
}

lapply(re, f)
}

If we now call the function with reorder = FALSE we get the original order of the rows:

ggCaterpillar(ranef(fit,condVar=TRUE), QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]]

Sample Image

You can reorder rows before plotting if you'd like, for example to reverse the order:

ref <- ranef(fit,condVar=TRUE)
ref$Subject <- ref$Subject[nrow(ref$Subject):1, ]
ggCaterpillar(ref, QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]]

Sample Image

Problem combining ranef plots of two models with grid.arrange()

(A reproducible example would be useful, I hope this example does what you want ...). I think the key here is to recognize that the dotplot.ranef.mer method returns a list of plots:

library(lme4)
fm1 <- lmer(angle ~ (1|recipe) + (1|recipe:replicate), cake, REML= FALSE)
dd <- dotplot(ranef(fm1))
length(dd) ## 2

They're not necessarily in the same order as in the formula:

names(dd) ## [1] "recipe:replicate" "recipe"          
print(dd[["recipe"]])
print(dd[["recipe:replicate"]])

So you would want something like

f <- function(m) dotplot(ranef(m))[["individual"]]
gridExtra::grid.arrange(f(lme1),f(lme2))

Plot Selected Random Effects Observations in Lattice

We need to manipulate the ranef.mer object.

library(lme4)
library(lattice)
data(sleepstudy)
fit <- lmer(Reaction ~ Days + (1 + Days|Subject), data=sleepstudy)

First, we store it.

r.int <- ranef(fit, condVar=TRUE)

Second, we create a vector of the desired subset row numbers.

s <- c(337, 310, 333, 349)

Third, inside a lapply function we subset both in the list, the data.frame and importantly the attributes, where the variances are hidden in an array.

r.int <- lapply(r.int, function(x) {
s2 <- which(rownames(x) %in% s)
x <- x[s2, ]
attributes(x)$postVar <- attributes(x)$postVar[, , s2]
return(x)
})

Fourth we hack the required class label.

class(r.int) <- "ranef.mer"

Et voilà, we finally can plot our desired selection.

dotplot(r.int)

Yields
Sample Image

How to plot some terms of a lmer model

For some reason, the terms option doesn't work. Should be reflected to the authors of sjplot. Here's two workarounds:

1) manually define the terms using ggplot2:

library(ggplot2)
library(sjPlot)

plot_model(model, type = "re") + scale_x_discrete(limits=c("A","D"))

Sample Image

You get the following warnings because you throw out data

Scale for 'x' is already present. Adding another scale for 'x', which
will replace the existing scale. Warning messages: 1: Removed 3 rows
containing missing values (geom_point). 2: Removed 3 rows containing
missing values (geom_errorbar).

2) Plot it from get_model()

    df = get_model_data(model,type="re")
df = subset(df,term %in% c("A","D"))
ggplot(df,aes(x=term,y=estimate,col=group)) +
geom_point(show.legend=FALSE) +
geom_segment(aes(xend=term,y=conf.low,yend=conf.high),show.legend=FALSE)+
scale_colour_brewer(palette = "Set1")+
coord_flip()+ggtitle("Random Effects")

Sample Image

How to only show fixed effect estimates of lmer model using sjPlot::plot_model

Have you tried reinstalling all the packages and restarting your RStudio session? It's possible that you overwrote some sort of default because I can't reproduce the images that you have attached. I'm using the most recent version of sjPlot from CRAN (2.8.9).

This is what I get when I run your code, slightly modified to call the plot to plot window:

library(lme4)
data(iris)
lmerMod <- lmer(Sepal.Length ~ Sepal.Width * Petal.Length + (1|Species), data=iris)

library(sjPlot)
plot1 <- plot_model(lmerMod,
type = "est",
show.intercept = TRUE,
title = "Estimates of fixed effects")
plot1

Sample Image

And then I run the second part:

glmMod <- glm(Sepal.Length ~ Sepal.Width * Petal.Length, data=iris)
plot2 <- plot_model(glmMod,
type = "est",
show.intercept = TRUE,
title = "Estimates of fixed effects")
plot2

Sample Image



Related Topics



Leave a reply



Submit