[R-sig-ME] significance of slope (different than zero) in triple interaction

[Lenth, Russell V]

I’m confused. Do you want to test the slopes against zero, or against 1? I ask because you note that most of the confidence intervals include 1. To test against 1, add `null = 1.0` to the summary() call. To perform an equivalence test that the slopes do not differ from 1 by more than a specified threshold, also add `delta = (desired threshold)` to the call. See help(“summary.emmGrid”) for details.

Russ

From: Guillaume Adeux <guillaumesimon.a2 using gmail.com> 
Sent: Wednesday, September 5, 2018 10:13 AM
To: Lenth, Russell V <russell-lenth using uiowa.edu>
Cc: R-mixed models mailing list <r-sig-mixed-models using r-project.org>
Subject: Re: [R-sig-ME] significance of slope (different than zero) in triple interaction

Hi Russell,
Thank you very much for your answer. It's very nice of you.
I ran what you recommended and it produced the following output:

> em=emtrends(mod1,~syst|timing,var="year")
> summary(em,infer=c(TRUE,TRUE))
timing = après2:
 syst year.trend         SE  df  asymp.LCL  asymp.UCL z.ratio p.value
 S1   -0.9946823 0.05952514 Inf -1.1113495 -0.8780152 -16.710  <.0001
 S2   -0.8997828 0.07019881 Inf -1.0373699 -0.7621956 -12.818  <.0001
 S3   -0.8885606 0.06511346 Inf -1.0161806 -0.7609405 -13.646  <.0001
 S4   -0.9976324 0.06288404 Inf -1.1208828 -0.8743819 -15.865  <.0001
 S5   -0.7435081 0.06728649 Inf -0.8753872 -0.6116291 -11.050  <.0001

timing = avant1:
 syst year.trend         SE  df  asymp.LCL  asymp.UCL z.ratio p.value
 S1   -0.9758510 0.05696640 Inf -1.0875031 -0.8641990 -17.130  <.0001
 S2   -0.9170666 0.06973012 Inf -1.0537351 -0.7803980 -13.152  <.0001
 S3   -0.8793391 0.06482135 Inf -1.0063867 -0.7522916 -13.566  <.0001
 S4   -1.0951932 0.06253345 Inf -1.2177565 -0.9726298 -17.514  <.0001
 S5   -0.8061455 0.06589310 Inf -0.9352936 -0.6769974 -12.234  <.0001

However, I find this highly surprising because I know the slopes are not different than 0 for at least 4 of my systems at level "après2" (which means after). We can actually see that the confidence interval for all levels of syst at "timing=après2" embrace 1.
Do you have any explanation to this?
Thanks again for your interest.
Guillaume ADEUX

=========================
Le mer. 5 sept. 2018 à 16:25, Lenth, Russell V <mailto:russell-lenth using uiowa.edu> a écrit :
I'm a little confused because you refer to predictors by different names in different places, and 'year' seems to be used as both a covariate and a grouping factor. But try something like this:

    library(emmeans)
    emt <- emtrends(mod, ~ syst:timing, var = "year")
    summary(emt, infer = c(TRUE, TRUE))

This will estimate the slope for year at each combination of the two factors syst and timing. (You may need to re-fit the model after creating an additional variable, say WEED$syear <- scale(WEED$year), and with syear in place of scale(year) in the model formula and the emtrends call. You may follow-up with call(s) to emmeans::contrast(emt, ...) to compare or contrast these slopes.

Hope that helps.
-- Russ

Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science   
The University of Iowa  -  Iowa City, IA 52242  USA   
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017



-----Original Message-----
Date: Wed, 5 Sep 2018 10:43:39 +0200
From: Guillaume Adeux <mailto:guillaumesimon.a2 using gmail.com>
To: R-mixed models mailing list <mailto:r-sig-mixed-models using r-project.org>
Subject: [R-sig-ME] significance of slope (different than zero) in
        triple interaction


Hi mixmoders,

I have the following model:

mod=glmer(Weed_density~block+scale(year)*syst*timing+(1|year)+(1|plot)+(1|plot:year)+(1|ID_quadrat)+(1|OLRE)+offset(log(size_quadrat)),family=poisson(link="log"),dat=WEED)

I have a significant triple interaction between time : treatment : season.

Time is continuous, syst(=treatment) has 5 levels and season(=sampling
session) has two levels.

Here is the model output:

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: WDall ~ block + scale(year) * syst * timing + (1 | year) + (1
|      plot) + (1 | plot:year) + (1 | ID_quadrat) + (1 | OLRE) +
offset(log(size_quadrat))
   Data: WEED_paired_2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))

     AIC      BIC   logLik deviance df.resid
 21206.3  21371.9 -10577.2  21154.3     4286

