Does this mixed model violate assumptions of independence?

A disturbance event caused damage to 5 streams (Set1). To quantify this damage, five additional unimpacted streams (Set2) were picked for comparison. During the selection process every effort was made for Set2 to be identical in every way except construction damage. Then 6 years later, a second, larger disturbance event affected both stream sets so they were resampled.

I am using mixed models to measure how animals have responded to the two disturbance events. I have these two Poisson models that I am picking from:

Model 1:

Count~Treatment+(1|StreamId)


Where Treatment has four levels: set1year1,set2year1,set1year2,set2year2

Model 2:

Count~StreamSet*Year+(1|StreamId)


Where StreamSet has two levels: Set1,Set2 and Year has two levels: Year1,Year2

Which to choose?

I recognize Model 2 as the traditional method for analyzing this type of experimental set up but a couple differences makes me lean toward Model 1:

• The resampling is 6 years apart
• This isn't technically a BACI design. We never recorded the "Before" measurement because the impact occurred before our first year of measurements. Instead it is more of a continuum of disturbance: i)unimpacted, ii)somewhat impacted, iii)more impacted, iv)impacted twice.

So my question is: Does using StreamId allow me to use Treatment without violating assumptions of independence?

The results you obtained for the two models should not be just similar, but exactly the same -- because both models have exactly the same fixed effects and the same random effects. The only difference is that the fixed effects portion is parameterized differently.

Here is an example using an available dataset to illustrate this. This is a split-plot experiment with a blocking structure on the whole plots, Variety as the whole-plot factor, and nitro as the split-plot factor (similar to repeated measures).

> data(Oats, package="nlme")
> Oats = transform(Oats, nv = interaction(nitro, Variety))


The factor nv represents the combinations of nitro and Variety. To make things more interesting, I'll throw out 12 observations at random:

> set.seed(1234)
> Oats = Oats[sample(1:72, 60), ]


I will fit equivalent models using nitro*Variety and nv as the fixed effects:

> library(lme4)
> modsep = lmer(yield ~ factor(nitro)*Variety + (1|Block/Variety), data = Oats)
> modint = lmer(yield ~ nv + (1|Block/Variety), data = Oats)


Now, look at the lsmeans results:

> lsmeans(modsep, ~ nitro*Variety)
nitro Variety        lsmean       SE    df  lower.CL  upper.CL
0.0 Golden Rain  81.67595 9.889542 24.82  61.30070 102.05121
0.2 Golden Rain  99.61801 9.216304 20.54  80.42536 118.81065
0.4 Golden Rain 117.59085 9.216626 20.53  98.39732 136.78437
0.6 Golden Rain 124.83333 8.742235 17.55 106.43258 143.23409
0.0 Marvellous   86.66667 8.742235 17.55  68.26591 105.06742
0.2 Marvellous  108.46275 9.905765 24.86  88.05572 128.86977
0.4 Marvellous  117.16667 8.742235 17.55  98.76591 135.56742
0.6 Marvellous  124.11076 9.244939 20.61 104.86242 143.35910
0.0 Victory      71.50000 8.742235 17.55  53.09924  89.90076
0.2 Victory      91.38669 9.250034 20.62  72.12846 110.64493
0.4 Victory     107.83807 9.937062 24.94  87.36987 128.30627
0.6 Victory     111.51815 9.920408 24.89  91.08223 131.95407

Confidence level used: 0.95

> lsmeans(modint, ~ nv)
nv                 lsmean       SE    df  lower.CL  upper.CL
0.Golden Rain    81.67595 9.889542 24.82  61.30070 102.05121
0.2.Golden Rain  99.61801 9.216304 20.54  80.42536 118.81065
0.4.Golden Rain 117.59085 9.216626 20.53  98.39732 136.78437
0.6.Golden Rain 124.83333 8.742235 17.55 106.43258 143.23409
0.Marvellous     86.66667 8.742235 17.55  68.26591 105.06742
0.2.Marvellous  108.46275 9.905765 24.86  88.05572 128.86977
0.4.Marvellous  117.16667 8.742235 17.55  98.76591 135.56742
0.6.Marvellous  124.11076 9.244939 20.61 104.86242 143.35910
0.Victory        71.50000 8.742235 17.55  53.09924  89.90076
0.2.Victory      91.38669 9.250033 20.62  72.12846 110.64493
0.4.Victory     107.83807 9.937062 24.94  87.36987 128.30627
0.6.Victory     111.51815 9.920408 24.89  91.08223 131.95407

Confidence level used: 0.95


The results are identical for the two fitted models -- estimates, SEs, df, and confidence limits.

I think some spurious error somewhere caused your two models to look different. I suggest that you go back to square one and re-fit both models from scratch, and I think you will find the same.

