custom plot of lmer random intercepts and slopes-CodePudding

I would like to find a way to modify the solution provided by Didzis Elferts in response to this post.

I have a model with a random intercept and random slope, whereas in the original question was about a model with just a random intercept. I can't figure out how to modify the code from the answer for random intercepts and slopes.

The reason I am not using dotplot() or sjPlot to do all this automatically is that I wanted to construct a more customized plot, with the random intercepts and slopes plotted in a more organized manner, under facets that reflect the fixed effect groups to which the random groups belong. The only way I see being able to do this is extracting the random slopes, intercepts, and their variances and having my fixed effect categories as the other columns in my data frame. So alternatively if someone has an easy way to customize a dotplot or other type of plot to sort by fixed effect categories, I'm all ears!

I've tweaked the reproducible code to show what I mean:

Add some variable to Dyestuff dataset, which will be used as the random slope variable

library("lme4")
Dyestuff$variable <- rep(c("X","Y"), each = 15)

Run a random intercept and random slope model (instead of just random intercept model)

fit1 <- lmer(Yield ~ variable   (variable|Batch), Dyestuff)

the rest of the code is the same as in the previous post, with the exception that I used condVar=TRUE because postVar has been depreciated since the original question was asked; note though, that the attribute still needs to be named postVar

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

str(qq)

rand.interc<-randoms$Batch



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

with both random intercept and slope in qq, sd.interc doesn't work. I get this message:

Error in qq[, , 1:length(qq)] : subscript out of bounds

I've played around with it but can't seem to get it right.

CodePudding user response：

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)