Least Square Points for Factor Level Groups

Question

I am looking for a simple way to extract (and construct) the least squares of the means of the indicated combinations of levels of one factor for each level of another factor.

Sample data:

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

Setting custom contrasts for the time factor:

 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) 

Setting the model and calculating the least squares means:

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

Drawing:

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

Now the second chart has the combinations that I want, but it shows the difference between the combined tools, and not the tools themselves.

I can get individual values ​​from tissue.lsm$lsmeans and figure out the combined remedies myself, but I have a nagging feeling that there is an easier way that I just don’t see. All data must be in lsmobj , after all.

 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 

I look forward to your comments or suggestions. Thanks!

Answer 1

One option is to calculate the contrast ratios for the group tools of interest, in addition to the contrast ratios for the differences in group averages that you calculated in custom_contrasts . For example, you can do this separately as 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 

Here are just outlets for the liver , which are the group you are behind.

 ... 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 

If you know that you want both group tools and differences in the group level, you can simply add to the matrix of contrast coefficients that you create using 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)) 

And then name and execute the .lsmc function.

Answer 2

Following @aosmith's example:

 custom.means <- 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.means) <- c("early", "late", "very late") custom.means.lsmc <- function(...) return(custom.means) tissue.means <- confint(lsmeans(tissue.model, custom.means ~ time | tissue)$contrasts) library("ggplot2") p <- ggplot(tissue.means, aes(x = contrast, y = estimate, ymin = lower.CL, ymax = upper.CL)) + geom_errorbar() + facet_wrap(~ tissue, ncol = 4) + xlab("time") print(p)