2017-04-13 125 views
1

我对4个时间点的个体进行了纵向重复测量。按照固定效应和随机斜率的时间混合模型分析,我已经使用平均值来估计每个时间点的平均值以及95%的置信区间。我现在想绘制带有时间点(x)和结果变量(y)与CI的平均值的线图。我可以使用例如ggplot来绘制我从lsmeans得到的结果?还是有另一种巧妙的方式来绘制这个?混合模型/ lsmeans结果的线图(使用ggplot?)

,我从LSMEANS得到,而我想积(lsmean,lower.CL,upperCL随着时间的推移)的结果是:

$lsmeans 
time lsmean  SE df lower.CL upper.CL 
0 21.967213 0.5374422 60 20.892169 23.04226 
1 16.069586 0.8392904 60 14.390755 17.74842 
2 13.486802 0.8335159 60 11.819522 15.15408 
3  9.495137 0.9854642 60 7.523915 11.46636 

Confidence level used: 0.95 

回答

0

这是你的意思?

# To convert from lsmeans output (d <- lsmeans(paramaters)) 
d <- summary(d)$lsmeans[c("lsmean", "lower.CL", "upper.CL")] 

library(ggplot2) 
ggplot(d, aes(time)) + 
    geom_line(aes(y = lsmean)) + 
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), 
        width = 0.2) + 
    geom_point(aes(y = lsmean), size = 3, 
       shape = 21, fill = "white") + 
    labs(x = "Time", y = "ls mean", 
     title = "ls mean result over time") + 
    theme_bw() 

Dirty solution

+0

是的,这也正是它。如何在上面的代码中将正确的数据转化为d?我已经尝试了 d < - summary(lsmeans(model,pairwise〜time,adjust =“tukey”)) 但是这会返回到ggplot无法使用的列表。 –

+0

试试这个:'library(broom); d < - tidy(lsmeans(model,pairwise〜time,adjust =“tukey”))' – PoGibas

+0

整洁的命令给我一个错误信息,无法识别这个列表 –

0

总之,整个代码,这将使你的估计和混合模型的情节是:

## random slope model 
summary(model <- lme(outcome ~ time, random = ~1+time|ID, data = data, 
na.action = na.exclude, method = "ML")) 

## pairwise comparisons of timepoints 
install.packages("lsmeans") 
library(lsmeans) 
lsmeans(model, pairwise~time, adjust="tukey") 

### Draw the picture 
d <- summary(lsmeans(model, ~time)) 

library(ggplot2) 
ggplot(d, aes(time)) + 
    geom_line(aes(y = lsmean, group = 1)) + 
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) + 
    geom_point(aes(y = lsmean), size = 3, shape = 21, fill = "white") + 
    labs(x = "Time", y = "ls mean", title = "ls mean result over time") + 
    theme_bw()