[SOLVED] Significant interaction in linear mixed effect model but plot shows overlapping confidence intervals?

I try to show you as much as possible of the structure of the data and the results I produced.

The structure of the data is the following:

   GroupID Person Factor2 Factor1 Rating
    <int>  <int>  <fctr>  <fctr>  <int>
1       2    109       2       0      1
2       2    109       2       1     -2
3       2    104       1       0      4
4       2    236       1       1      1
5       2    279       1       1      2
6       2    179       2       1      0

Person is the participant ID, GroupID is the kind of stimulus rated, Factor 1 (levels 0 and 1) and Factor 2 (levels 1 and 2) are fixed factors and the Ratings are the outcome variables.

I am trying to print a plot for a significant interaction in a linear mixed effect model. I used the packages lme4 and lmerTest to analyze the data.

This is the model we ran:

> model_interaction <- lmer(Rating ~ Factor1 * Factor2 + ( 1 | Person) +
(1 | GroupID), data)
> model_interaction

Linear mixed model fit by REML ['merModLmerTest']
Formula: Rating ~ Factor1 * Factor2 + (1 | Person) + (1 | GroupID)
Data: data
REML criterion at convergence: 207223.9

Random effects:
 Groups   Name        Std.Dev.
 Person   (Intercept) 1.036   
 GroupID  (Intercept) 1.786   
 Residual             1.880   
Number of obs: 50240, groups:  Person, 157; GroupID, 80
Fixed Effects:
  (Intercept)           Factor11           Factor22  Factor11:Factor22  
     -0.43823            0.01313            0.08568            0.12440

When I use the summary() function R returns the following output

> summary(model_interaction)

Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom
['lmerMod']
Formula: Rating ~ Factor1 * Factor2 + (1 | Person) + (1 | GroupID)
   Data: data

REML criterion at convergence: 207223.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.8476 -0.6546 -0.0213  0.6516  4.2284 

Random effects:
 Groups   Name        Variance Std.Dev.
 Person   (Intercept) 1.074    1.036   
 GroupID  (Intercept) 3.191    1.786   
 Residual             3.533    1.880   
Number of obs: 50240, groups:  Person, 157; GroupID, 80

Fixed effects:
                    Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)       -4.382e-01  2.185e-01  1.110e+02  -2.006 0.047336 *  
Factor11           1.313e-02  2.332e-02  5.004e+04   0.563 0.573419    
Factor22           8.568e-02  6.275e-02  9.793e+03   1.365 0.172138    
Factor11:Factor22  1.244e-01  3.385e-02  5.002e+04   3.675 0.000238 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Fctr11 Fctr22
Factor11    -0.047              
Factor22    -0.135  0.141       
Fctr11:Fc22  0.034 -0.694 -0.249

I know it is not possible to interpret p-Values for linear mixed effects model. So I ran an additional anova, comparing the interaction model to a model with just the main effects of Factor1 and Factor2

> model_Factor1_Factor2 = lmer(Rating ~ Factor1 + Factor2 +
  ( 1 | Person) + (1 | GroupID), data)
> anova(model_Factor1_Factor2, model_interaction)

Data: data
Models:
object: Rating ~ Factor1 + Factor2 + (1 | Person) + (1 | GroupID)
..1: Rating ~ Factor1 * Factor2 + (1 | Person) + (1 | GroupID)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
object  6 207233 207286 -103611   207221                             
..1     7 207222 207283 -103604   207208 13.502      1  0.0002384 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I interpreted this Output as: the Interaction of Factor1 and Factor2 explains additional variance in my outcome measurement compared to the model with only the main effects of Factor1 and Factor2.

Since interpreting output for linear mixed effects models is hard I would like to print a graph showing the interaction of Factor1 and Factor2. I did so using lsmeans package (first I used the plot(allEffects) but after reading this How to get coefficients and their confidence intervals in mixed effects models? question I realized that it was not the right way to print graphs for linear mixed effect models).

So this is what I did (following this Website http://rcompanion.org/handbook/G_06.html)

> leastsquare = lsmeans(model_interaction, pairwise ~ Factor2:Factor1,
 adjust="bon")
> CLD = cld(leastsquare, alpha=0.05, Letters=letters, adjust="bon")
> CLD$.group=gsub(" ", "", CLD$.group)
> CLD
 Factor2 Factor1     lsmean        SE     df   lower.CL  upper.CL .group
 1       0       -0.4382331 0.2185106 111.05 -0.9930408 0.1165746  a    
 1       1       -0.4251015 0.2186664 111.36 -0.9803048 0.1301018  a    
 2       0       -0.3525561 0.2190264 112.09 -0.9086735 0.2035612  a    
 2       1       -0.2150234 0.2189592 111.95 -0.7709700 0.3409233   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 4 estimates 
P value adjustment: bonferroni method for 6 tests 
significance level used: alpha = 0.05

This is the plotting funtion I used

> ggplot(CLD, aes(`Factor1`, y = lsmean, ymax = upper.CL,
  ymin = lower.CL, colour = `Factor2`, group = `Factor2`)) +
  geom_pointrange(stat = "identity",
  position = position_dodge(width = 0.1)) +
  geom_line(position = position_dodge(width = 0.1))

The plot can be found using this link (I am not allowed to post images yet, please excuse the workaround)

Interaction of Factor1 and Factor2

Now my question is the following: Why do I have a significant interaction and a significant amount of explained variance by this interaction but my confidence intervals in the plot overlap? I guess I did something wrong with the confidence intervals? Or is it because it is just not possible to interpret the significance indices for linear mixed effects models?