r - 绘制因子水平组的最小二乘均值

Question

我正在寻找一种简单的方法来提取(并绘制)一个因素水平的指定组合对于另一个因素的每个水平的最小二乘均值。

示例数据:

set.seed(1)
model.data <- data.frame(time = factor(paste0("day", rep(1:8, each = 16))),
  animal = factor(rep(1:16, each = 8)),
  tissue = factor(c("blood", "liver", "kidney", "brain")),
  value = runif(128)
  )

为因素“时间”设置自定义对比:

library("phia")
custom.contrasts <- as.data.frame(contrastCoefficients(
   time ~ (day1+day2+day3)/3 - (day4+day5+day6)/3,
   time ~ (day1+day2+day3)/3 - (day7+day8)/2,
   time ~ (day4+day5+day6)/3 - (day7+day8)/2,
   data = model.data, normalize = FALSE))

colnames(custom.contrasts) <- c("early - late",
  "early - very late",
  "late  - very late")

custom.contrasts.lsmc <- function(...) return(custom.contrasts)

拟合模型并计算最小二乘均值:

library("lme4")
tissue.model <- lmer(value ~ time * tissue + (1|animal), model.data)
library("lsmeans")
tissue.lsm <- lsmeans(tissue.model, custom.contrasts ~ time | tissue)

绘图:

plot(tissue.lsm$lsmeans)
dev.new()
plot(tissue.lsm$contrasts)

现在，第二个图具有我想要的组合，但它显示了组合平均值之间的差异，而不是平均值本身。

我可以从 tissue.lsm$lsmeans 获取各个值并自己计算组合平均值，但我总觉得有一种我看不到的更简单的方法。毕竟，所有数据都应该位于 lsmobj 中。

early.mean.liver = mean(model.data$value[model.data$tissue == "liver" & 
  model.data$time %in% c("day1", "day2", "day3")])
late.mean.liver = mean(model.data$value[model.data$tissue == "liver" & 
  model.data$time %in% c("day4", "day5", "day6")])
vlate.mean.liver = mean(model.data$value[model.data$tissue == "liver" & 
  model.data$time %in% c("day7", "day8")])
# ... for each level of "tissue"


#compare to tissue.lsm$contrasts
early.mean.liver - late.mean.liver 
early.mean.liver - vlate.mean.liver
late.mean.liver - vlate.mean.liver

我期待听到您的意见或建议。谢谢!

Answer 1

另一种替代方法是，除了在 custom_contrasts 中计算的组均值差异的对比系数之外，还计算感兴趣的组均值的对比系数。例如，您可以单独执行此操作:custom.contrasts2。

custom.contrasts2 <- as.data.frame(contrastCoefficients(
    time ~ (day1+day2+day3)/3,
    time ~ (day4+day5+day6)/3,
    time ~ (day7+day8)/2,
    data = model.data, normalize = FALSE))

colnames(custom.contrasts2) <- c("early",
                          "late",
                          "very late")

custom.contrasts2.lsmc <- function(...) return(custom.contrasts2)

lsmeans(tissue.model, custom.contrasts2 ~ time | tissue)$contrasts

这里只是 liver 的输出，这是您所追求的组。

...
 tissue = liver:
 contrast    estimate         SE   df t.ratio p.value
 early      0.4481244 0.07902715 70.4   5.671  <.0001
 late       0.4618041 0.07902715 70.4   5.844  <.0001
 lvery late 0.3824247 0.09678810 70.4   3.951  0.0002

如果您知道需要组均值和组均值差异，则只需添加到通过contrastCoefficients 创建的对比系数矩阵即可。

custom.contrasts <- as.data.frame(contrastCoefficients(
    time ~ (day1+day2+day3)/3,
    time ~ (day4+day5+day6)/3,
    time ~ (day7+day8)/2,
    time ~ (day1+day2+day3)/3 - (day4+day5+day6)/3,
    time ~ (day1+day2+day3)/3 - (day7+day8)/2,
    time ~ (day4+day5+day6)/3 - (day7+day8)/2,
    data = model.data, normalize = FALSE))

然后相应地命名并创建 .lsmc 函数。