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

Ben Bolker bbolker @ending from gm@il@com
Wed Sep 5 18:50:58 CEST 2018 

 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
>

More information about the R-sig-mixed-models mailing list