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

Discussion:

Guillaume Adeux

2018-09-05 08:43:39 UTC

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

[[alternative HTML version deleted]]

Lenth, Russell V

2018-09-05 14:25:13 UTC

Permalink

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 <***@gmail.com>
To: R-mixed models mailing list <r-sig-mixed-***@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

Lenth, Russell V

2018-09-05 16:38:44 UTC

Permalink

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 <***@gmail.com>
Sent: Wednesday, September 5, 2018 10:13 AM
To: Lenth, Russell V <russell-***@uiowa.edu>
Cc: R-mixed models mailing list <r-sig-mixed-***@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.

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-***@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:***@gmail.com>
To: R-mixed models mailing list <mailto:r-sig-mixed-***@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

Ben Bolker

2018-09-05 16:50:58 UTC

Permalink

Even more confusing, your slope CI for après2 seem to (mostly) include
-1, not 1 (except for S5). Maybe use `null = -1.0` ?

Post by 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
Sent: Wednesday, September 5, 2018 10:13 AM
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.

em=emtrends(mod1,~syst|timing,var="year")
summary(em,infer=c(TRUE,TRUE))

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
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
=========================
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
Subject: [R-sig-ME] significance of slope (different than zero) in
triple interaction
Hi mixmoders,
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.
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
Min 1Q Median 3Q Max
-1.6531 -0.4373 -0.1646 0.1426 2.6313
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
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.
testInteractions(mod1,custom=list(syst=c(1,0,0,0,0),timing=c(1,0)),
slope="scale(year)", adjustment="none")
Adjusted slope for scale(year)
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
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Guillaume Adeux

2018-09-05 17:08:23 UTC

Permalink

Sorry for the trouble.

What I want to test is that the slopes are effectively increasing or
decreasing; that back on the original scale, the regression line is not
parralel to the x-axis. This must be possible but I don't even know if this
is the way to go.

Thank you very much for you interest,

Guillaume ADEUX

Post by Ben Bolker
Even more confusing, your slope CI for après2 seem to (mostly) include
-1, not 1 (except for S5). Maybe use `null = -1.0` ?

Post by 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.

Post by Lenth, Russell V
Russ
Sent: Wednesday, September 5, 2018 10:13 AM
Subject: Re: [R-sig-ME] significance of slope (different than zero) in

triple interaction

Post by Lenth, Russell V
Hi Russell,
Thank you very much for your answer. It's very nice of you.

em=emtrends(mod1,~syst|timing,var="year")
summary(em,infer=c(TRUE,TRUE))

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.

Post by Lenth, Russell V
Do you have any explanation to this?
Thanks again for your interest.
Guillaume ADEUX
=========================
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

Post by Lenth, Russell V
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.

Post by Lenth, Russell V
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
Subject: [R-sig-ME] significance of slope (different than zero) in
triple interaction
Hi mixmoders,

Post by Lenth, Russell V
Time is continuous, syst(=treatment) has 5 levels and season(=sampling
session) has two levels.
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))

Post by Lenth, Russell V
AIC BIC logLik deviance df.resid
21206.3 21371.9 -10577.2 21154.3 4286
Min 1Q Median 3Q Max
-1.6531 -0.4373 -0.1646 0.1426 2.6313
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

Post by Lenth, Russell V
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

Post by Lenth, Russell V
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.

Post by Lenth, Russell V
It seems like we can do this with testInteractions but I'm not sure if
testInteractions(mod1,custom=list(syst=c(1,0,0,0,0),timing=c(1,0)),
slope="scale(year)", adjustment="none")
Adjusted slope for scale(year)
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".

Post by Lenth, Russell V
Could anyone shed their light?
Thank you very much!
Guillaume ADEUX
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

Guillaume Adeux

2018-09-06 10:03:23 UTC

Permalink

Ok i found out that was messing up everything.

When an offset is included in the model, the slopes produced by
emtrends(model,~factor1|factor2,"num_var") are very different than the ones
that can be computed with the model output (the coefficients for the slopes
in the output need to be divided by the standard deviation of "num_var" in
case of a scaled(num_var) in order to be directly comparable with emtrends
output). I don't know what emtrends is doing in this case.

However, the slopes produced by emtrends(model,~factor1|factor2,"num_var")
are equivalent to the back transformed scaled coefficients of slopes of the
summary output when the offset is taken out. Is there any way to keep the
offset and have the correct values of slope? I do not know.

To come back to my original questions, I was indeed looking for the "null"
argument which I need to set to 0 because the tests are done on the scale
of the linear predictor and a slope of 0 on the log scale is equivalent to
a multiplicative factor of 1 on the original scale (which means a flat
regression line - no effect).

Don't hesitate to give me your feedback if you believe something I have
said above is incorrect.

Thank you a thousand time for you help.