• Awesome, thanks for the repeatable code showing you're intuition was correct. I tracked my issue down to using glmmadmb, instead of lmer. I am not sure why the two glmmadmb models don't agree. I tried to run your example using glmmadmb but the models failed to converge... I am updating my first "answer" with this new information. – Adam C Jun 26 '15 at 18:59
• They failed to converge due to me using the wrong family. Here is the code to add to your example showing glmmadmb producing different model results: modsep2 = glmmadmb(yield ~ factor(nitro)*Variety , data = Oats, family="gaussian") modint2 = glmmadmb(yield ~ nv , data = Oats, family="gaussian") – Adam C Jun 26 '15 at 19:46

After the very helpful answer from @rvl my understanding is improved significantly, but the issue remains unexplained. It turns out the packages glmmADMB solves these two models slightly differently from each other, while the package lme4 gives the exact same answer. See below.

> mod1<-glmmadmb(count~Treatment+offset(area)+1|StreamId),
data=df,family="poisson")
> lsmeans(mod1,pairwise ~ Treatment, type="response")$lsmeans Treatment rate SE df asymp.LCL asymp.UCL Set1Year1 1.5157807 0.35555667 NA 0.95712798 2.4005055 Set2Year1 0.5412894 0.23822837 NA 0.22845719 1.2824906 Set1Year2 0.2816674 0.07497567 NA 0.16717037 0.4745847 Set2Year2 0.1107805 0.05169448 NA 0.04438739 0.2764821 Confidence level used: 0.95 > mod2<-glmmadmb(count~StreamSet*Year+offset(area)+1|StreamId), data=df,family="poisson") > lsmeans(mod2,pairwise ~ Impact*Year, type="response")$lsmeans
StreamSet Year      rate        SE   df  asymp.LCL asymp.UCL
Set1     Year1  1.5157779 0.63536250 NA 0.66656250 3.4469127
Set2     Year1  0.5412885 0.14797203 NA 0.31676400 0.9249574
Set1     Year2  0.2816687 0.14814817 NA 0.10046981 0.7896628
Set2     Year2  0.1107806 0.04002068 NA 0.05457054 0.2248896


lme4 models:

> mod1b<-glmer(count~Treatment+offset(area)+1|StreamId),
data=df,family=poisson)
> lsmeans(mod1b,pairwise ~ Treatment, type="response")$lsmeans Treatment rate SE df asymp.LCL asymp.UCL Set1Year1 1.4423756 0.35036760 NA 0.89601011 2.3219016 Set2Year1 0.5348784 0.14876125 NA 0.31011127 0.9225558 Set1Year2 0.3277244 0.08150271 NA 0.20128947 0.5335762 Set2Year2 0.1095906 0.03209355 NA 0.06173048 0.1945572 Confidence level used: 0.95 > mod2b<-glmer(count~StreamSet*Year+offset(area)+1|StreamId), data=df,family=poisson) > lsmeans(mod2b,pairwise ~ StreamSet*Year, type="response")$lsmeans
StreamSet  Year      rate         SE df asymp.LCL asymp.UCL
Set1     Year1 1.4423708 0.35036019 NA 0.8960148 2.3218742
Set2     Year1 0.5348812 0.14876006 NA 0.3101151 0.9225539
Set1     Year2 0.3277193 0.08150026 NA 0.2012877 0.5335642
Set2     Year2 0.1095903 0.03209306 NA 0.0617307 0.1945552


As you can see, Lme4 give identical results while glmmADMB gives different SE for the two models. Here is a reproducible example:

library(glmmADMB); library(lsmeans); library(lme4)
Owls <- transform(Owls, FoodSex = interaction(FoodTreatment, SexParent),
NegCount=as.integer(NegPerChick*BroodSize))

mod1<-glmer(NegCount~FoodTreatment*SexParent+(1|Nest),data=Owls,family=poisson)
mod2<-glmer(NegCount~FoodSex+(1|Nest),data=Owls,family=poisson)
lsmeans(mod1,~FoodTreatment*SexParent)
lsmeans(mod2,~FoodSex)


• Is Impact in this answer identical to StreamSet in the question? I guess I'm dumb, but I see no reason why there should be any difference at all between these two models. All you've done is parametrized the 3df for fixed effects in two different ways. So I am suspicious of your results -- and all the more so since the variable name changed. – Russ Lenth Jun 24 '15 at 22:26
• I think you edited the output to make the variable names match-- that is not the same as re-doing the analysis. Also, use the table function to confirm that your Treatment factor is defined correctly. I just posted an "answer" to this question with an example that illustrates a similar situation. – Russ Lenth Jun 26 '15 at 17:52
• I don't understand the discrepancy with glmmsdmb either. I considered perhaps that there is a bug in lsmeans support for those objects, but I did confirm with your examples that I obtain the same estimates and SEs as lsmeans does via the predict function with corresponding newdata. – Russ Lenth Jun 26 '15 at 21:45