lsmeans analysis is creating data that does not exist in my input file

Question

Good afternoon,

I am trying to see if the expression of a given gene (conc_ul) is explained by 3 categories (Sample, background and Experiment).

Initially I am trying to check which is the most simple model to explain my data by doing the following:

---

m1<-glm(conc_ul ~ Sample*background+Experiment+Experiment:Sample+Experiment:background,family="quasipoisson", data=mydata)

m2<-glm(conc_ul ~ Sample*background+Experiment+Experiment:Sample,family="quasipoisson", data=mydata)

anova(m1,m2,test="Chisq")

(...)

and so on until I cannot simplify deeper my model

---

Once I have the categories and interactions that explain the conc_ul (simpler model) I wanted to do the multiple comparisons between the levels of each category using:

---

(m6 is the more simplified model for a given gene)

m6<-glm(conc_ul ~ Sample+Experiment,family="quasipoisson", data=mydata)

lsm_m6 <- lsmeans(m6, ~ Sample | Experiment, test="Tukey")

plot(lsm_m6)

cld(lsm_m6, alpha = 0.01)

---

Here it comes the problem.

My "Sample" category has 3 levels; "Experiment" has 2. Only 1 of the experiments has the 3 levels of Sample, the other one has 2.

For some of my genes when I do the multiple comparisons using lsmeans and cld the 3 levels of Sample appear for BOTH levels of "Experiment". I already checked if there is something wrong in my input file which is not the case.

I don't know what I am doing wrong or if I am using incorrectly lsmeans for this analysis.

I would really appreciate if someone could help me with this! :)

Many thanks and all the best, Marta

Answer 1

I'll attempt to explain my comment further, via one of the example datasets in R. The warpbreaks dataset has data on 2 wools and 3 tensions. I'll subset the data so that there are only two tensions with the second wool. This is similar to the situation you describe. I'll fit a model; and obtain the lsmeans:

> wb = warpbreaks[1:45,]
> wb.wt = lm(breaks ~ wool*tension, data = wb)

> library(lsmeans)
> lsmeans(wb.wt, ~ tension | wool)
wool = A:
 tension   lsmean       SE df lower.CL upper.CL
 L       44.55556 3.927664 40 36.61745 52.49366
 M       24.00000 3.927664 40 16.06190 31.93810
 H       24.55556 3.927664 40 16.61745 32.49366

wool = B:
 tension   lsmean       SE df lower.CL upper.CL
 L       28.22222 3.927664 40 20.28412 36.16033
 M       28.77778 3.927664 40 20.83967 36.71588
 H             NA       NA NA       NA       NA

Confidence level used: 0.95 

The case with wool B and tension H has NAs because there are no data in that cell. lsmeans always constructs a regular grid of the factor combinations, viewed as crossed factors.

If you don't want a grid like that, you can use the interaction function to construct a single factor having as its levels the factor combinations actually observed, and fit a model using that factor:

> wb = transform(wb, trt = interaction(wool, tension, lex.order = TRUE))
> wb.trt = lm(breaks ~ trt, data = wb)
> lsmeans(wb.trt, ~ trt)
 trt   lsmean       SE df lower.CL upper.CL
 A.L 44.55556 3.927664 40 36.61745 52.49366
 A.M 24.00000 3.927664 40 16.06190 31.93810
 A.H 24.55556 3.927664 40 16.61745 32.49366
 B.L 28.22222 3.927664 40 20.28412 36.16033
 B.M 28.77778 3.927664 40 20.83967 36.71588

Confidence level used: 0.95

Note that now we only have the five means where data actually exist. If you do pairwise comparisons of these, there will be ${5 \choose 2}=10$ comparisons. But ifyou want comparisons only within levels of wool, you need custom contrasts, like this:

> con = list(`M-L|A` = c(-1,1,0, 0,0), `H-L|A` = c(-1,0,1, 0,0), 
+            `H-M|A` = c(0,-1,1, 0,0), `M-L|B` = c(0,0,0, -1,1))

> contrast(lsmeans(wb.trt, ~ trt), con, adjust = "mvt")
 contrast    estimate       SE df t.ratio p.value
 M.L.A    -20.5555556 5.554555 40  -3.701  0.0025
 H.L.A    -20.0000000 5.554555 40  -3.601  0.0033
 H.M.A      0.5555556 5.554555 40   0.100  0.9996
 M.L.B      0.5555556 5.554555 40   0.100  0.9996

P value adjustment: mvt method for 4 tests

I chose the "mvt" (multivariate $t$) method for adjusting these. It accounts for the irregular correlation structure associated with these selected contrasts.

(I'm not sure why the labels got altered in the output. I'll try to fix that in the next update.)