[R-sig-ME] significance of slope (different than zero) in triple interaction
Guillaume Adeux
guill@ume@imon@@2 @ending from gm@il@com
Wed Sep 5 19:08:23 CEST 2018
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