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)



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


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

  • $\begingroup$ Unfortunately I still didn't solve the problem. I would be very very very thankful if someone could give me an hint on this or possible alternatives if possible... Thank you! $\endgroup$ Jul 15, 2015 at 22:01
  • $\begingroup$ From ?subset "Factors may have empty levels after subsetting; unused levels are not automatically removed. See droplevels for a way to drop all unused levels from a data frame." $\endgroup$
    – user20637
    Jul 16, 2015 at 10:05
  • $\begingroup$ Well, it's part of the way lsmeans is designed. See the documentation and the vignette where a reference grid is defined. If you truly want to compare just the cells that have data, then you are implicitly saying you are including the interaction between the factors. If you use that model, it will find one of the cells non-estimable and so set its lsmeans to NA (though the Tukey adjustment will still come out wrong). If you like the additive model, then that model is able to predict all 6 factor combinations, and your lsmeans are those predictions. $\endgroup$
    – Russ Lenth
    Jul 16, 2015 at 20:12

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

  • $\begingroup$ Thank you so much for your reply! I was using the additive model.. but now I can understand better what exactly I am doing. Thank you so much for your time!! It was really helpful.. :) $\endgroup$ Jul 23, 2015 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.