sjPlotのやりとりにおける平均モデレータオプションの適合

Question

私は、sjpPlot、sjp.int関数を使用して、lmeのインタラクションをプロットしています。 モデレータ値のオプションは、平均+/- sd、四分位数、すべて、最大/最小です。 +/- 2sdの平均をプロットする方法はありますか？

通常、それはこのようになります：

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

感謝

再現例：sjPlotで

#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

、それは現時点ではできません。しかし、私は特別な計算を行い、限界効果をプロットするパッケージを書いています：ggeffects。このパッケージは、（マージンエフェクトプロット用に）少し柔軟です。

ggeffects -packageには、特定の値で限界効果を計算できるggpredict()関数があります。あなたが問題になっているモデルの用語のSDを知ったら、あなたの相互作用をプロットする関数呼び出しでこれらの値を指定することができます。

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) 

いくつかの例がpackage-vignetteであり、特にthis sectionを参照してください。ここを編集

はあなたの再現性の例をもとに、結果である（GitHubの-バージョンが現在必要とされていることをノート！）：

# 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-プロット

# 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プロット

# 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プロット