0

我正在使用 sjPlot(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")

非常感谢

可重现的例子:

#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)

在此处输入图像描述

4

1 回答 1

2

使用sjPlot,目前是不可能的。但是,我编写了一个专门用于计算和绘制边际效应的包:ggeffects。这个包更灵活一些(用于边际效应图)。

ggeffects -package中有一个-function 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中有一些示例,请特别参阅本节

编辑

以下是基于您的可重现示例的结果(请注意,当前需要 GitHub-Version!):

# 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-图

在此处输入图像描述

在此处输入图像描述

于 2017-04-30T20:24:44.460 回答