Home › forums › Mixed Models › Printing tables with sjt.lmer function
Tagged: afex, mixed models, sjPlot, tables
 This topic has 10 replies, 2 voices, and was last updated 2 years, 7 months ago by henrik.

AuthorPosts


February 8, 2018 at 07:35 GMT+0000 #191Jarkko HautalaParticipant
I’ve managed to run my analyses with afex mixedfunction.
Now I want to print result tables and the most convenient way to do is to use sjt.lmer function from sjPlotpackage.
However, it seem not to find the TYPE III test pvalues but provides different pvalues and possibly confidence intervals (Wald).
Is this something which may be solved from afex or is solely the issue of sjPlot package handling the object returned from mixedfunction? 
February 8, 2018 at 08:50 GMT+0000 #192henrikKeymaster
I do not think
afex
andsjt.
orsjp.
functions are compatible. The main reason is that thesj
functions focus on the standardlmer
output and their fixedeffects parameters. In contrast,afex
focuses on tests of effects. In many cases, effects involve more than one parameter.However,
afex
involves many functions that support nice printing of the"mixed"
objects returned frommixed
. Mainlynice()
andanova()
. See the following example:library("afex") data("Machines", package = "MEMSS") # simple model with randomslopes for repeatedmeasures factor m1 < mixed(score ~ Machine + (MachineWorker), data=Machines) m1 # Mixed Model Anova Table (Type 3 tests, KRmethod) # # Model: score ~ Machine + (Machine  Worker) # Data: Machines # Effect df F p.value # 1 Machine 2, 4 32.80 ** .003 #  # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 nice(m1) # Mixed Model Anova Table (Type 3 tests, KRmethod) # # Model: score ~ Machine + (Machine  Worker) # Data: Machines # Effect df F p.value # 1 Machine 2, 4 32.80 ** .003 #  # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 knitr::kable(nice(m1)) # Effect df F p.value  # :::: # Machine 2, 4 32.80 ** .003  anova(m1) # Mixed Model Anova Table (Type 3 tests, KRmethod) # # Model: score ~ Machine + (Machine  Worker) # Data: Machines # num Df den Df F Pr(>F) # Machine 2 4 32.803 0.003302 ** #  # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 knitr::kable(anova(m1)) #   num Df den Df F Pr(>F) # ::::: # Machine  2 4 32.80302 0.0033024
Please let me know if this is not enough. If you definitely want to use the
sj
family of functions you need to contact the maintainer and ask if he is willing to supportafex
functions. 
February 8, 2018 at 09:04 GMT+0000 #193Jarkko HautalaParticipant
Thanks, I opened an issue also on sjPlot GIT page.
I believe I need to few other values, at least confidence intervals and variance explained for the fixed and random parts of the model.
What would be best way to get these?I’ve used:
>r.squaredGLMM(fit$full_model)
>confint(fit$full_model)but at least the confint seem to take a lot of computation time at least if continuous predictors are involved.

February 8, 2018 at 09:33 GMT+0000 #194henrikKeymaster
I do not know exactly which package
r.squaredGLMM
comes from, but I expect that, similar toconfint
, it also works on the level of the model parameters and not the level of the model effects (or ‘model terms’). So they are in principle incompatible withafex
. If in your case all effects only consist of one parameter, then you can of course use them.In such a case, the simplest way would be to add these values to the
data.frame
s returned by eithernice()
oranova()
. Continuing from the previous example:df1 < nice(m1) df1$new_col < 2 df1 # Mixed Model Anova Table (Type 3 tests, KRmethod) # # Model: score ~ Machine + (Machine  Worker) # Data: Machines # Effect df F p.value new_col # 1 Machine 2, 4 32.80 ** .003 2 #  # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 df2 < anova(m1) df2$new_col < 2 df2 # Mixed Model Anova Table (Type 3 tests, KRmethod) # # Model: score ~ Machine + (Machine  Worker) # Data: Machines # num Df den Df F Pr(>F) new_col # Machine 2 4 32.803 0.0033024 2
As you can see, both functions simply return a
data.frame
for which specificprint
methods exists. But you can simply add additional columns. If you want to remove the specificprint
methods, simply wrap the returned object in aas.data.frame()
call. For example:as.data.frame(df2) # num Df den Df F Pr(>F) new_col # Machine 2 4 32.80302 0.003302373 2
Note that you can also get confidence intervals via the methods implemented in
emmeans
, which work nicely withafex
functions, see: https://cran.rproject.org/package=emmeans 
February 8, 2018 at 09:40 GMT+0000 #195Jarkko HautalaParticipant
Ok, so the emmeans seems to be the way to go as I have a rather complex models to report, and in some models I have factors with three levels (=multiple parameters?).

February 8, 2018 at 09:48 GMT+0000 #196henrikKeymaster
In general, I would give up on the idea of having something like a standardized effect size on the level of the effects (i.e. model terms) in mixed models. They generally do not appear to work as soon as a model includes random slopes or multiple randomeffects grouping factors. Same applies for something like Rsquared.
Confidence intervals can be obtained on the level of the cell means (i.e., factor level or linear combinations thereof) using
emmeans
. For linear mixed modelsemm_options(lmer.df = ...)
allows you to choose how the degrees of freedoms are calculated which can have dramatic consequences on the speed. Three options exist:lmer.df = c("kenwardroger", "satterthwaite", "asymptotic")
.It is not always trivial how to combine the table of effects with a measure of effect size, but one can sometimes find ways. For example see Table 2 in: http://singmann.org/download/publications/2015_kellen_singmann_vogt_klauer.pdf

