[R-sig-ME] significance of slope (different than zero) in triple interaction
Guillaume Adeux
guill@ume@imon@@2 @ending from gm@il@com
Thu Sep 6 12:03:23 CEST 2018
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 <guillaumesimon.a2 using 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
>
> Le mer. 5 sept. 2018 à 18:54, Ben Bolker <bbolker using gmail.com> a écrit :
>
>>
>> Even more confusing, your slope CI for après2 seem to (mostly) include
>> -1, not 1 (except for S5). Maybe use `null = -1.0` ?
>>
>> On 2018-09-05 12:38 PM, Lenth, Russell V wrote:
>> > 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
>> >
>> > _______________________________________________
>> > R-sig-mixed-models using r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list