Guillaume ADEUX

Post by Guillaume Adeux
Sorry for the trouble.
What I want to test is that the slopes are effectively increasing or
decreasing; that back on the original scale, the regression line is not
parralel to the x-axis. This must be possible but I don't even know if this
is the way to go.
Thank you very much for you interest,
Guillaume ADEUX

Post by Ben Bolker
Even more confusing, your slope CI for après2 seem to (mostly) include
-1, not 1 (except for S5). Maybe use `null = -1.0` ?

Post by 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.

Post by Lenth, Russell V
Russ
Sent: Wednesday, September 5, 2018 10:13 AM
Subject: Re: [R-sig-ME] significance of slope (different than zero) in

triple interaction

Post by Lenth, Russell V
Hi Russell,
Thank you very much for your answer. It's very nice of you.

em=emtrends(mod1,~syst|timing,var="year")
summary(em,infer=c(TRUE,TRUE))

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.

names in different places, and 'year' seems to be used as both a covariate

Post by Lenth, Russell V
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

2e+05))

1.429e+00 1.195e+00

Post by Lenth, Russell V
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

differ from zero at level 1 of timing.

the first level of "syst" at the first level of "timing".

Post by Lenth, Russell V
Could anyone shed their light?
Thank you very much!
Guillaume ADEUX
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

Lenth, Russell V

2018-09-06 15:05:24 UTC

Permalink

The emtrends function computes difference quotients using the 'var' argument, which is actually processed as an expression. That makes it possible to scale the slopes. For example (with a different model):

require(emmeans)
fiber.lm <- lm(strength ~ diameter*machine, data=fiber)

summary(emtrends(fiber.lm, "machine", var = "diameter"), infer = TRUE)
## machine diameter.trend SE df lower.CL upper.CL t.ratio p.value
## A 1.1042781 0.1936634 9 0.6661810 1.542375 5.702 0.0003
## B 0.8571429 0.2238228 9 0.3508205 1.363465 3.830 0.0040
## C 0.8641975 0.2080707 9 0.3935090 1.334886 4.153 0.0025
## Confidence level used: 0.95

summary(emtrends(fiber.lm, "machine", var = "2*diameter"), infer = TRUE)
## machine 2*diameter.trend SE df lower.CL upper.CL t.ratio p.value
## A 0.5521390 0.0968317 9 0.3330905 0.7711876 5.702 0.0003
## B 0.4285714 0.1119114 9 0.1754102 0.6817326 3.830 0.0040
## C 0.4320988 0.1040353 9 0.1967545 0.6674430 4.153 0.0025
## Confidence level used: 0.95## Confidence level used: 0.95

The second table has estimates and SEs half as large, because we are differentiating with respect to 2*diameter. Note however that scaling the slopes has no effect on the t (or z) ratios or on the P values. The tests given, in both cases, are tests against the slope being zero. In your own results earlier in this discussion, all the slopes are decidedly negative, no matter how you scale them.

Russ

PS to Ben: Good call -- I should have said -1.

======================
From: Guillaume Adeux <***@gmail.com>
Sent: Thursday, September 6, 2018 5:03 AM
To: Ben Bolker <***@gmail.com>; Lenth, Russell V <russell-***@uiowa.edu>
Cc: R-mixed models mailing list <r-sig-mixed-***@r-project.org>
Subject: Re: [R-sig-ME] significance of slope (different than zero) in triple interaction

Ok i found out that was messing up everything.
When an offset is included in the model, the slopes produced by emtrends(model,~factor1|factor2,"num_var") are very different than the ones that can be computed with the model output (the coefficients for the slopes in the output need to be divided by the standard deviation of "num_var" in case of a scaled(num_var) in order to be directly comparable with emtrends output). I don't know what emtrends is doing in this case.
However, the slopes produced by emtrends(model,~factor1|factor2,"num_var") are equivalent to the back transformed scaled coefficients of slopes of the summary output when the offset is taken out. Is there any way to keep the offset and have the correct values of slope? I do not know.
To come back to my original questions, I was indeed looking for the "null" argument which I need to set to 0 because the tests are done on the scale of the linear predictor and a slope of 0 on the log scale is equivalent to a multiplicative factor of 1 on the original scale (which means a flat regression line - no effect).
Don't hesitate to give me your feedback if you believe something I have said above is incorrect.
Thank you a thousand time for you help.
Guillaume ADEUX

Le mer. 5 sept. 2018 à 19:08, Guillaume Adeux <mailto:***@gmail.com> a écrit :
Sorry for the trouble.

What I want to test is that the slopes are effectively increasing or decreasing; that back on the original scale, the regression line is not parralel to the x-axis. This must be possible but I don't even know if this is the way to go.

Thank you very much for you interest,

Guillaume ADEUX