4

I have a question about running post hoc tests on linear mixed models:

I am running a linear mixed model in lme4 with 3 groups, 5 snakes per group, each group at a different ventilation rate (Vent), measurements taken at different time points (Time), with snake specified as a random effect (ID)

Data subset below:

subset1 <- structure(list(ID = structure(c(5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 
9L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 10L, 10L, 
10L, 10L, 10L, 20L, 20L, 20L, 20L, 20L, 4L, 4L, 4L, 4L, 4L, 11L, 
11L, 11L, 11L, 11L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 12L, 
13L, 14L, 15L, 16L, 17L, 17L, 17L, 17L, 17L), .Label = c("", 
"1_1_2", "10", "10_1_1", "13_1_4", "14_2_4", "15_3_4", "16_1_4", 
"17_2_4", "2_2_1", "5", "5_2_2", "5_2_3", "5_2_4", "5_2_5", "5_2_6", 
"7_1_2", "8", "9", "9_3_1"), class = "factor"), Vent = c(30L, 
30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 
30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 125L, 
125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 
125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 
125L, 125L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 
250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 
250L, 250L, 250L, 250L, 250L), Time = c(60L, 80L, 180L, 720L, 
1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 
60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 
80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 
180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 
720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 
1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 
60L, 80L, 180L, 720L, 1440L), corr.pO2 = c(224.1388673, 233.9456874, 
239.1553778, 107.2373336, 76.71835625, 164.6293748, 243.8501858, 
234.8205544, 71.74240501, 62.23789874, 69.69478654, 62.23789874, 
152.1142885, 79.61325688, 63.33285001, 240.8713061, 231.304842, 
222.7743953, 95.7912966, 64.41744793, 241.7255035, 238.2936023, 
138.1188987, 43.00663696, 50.64392111, 265.4973967, 274.0599252, 
285.0144919, 83.37647392, NA, 292.3660214, 281.6533627, 275.9747984, 
63.33285001, 56.59660394, 254.2521631, 222.3180596, 208.736288, 
88.83223104, 114.1782867, 208.255285, 232.1878564, 193.3861802, 
72.75355024, 60.01517133, 209.6956308, 245.9596884, 200.4342522, 
75.73874562, 67.61194011, 240.0146049, 261.1278627, 166.9318704, 
74.75152919, 73.75652657, 270.9724687, 251.7882317, 245.9596884, 
147.1396383, 50.64392111, 294.179467, 296.3431178, 284.6426934, 
73.75652657, 75.73874562, 233.0681297, 234.3834557, 143.3247511, 
73.75652657, 66.55672391, 245.9596884, 249.3041163, 223.6847954, 
92.35383362, 78.65544784)), .Names = c("ID", "Vent", "Time", 
"corr.pO2"), row.names = c(NA, 75L), class = "data.frame")

Code:

attach(subset1)

require(lme4)

with.vent = lmer(corr.pO2 ~ Vent * Time + (1|ID),REML = FALSE, data = subset1)

with.vent.no.int = lmer(corr.pO2 ~ Vent + Time + (1|ID),REML = FALSE, data = subset1)

anova(with.vent, with.vent.no.int)
#no significant interaction

Test for effect of vent:

without.vent = lmer(corr.pO2 ~ Time + (1|ID), REML = FALSE, data = subset1)

Compare with with vent:

anova(with.vent.no.int, without.vent)
#no significant effect of ventilation treatment p=0.09199

Test for effect of time:

without.time = lmer(corr.pO2 ~ Vent + (1|ID), data = subset1)

anova(with.vent.no.int, without.time)
# highly significant effect of time on pO2 < 2.2e-16 ***

So try post hoc test:

require(lsmeans)
lsmeans(with.vent.no.int, pairwise ~ Time, adjust = "tukey", data = subset1)

This is where I get the error which reads:

Error in solve.default(L %*% V0 %*% t(L), L) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

I can run pairwise tests with a correction using:

pairwise.t.test(corr.pO2, Time, p.adj = "BH", paired = T)

but know this will not work where other variables have an interaction (as is the case with my other data), and where I would like pairwise comparisons between each time point and ventilation regime. Is this possible with lsmeans()?

Thank you for your input, I know likelihood ratio tests are in themselves controversial. I have considered a mixed effects ANOVA, but there are some missing data points which make this impossible. The data was previously analysed by another student as a two way anova, no repeated measures, but my feeling was that this was inappropriate as each snakes was measured at repeated time points

Catherine Williams
  • 41
  • 2
  • 4
    I'm voting to close this question as off-topic because it is about R error messages without a reproducible example. – gung - Reinstate Monica Jan 07 '16 at 11:24
  • 1
    thank you @gung for the link, I have tried to add the data in a way that is reproducible. –  Jan 07 '16 at 13:09
  • The "l" in `lmer` stands for "linear" and a linearity assumption is often not really valid if time is modeled as a continuous predictor. – Roland Jan 08 '16 at 08:33

