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)
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)
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"]]
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"]]
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
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"))
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")
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
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
Related Topics
Do I Need to Normalize (Or Scale) Data for Randomforest (R Package)
R Not Finding Package Even After Package Installation
Reading Hdf Files into R and Converting Them to Geotiff Rasters
Conditional Assignment of One Variable to the Value of One of Two Other Variables
Could Not Find Function Inside Foreach Loop
Creating Legend with Circles Leaflet R
Using a Static (Prebuilt) PDF Vignette in R Package
Knitr (R) - How Not to Embed Images in the HTML File
Ggplot Year by Year Comparison
How to Read \" Double-Quote Escaped Values with Read.Table in R
Shift Legend into Empty Facets of a Faceted Plot in Ggplot2
Marking Specific Tiles in Geom_Tile()/Geom_Raster()
Accessing Excel File from Sharepoint with R
Display and Save the Plot Simultaneously in R, Rstudio
How to Create 3D - Matlab Style - Surface Plots in R
What Does the Double Percentage Sign (%%) Mean
R Error in Unique.Default(X) Unique() Applies Only to Vectors