Scaled residuals:
    Min      1Q  Median      3Q     Max
-1.6531 -0.4373 -0.1646  0.1426  2.6313

Random effects:
 Groups     Name        Variance  Std.Dev.
 OLRE       (Intercept) 4.456e-01 6.675e-01
 ID_quadrat (Intercept) 1.011e+00 1.006e+00  plot:year  (Intercept) 1.429e+00 1.195e+00
 year       (Intercept) 5.635e-15 7.506e-08
 plot       (Intercept) 0.000e+00 0.000e+00
Number of obs: 4312, groups:  OLRE, 4312; ID_quadrat, 2156; plot:year, 86; year, 17; plot, 10

Fixed effects:
                                Estimate Std. Error z value Pr(>|z|)
(Intercept)                     -0.84765    0.33352  -2.542 0.011036 *
blockD                          -0.28663    0.27596  -1.039 0.298971
scale(year)                      0.11385    0.25128   0.453 0.650500
systS2                           2.21797    0.43765   5.068 4.02e-07 ***
systS3                           2.97934    0.42857   6.952 3.61e-12 ***
systS4                           2.64787    0.43488   6.089 1.14e-09 ***
systS5                           0.55059    0.45565   1.208 0.226912
timingavant1                     1.87971    0.10286  18.275  < 2e-16 ***
scale(year):systS2               0.40061    0.38882   1.030 0.302863
scale(year):systS3               0.44798    0.37297   1.201 0.229698
scale(year):systS4              -0.01245    0.36549  -0.034 0.972819
scale(year):systS5               1.06031    0.37957   2.793 0.005215 **
scale(year):timingavant1         0.07949    0.09954   0.799 0.424489
systS2:timingavant1             -0.36039    0.12128  -2.972 0.002963 **
systS3:timingavant1             -0.56704    0.11777  -4.815 1.47e-06 ***
systS4:timingavant1             -0.39785    0.11984  -3.320 0.000901 ***
systS5:timingavant1             -0.06724    0.14990  -0.449 0.653770
scale(year):systS2:timingavant1 -0.15246    0.11992  -1.271 0.203628
scale(year):systS3:timingavant1 -0.04057    0.11556  -0.351 0.725543
scale(year):systS4:timingavant1 -0.49134    0.11614  -4.231 2.33e-05 ***
scale(year):systS5:timingavant1 -0.34391    0.13427  -2.561 0.010429 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


I wish to set up constrats to test if the slopes for scale(year):syst differ from zero at level 1 of timing.

It seems like we can do this with testInteractions but I'm not sure if my set up is correct:

testInteractions(mod1,custom=list(syst=c(1,0,0,0,0),timing=c(1,0)),
slope="scale(year)", adjustment="none")

The preceding code yields the following:

Adjusted slope for scale(year)
Chisq Test:
P-value adjustment method: none
                   Value Df  Chisq Pr(>Chisq)
syst1 : timing1 -0.82831  1 0.6464     0.4214

This doesn't seem correct because Value doesn't represent the slope for the first level of "syst" at the first level of "timing".

Could anyone shed their light?

Thank you very much!

Guillaume ADEUX