For my experiment, I clipped plants and measured their responses, such as leaf mass produced, at the end of the season. I manipulated both clipping intensity and clipping time and crossed these two treatments. I also included a control clipped treatment resulting in 5 different clipping treatment combinations. With 12 plants per treatment there is a total of 60 plants which I followed over the course of two years. That is, I collected measurements on these 60 plants in year 1 and the same plants again in year 2.

Here is my design, where "never" under timing and "zero" under intensity arbitrarily replaced "control" treatment:

```
 Year   Timing  intensity   treatments
 2015   early   high       early-high
 2015   early   low        early-low
 2015   late    high       late-high
 2015   late    low        late-low
 2015   never   zero       control
 2014   early   high       early-high
 2014   early   low        early-low
 2014   late    high       late-high
 2014   late    low        late-low
 2014   never   zero       control

```

I followed one suggestion by Ben Bolker to ignore the warnings from running lme4 and afterwards run an F-test on the model ([R- analyzing repeated measures unbalanced design with lme4?](https://stackoverflow.com/questions/39135939/r-analyzing-repeated-measures-unbalanced-design-with-lme4)):

```
m1<-lmer(log(plant.leaf.g)~timing*intensity*year+(1|id), data=cmv)

Anova(m1, type="III", test="F")

```

The anova output gave me significant interactions between timing and intensity (p=0.006), and I followed up with a multiple comparison test using:

```
cmv$SHD<-interaction(cmv$timing, cmv$intensity)
m2<-lmer(log(plant.leaf.g)~-1+SHD+(1|id),data=cmv,  na.action=na.exclude)
summary(glht(m2, linfct=mcp(SHD="Tukey")))

```

Here is a clip of my output with the only significant pair where p=0.08:

```
                                Estimate Std. Error z value Pr(>|z|)  
late.2014 - early.2014 == 0   -0.6584     0.3448  -1.910   0.3844  
never.2014 - early.2014 == 0   0.1450     0.4102   0.354   0.9992  
early.2015 - early.2014 == 0  -0.4906     0.2786  -1.761   0.4788  
late.2015 - early.2014 == 0   -0.1687     0.3494  -0.483   0.9965  
never.2015 - early.2014 == 0   0.4201     0.4079   1.030   0.9032  
never.2014 - late.2014 == 0    0.8034     0.4119   1.951   0.3597  
early.2015 - late.2014 == 0    0.1678     0.3419   0.491   0.9963  
late.2015 - late.2014 == 0     0.4897     0.2724   1.797   0.4553  
never.2015 - late.2014 == 0    1.0785     0.4119   2.618   0.0885 .
early.2015 - never.2014 == 0  -0.6356     0.4074  -1.560   0.6133 

```

Why did Anova deem timing\*intensity to be highly significant, but no significance shows up in my multiple comparison test? Is there another way I should be doing the multiple comparison?

In other multiple comparison outputs I get p-values as high as 1.00000, is this normal?

```
data<-structure(list(id = c(91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 
99L, 100L, 101L, 102L, 103L, 105L, 106L, 107L, 108L, 109L, 110L, 
111L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L, 122L, 
123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 133L, 
134L, 135L, 136L, 137L, 138L, 139L, 140L, 141L, 142L, 143L, 144L, 
146L, 147L, 148L, 149L, 150L, 91L, 92L, 93L, 94L, 95L, 96L, 97L, 
98L, 99L, 100L, 101L, 102L, 103L, 105L, 106L, 107L, 108L, 109L, 
110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L, 
122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 
133L, 134L, 135L, 136L, 137L, 138L, 139L, 140L, 141L, 142L, 143L, 
144L, 146L, 147L, 148L, 149L, 150L), quad = c(2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), year = c(2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L), timing = structure(c(1L, 
3L, 2L, 1L, 1L, 2L, 3L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 3L, 1L, 1L, 
3L, 2L, 3L, 1L, 3L, 2L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 
1L, 3L, 3L, 2L, 2L, 1L, 2L, 3L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 2L, 
2L, 2L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 3L, 2L, 1L, 1L, 2L, 3L, 2L, 
2L, 2L, 2L, 1L, 1L, 2L, 3L, 1L, 1L, 3L, 2L, 3L, 1L, 3L, 2L, 3L, 
1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 2L, 
3L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 3L, 
1L), .Label = c("early", "late", "never"), class = "factor"), 
intensity = structure(c(2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 
1L, 2L, 1L, 1L, 2L, 3L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 
1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 3L, 3L, 2L, 1L, 1L, 
1L, 3L, 1L, 1L, 2L, 2L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 
2L, 1L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L, 
1L, 2L, 3L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 2L, 1L, 
1L, 2L, 2L, 2L, 2L, 1L, 2L, 3L, 3L, 2L, 1L, 1L, 1L, 3L, 1L, 
1L, 2L, 2L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 3L, 2L
), .Label = c("high", "low", "zero"), class = "factor"), 
treatment = structure(c(3L, 1L, 4L, 3L, 2L, 5L, 1L, 4L, 5L, 
4L, 5L, 2L, 2L, 5L, 1L, 3L, 2L, 1L, 4L, 1L, 2L, 1L, 4L, 1L, 
2L, 3L, 2L, 4L, 3L, 5L, 5L, 3L, 2L, 3L, 1L, 1L, 5L, 4L, 2L, 
4L, 1L, 4L, 2L, 3L, 5L, 4L, 1L, 3L, 4L, 5L, 4L, 2L, 3L, 5L, 
3L, 2L, 1L, 3L, 1L, 4L, 3L, 2L, 5L, 1L, 4L, 5L, 4L, 5L, 2L, 
2L, 5L, 1L, 3L, 2L, 1L, 4L, 1L, 2L, 1L, 4L, 1L, 2L, 3L, 2L, 
4L, 3L, 5L, 5L, 3L, 2L, 3L, 1L, 1L, 5L, 4L, 2L, 4L, 1L, 4L, 
2L, 3L, 5L, 4L, 1L, 3L, 4L, 5L, 4L, 4L, 3L, 5L, 2L, 1L, 3L
), .Label = c("control", "early-high", "early-low", "late-high", 
"late-low"), class = "factor"), plant.leaf.g = c(846.216, 
382.704, 2393.088, 61.832, 1315.86, 275.816, 3705.862, 3500.52, 
67.482, 432, 487.492, 1228.618, 776.16, 1575, 735.9, 2417.75, 
1342.92, 2359.046, 686.726, 1385.856, 343.684, 2277.312, 
465.528, 2314.584, 508.4, 1243.644, 1064.448, 1020.646, NA, 
494.832, 1318.248, 1516.4, 1271.218, 512.512, 157.878, 3753.992, 
586.032, 1042.176, 889.632, 651.052, 498.042, 625.872, 16.28, 
497.51, 593.75, 706.84, 2238.742, 232.584, 671.532, 90.72, 
1412.442, 902.728, 3077.184, 619.106, 0.576, 400.452, 684.522, 
849.852, 152.76, 1280.448, 274.47, 387.614, 98.496, 2304.504, 
644.952, 35.392, 250.56, 267.33, 2212.08, 2392.596, 751.944, 
629.418, 731.544, 1013.196, 1516.4, 130.536, 2910.6, 554.4, 
2163.35, 223.86, 2369.376, 551.976, 985.6, 1482.24, 815.386, 
1664.132, 596.376, 1581.432, 217.128, 1041.656, 951.168, 
256.172, 1587.148, 359.448, 546.48, 1226.544, 371.64, 293.504, 
177.726, 343.26, 691.24, 207.604, 588.924, 1405.258, 136.17, 
451.432, 576.18, 424.804, 884.534, 2466.45, 1524.432, 973.208, 
369.474, 410.048)), .Names = c("id", "quad", "year", "timing", 
"intensity", "treatment", "plant.leaf.g"), class = "data.frame", row.names = c(NA, 
-114L))

```

PS. I could not for the life of me get lsmeans to work with this unbalanced design. A lot of NAs are reported in the output.