February 8, 2018 at 11:44 GMT+0000 #197Jarkko HautalaParticipant
Many thanks I was now being able to plot figures:
>R=ref_grid(fit, type=”response”, at = list(a = c(2, 1, 0, 1, 2)))
>A2=emmeans(R, specs=c(“a”,”b”, “c”))
>emmip(A2, b~ac, CIs=T)And also confidence intervals for a result table, e.g:
>R1=ref_grid(fit, type=”response”)
>A3=emmeans(R1, specs=c(“a”))Here a=continuous predictor, b, and c categorical predictors.
Its a pity if no variance explained or equivalent value cannot be computed, that would be highly appreciated by researchers.
Of course we get sense of how strong the predictors are by looking at the figures. 
February 8, 2018 at 12:33 GMT+0000 #198Jarkko HautalaParticipant
One thing more:
Is there a way to backtransform the values returned by summaryfunction for a mixed object?

February 8, 2018 at 12:42 GMT+0000 #199henrikKeymaster
I am sorry, but I do not understand what you want to backtransform into what. In general, a reproducible example (e.g., http://stackoverflow.com/q/5963269/289572) helps a lot in this regard. But you at least have to show the output you have and then tell me exactly what you want.

February 8, 2018 at 13:38 GMT+0000 #200Jarkko HautalaParticipant
>fit=mixed(log(A)~B+C+(1+Bid)+(1+Bitem), method=”S”, data=dat2, progress=TRUE, cl=cl)
>summary(fit)
here i would like to get regression estimates and random effect variance on the original values, not on logtransformedanother problem I’m encountering is that when I try to use emmeans emmip or summary functions it fails when the model includes the expand_re=TRUE specification:
>fit=mixed(log(A)~B+C+(1+Bid)+(1+Bitem), method=”S”, data=dat2, progress=TRUE, expand_re=TRUE, cl=cl)
>Grid=ref_grid(fit, at=list(C=c(2, 1, 0, 1, 2)))
>EMM=emmeans(Grid, specs=c(“C”, “B”))
>emmip(EMM, type=”response”, B~C, CIs=TRUE)Error in eval(predvars, data, env) : object ‘re1.B1’ not found

February 9, 2018 at 11:24 GMT+0000 #203henrikKeymaster
The output from
summary()
formixed()
objects is identical to the output ofsummary()
forlmer()
objects (this is why you get output on the level of the parameters and not at the level of model terms/effects). So if this is possible for normallmer()
models it is possible here as well (if it does make sense in your case is not something I can answer). If you want to get at the contents of the summary output, you can save it and work with it:summ < summary(fit) as.data.frame(summ$varcor) summ$coefficients
The problem when using
expand_re = TRUE
is that it creates a newdata.frame
that is used for fitting. When trying to access this, this fails. There are several ways around this.
1. You could avoid accessing it.emmeans
only tries to access it, if thelmer.df
are"satterthwaite"
:library("afex") data("Machines", package = "MEMSS") m1 < mixed(score ~ Machine + (Machine  Worker), data=Machines, expand_re = TRUE) emm_options(lmer.df = "kenwardroger") emmeans(m1, "Machine") #works # Machine emmean SE df lower.CL upper.CL # A 52.35556 2.397324 7.35 46.74087 57.97024 # B 60.32222 2.633283 9.04 54.36897 66.27548 # C 66.27222 2.836796 11.60 60.06746 72.47698 # # Degreesoffreedom method: kenwardroger # Confidence level used: 0.95 emm_options(lmer.df = "satterthwaite") emmeans(m1, "Machine") #fails # Error in eval(predvars, data, env) : object 're1.Machine1' not found emm_options(lmer.df = "asymptotic") emmeans(m1, "Machine") #works # Machine emmean SE df asymp.LCL asymp.UCL # A 52.35556 2.397324 Inf 47.65689 57.05422 # B 60.32222 2.633283 Inf 55.16108 65.48336 # C 66.27222 2.836796 Inf 60.71220 71.83224 # # Degreesoffreedom method: asymptotic # Confidence level used: 0.95
2. You can overwrite your existing
data.frame
with one that contains the missing columns (note that this might change the number of rows if you have missing data). I would only do this on a copy of my real data:m1 < mixed(score ~ Machine + (Machine  Worker), data=Machines, expand_re = TRUE) Machines < mixed(score ~ Machine + (Machine  Worker), data=Machines, expand_re = TRUE, return = "data") emm_options(lmer.df = "satterthwaite") emmeans(m1, "Machine") # works now # Machine emmean SE df lower.CL upper.CL # A 52.35556 2.397324 7.35 46.74094 57.97017 # B 60.32222 2.633283 9.04 54.36898 66.27546 # C 66.27222 2.836796 11.56 60.06544 72.47901 # # Degreesoffreedom method: satterthwaite # Confidence level used: 0.95


AuthorPosts
 You must be logged in to reply to this topic.