I am fitting a linear mixed effect model to study the interaction of two independent variables, a covariate time and a factor m (levels "R" and "V"), on the outcome variable var. Data are grouped by the variable id, and I consider the intercept for each group as a random variable. In particular, I am interested in whether var assumes different values depending on the levels of m at different time points.
Here is how I am fitting the model:
> options( contrasts = c( factor = "contr.treatment",ordered = "contr.poly")) > var.lm <- lme(var ~ m*time, random = ~1|id, data=dat, na.action=na.omit, method = "ML") > summary(var.lm) Linear mixed-effects model fit by maximum likelihood Data: dat AIC BIC logLik 2091.779 2117.895 -1039.89 Random effects: Formula: ~1 | id (Intercept) Residual StdDev: 4.534834e-05 1.480997 Fixed effects: var ~ m * time Value Std.Error DF t-value p-value (Intercept) 2.5819962 0.10587243 570 24.387806 0.0000 mV -0.8477854 0.14972621 570 -5.662238 0.0000 time 0.0103348 0.00490927 570 2.105150 0.0357 mV:time -0.0161672 0.00694276 570 -2.328636 0.0202 Correlation: (Intr) mV time mV -0.707 time -0.560 0.396 mV:time 0.396 -0.560 -0.707 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.6617004 -0.6537707 -0.2195518 0.4517715 3.9370134 Number of Observations: 574 Number of Groups: 1
So there seems to be a significant interaction between m and time!
The data are normalized with respect to their mean at time=-1 for each level of m, such that the expected value $$E(var|time=-1,m=R) = E(var|time=-1,m=V) = 1$$
This is confirmed as follows:
> mean(dat[dat$time==-1 & dat$m=="R",]$var, na.rm=TRUE)  1 > mean(dat[dat$time==-1 & dat$m=="V",]$var, na.rm=TRUE)  1 > boxplot(var~time*m, data=dat[dat$time %in% c(-1),])
I therefore expect, by construction, that the post-hoc analysis to test the null hypothesis E(var|time=-1,m=R) - E(var|time=-1,m=V) = 0 is not statistical significant. However, it turns out to be, and I do not understand why. Most probably I am doing something wrong.
Since I am using dummy coding for the factor m with the level "R" as the reference value (meaning that mV=1 if m==V), the null hypothesis should translate into the formula mV-mV:time=0. Therefore, I perform this tests as follows:
> ph_conditional <- c("mV - mV:time = 0"); > var.ph <- glht(var.lm, linfct = ph_conditional); > summary(var.ph) Simultaneous Tests for General Linear Hypotheses Fit: lme.formula(fixed = var ~ m * time, data = dat, random = ~1 | id, method = "ML", na.action = na.omit) Linear Hypotheses: Estimate Std. Error z value Pr(>|z|) mV - mV:time == 0 -0.8316 0.1532 -5.429 5.67e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method)
The null-hypothesis is rejected. There must be a problem and I cannot find it.
To make sure that my expectations are correct, I reduced the dataset to only the data at time=-1, and fit a model with no interactions. I expect the factor mV to be non statistically significant, meaning that the values of var at time=-1 and m=R are not statistically significantly different than those at time=-1 and m=V.
> dat_small <- dat[dat$time==-1,] > options( contrasts = c( factor = "contr.treatment",ordered = "contr.poly")) > var_small.lm <- lme(var ~ m, random = ~1|id, data=dat_small, na.action=na.omit, method = "ML") > summary(var_small.lm) Linear mixed-effects model fit by maximum likelihood Data: dat_small AIC BIC logLik 306.861 319.2114 -149.4305 Random effects: Formula: ~1 | id (Intercept) Residual StdDev: 1.242641e-05 0.6086403 Fixed effects: var ~ m Value Std.Error DF t-value p-value (Intercept) 1 0.06804805 160 14.6955 0 mV 0 0.09623447 160 0.0000 1 Correlation: (Intr) mV -0.707 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.4611461 -0.7789445 -0.2552357 0.6026049 4.6896472 Number of Observations: 162 Number of Groups: 1
As expected, mV is not significant. Therefore, I am doing some mistakes during the post-hoc analysis of the model with interaction. Any help is very much appreciated.
You can find the dataset as well as the code in this link. This is the dataset from which I have computed the results in this post, but it is not the complete dataset of my project. That's why I am using a random effect on the variable id even if in this reduced dataset there is only one group. I get very similar results by using the complete dataset. Thanks a lot!