2 Answers2

4

The answer turns out to be fairly simple: you need to make sure that your Vent and Time predictors are factors. Otherwise lsmeans gets confused about the meaning of a pairwise test. (There is a slightly longer conversation to be had about whether you really wanted to analyze this model with continuous predictor, i.e. as a response surface design rather than 2-way ANOVA ...) Here's a slightly more compact version of your analysis:

subset1 <- transform(subset1,Vent=factor(Vent), Time=factor(Time))
require(lme4)
with.vent = lmer(corr.pO2 ~ Vent * Time + (1|ID),
       REML = FALSE, data = subset1)
drop1(with.vent,test="Chisq")  ## test interaction
with.vent.no.int = update(with.vent, . ~ . - Vent:Time)
drop1(with.vent.no.int,test="Chisq")  ## test main effects
require(lsmeans)
lsmeans(with.vent.no.int, pairwise ~ Time)

Subset of output:

$contrasts
 contrast    estimate       SE    df t.ratio p.value
 60 - 80     -6.99222 12.76886 63.45  -0.548  0.9819
 60 - 180    14.74281 12.76886 63.45   1.155  0.7768
 60 - 720   147.27139 12.76886 63.45  11.534  <.0001
...

I do agree that the error message is inscrutable. It might be worth mentioning this to the lsmeans maintainer to see if it's possible to detect and flag this (really very common) error.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • 2
    I'll take a look at it and see if I can trap it. Note the `adjust` and `data` arguments in the `lsmeans` call are unnecessary. – Russ Lenth Jan 08 '16 at 02:49
  • I think I would make at least `Time` an ordered factor, probably `Vent` too (although that shouldn't impact the `drop1` results). – Roland Jan 08 '16 at 08:32
  • Thank you very much @Ben Bolker, this was really helpful. The question on whether I am taking the analysis the right way would also be useful for me, I had thought that a 2 way ANOVA was inappropriate as repeated samples were taken from the same snakes, but also that repeated measures anova would not work as there are a small number of missing data points. I would welcome your opinion on this, as I am trying to become more stats literate. – Catherine Williams Jan 08 '16 at 10:05
  • Thank you @Roland I am looking into this – Catherine Williams Jan 08 '16 at 10:05
  • By "2-way ANOVA" I just meant a linear mixed model with two fixed categorical predictors (as opposed to two continuous predictors) – Ben Bolker Jan 08 '16 at 14:19
  • Ah OK, thanks @ben bolker, I will read some more about response surface designs. – Catherine Williams Jan 08 '16 at 17:27
3

The next update for lsmeans (probably around Feb 1, 2016) will catch this kind of error:

> lsmeans(with.vent.no.int, pairwise ~ Vent)

$lsmeans
     Vent   lsmean       SE    df lower.CL upper.CL
 135.1351 167.4871 6.859275 18.63 153.1111  181.863

Confidence level used: 0.95 

$contrasts


Warning message:
In contrast.ref.grid(result, method = contr, by, ...) :
  No contrasts were generated! Perhaps only one lsmean is involved.
  This can happen, for example, when your predictors are not factors.

The ref.grid function is handy for understanding what you've got:

> ref.grid(with.vent.no.int)
'ref.grid' object with variables:
    Vent = 135.14
    Time = 483.24

Both Vent and Time are covariates, so by default their averages are used. To change this, you do not necessarily have to alter the dataset; you can just coerce the predictors to factors in the model:

> repaired = lmer(corr.pO2 ~ factor(Vent) + factor(Time) + (1|ID), 
                  REML = FALSE, data = subset1)
> ref.grid(repaired)
'ref.grid' object with variables:
    Vent =  30, 125, 250
    Time =   60,   80,  180,  720, 1440

> lsmeans(repaired, pairwise ~ Vent)
$lsmeans
 Vent   lsmean       SE    df lower.CL upper.CL
   30 146.0967 12.19373 18.16 120.4952 171.6981
  125 177.0917 12.29417 18.66 151.3274 202.8559
  250 173.2568 11.12879 26.72 150.4111 196.1024

Results are averaged over the levels of: Time 
Confidence level used: 0.95 

$contrasts
 contrast    estimate       SE    df t.ratio p.value
 30 - 125  -30.994975 17.31570 18.41  -1.790  0.2005
 30 - 250  -27.160077 16.50870 21.52  -1.645  0.2490
 125 - 250   3.834898 16.58302 21.81   0.231  0.9710

Results are averaged over the levels of: Time 
P value adjustment: tukey method for comparing a family of 3 estimates
Russ Lenth
  • 5,922
  • 2
  • 13
  • 21