Adapting the meansd moderator option in sjPlot interaction

Question

I am using sjPlot, the sjp.int function, to plot an interaction of an lme. The options for the moderator values are means +/- sd, quartiles, all, max/min. Is there a way to plot the mean +/- 2sd?

Typically it would be like this:

 model <- lme(outcome ~ var1+var2*time, random=~1|ID, data=mydata, na.action="na.omit")
 sjp.int(model, show.ci=T, mdrt.values="meansd")

Many thanks

Reproducible example:

#create data
mydata <- data.frame( SID=sample(1:150,400,replace=TRUE),age=sample(50:70,400,replace=TRUE), sex=sample(c("Male","Female"),200, replace=TRUE),time= seq(0.7, 6.2, length.out=400), Vol =rnorm(400),HCD =rnorm(400))  
mydata$time <- as.numeric(mydata$time)

 #insert random NAs
  NAins <-  NAinsert <- function(df, prop = .1){
n <- nrow(df)
m <- ncol(df)
num.to.na <- ceiling(prop*n*m)
id <- sample(0:(m*n-1), num.to.na, replace = FALSE)
rows <- id %/% m + 1
cols <- id %% m + 1
sapply(seq(num.to.na), function(x){
    df[rows[x], cols[x]] <<- NA
}
)
return(df)
}


mydata2 <- NAins(mydata,0.1)

#run the lme which gives error message
model = lme(Vol ~ age+sex*time+time* HCD, random=~time|SID,na.action="na.omit",data=mydata2);summary(model)

mydf <- ggpredict(model, terms=c("time","HCD [-2.5, -0.5, 2.0]"))

#lmer works
 model2 = lmer(Vol ~ age+sex*time+time* HCD+(time|SID),control=lmerControl(check.nobs.vs.nlev = "ignore",check.nobs.vs.rankZ = "ignore", check.nobs.vs.nRE="ignore"), na.action="na.omit",data=mydata2);summary(model)
 mydf <- ggpredict(model2, terms=c("time","HCD [-2.5, -0.5, 2.0]"))

#plotting gives problems (jittered lines)
plot(mydf)

Answer 1

With sjPlot, it's currently not possible. However, I have written a package especially dedicated to compute and plot marginal effects: ggeffects. This package is a bit more flexible (for marginal effects plots).

In the ggeffects-package, there's a ggpredict()-function, where you can compute marginal effects at specific values. Once you know the sd of your model term in question, you can specify these values in the function call to plot your interaction:

library(ggeffects)
# plot interaction for time and var2, for values
# 10, 30 and 50 of var2
mydf <- ggpredict(model, terms = c("time", "var2 [10,30,50]"))
plot(mydf)

There are some examples in the package-vignette, see especially this section.

Edit

Here are the results, based on your reproducible example (note that GitHub-Version is currently required!):

# requires at least the GitHub-Versiob 0.1.0.9000!
library(ggeffects)
library(nlme)
library(lme4)
library(glmmTMB)

#create data
mydata <-
  data.frame(
    SID = sample(1:150, 400, replace = TRUE),
    age = sample(50:70, 400, replace = TRUE),
    sex = sample(c("Male", "Female"), 200, replace = TRUE),
    time = seq(0.7, 6.2, length.out = 400),
    Vol = rnorm(400),
    HCD = rnorm(400)
  )
mydata$time <- as.numeric(mydata$time)

#insert random NAs
NAins <-  NAinsert <- function(df, prop = .1) {
  n <- nrow(df)
  m <- ncol(df)
  num.to.na <- ceiling(prop * n * m)
  id <- sample(0:(m * n - 1), num.to.na, replace = FALSE)
  rows <- id %/% m + 1
  cols <- id %% m + 1
  sapply(seq(num.to.na), function(x) {
    df[rows[x], cols[x]] <<- NA
  })
  return(df)
}

mydata2 <- NAins(mydata, 0.1)

# run the lme, works now
model = lme(
  Vol ~ age + sex * time + time * HCD,
  random =  ~ time |
    SID,
  na.action = "na.omit",
  data = mydata2
)
summary(model)

mydf <- ggpredict(model, terms = c("time", "HCD [-2.5, -0.5, 2.0]"))
plot(mydf)

lme-plot

# lmer also works
model2 <- lmer(
  Vol ~ age + sex * time + time * HCD + (time |
                                           SID),
  control = lmerControl(
    check.nobs.vs.nlev = "ignore",
    check.nobs.vs.rankZ = "ignore",
    check.nobs.vs.nRE = "ignore"
  ),
  na.action = "na.omit",
  data = mydata2
)
summary(model)
mydf <- ggpredict(model2, terms = c("time", "HCD [-2.5, -0.5, 2.0]"), ci.lvl = NA)

# plotting works, but only w/o CI
plot(mydf)

lmer-plot

# lmer also works
model3 <- glmmTMB(
  Vol ~ age + sex * time + time * HCD + (time | SID),
  data = mydata2
)
summary(model)
mydf <- ggpredict(model3, terms = c("time", "HCD [-2.5, -0.5, 2.0]"))
plot(mydf)
plot(mydf, facets = T)

glmmTMB-plots