# Mixed-model in R - Am I doing this right? Why is it always significant?

Here is my experiment: I evaluated the same individuals before applying treatment1, treatment2, a combination of both or none of them (pre_treatment_value); I then treated them as described before, and then evaluated them all 8 times after that (evaluation_number). I want to test whether or not there is an effect of any of of these treatments or of their interaction.

Can you confirm that I am using the right methodology to get some p-values to report:

library(lme4)

#simulating dataset
id <- rep(1:8,8)
evaluation_number <- rep(c(1:8), each=8)
treatment1 <- rep(c("t1","t1","t1","t1","control","control","control","control"),8)
treatment2 <- rep(c("t2","t2","control","control","control","control","t2","t2"),8)
dependant_variable <- rnorm(64,5,10)
pre_treatment_value <- rnorm(8,3,5)

evaluation_number

df <- as.data.frame(cbind(id,evaluation_number,treatment1,treatment2,dependant_variable,pre_treatment_value))
df$dependant_variable <- as.numeric(as.character(df$dependant_variable))
df$pre_treatment_value <- as.numeric(as.character(df$pre_treatment_value))

#analyses

mod1 <- lmer(dependant_variable ~
treatment1
+treatment2
+treatment1:evaluation_number
+treatment2:evaluation_number
+evaluation_number:treatment1:treatment2
+pre_treatment_value
+ (1|id),
data=df, REML=T)

mod2 <- lmer(dependant_variable ~
treatment1
+treatment2
+treatment1:evaluation_number
+treatment2:evaluation_number
+evaluation_number:treatment1:treatment2
+ (1|id),
data=df, REML=T)

mod3 <- lmer(dependant_variable ~
treatment1
+treatment2
+treatment1:evaluation_number
+treatment2:evaluation_number
+ (1|id),
data=df, REML=T)

mod4 <- lmer(dependant_variable ~
treatment1
+treatment2
+treatment1:evaluation_number
+ (1|id),
data=df, REML=T) ###and so on

anova(mod1,mod2,refit=F,test="Chisq")
anova(mod2,mod3,refit=F,test="Chisq")
anova(mod3,mod4,refit=F,test="Chisq") ###and so on


I noticed while simulating these results randomly that in almost every simulation I would get significant results. Either I am doing this wrong or this analysis is not conservative enough.

• You have a very common, but very bad bug in your code. as.numeric for factors does not do what most people think it does. For factors created from numeric values, as.numeric does not return the numeric values, but an internal representation. Running your example I get 12 for the first value, which becomes 31 after as.numeric (31st value after ordering numbers as strings). See cran.r-project.org/doc/FAQ/… and stackoverflow.com/questions/3418128/… – LiKao Jan 30 '18 at 10:06
• Thank you very much. I just checked my code and I only used as.numeric to create the repeatable dataset, my actual values were correctly included in the real model. – Nakx Jan 30 '18 at 10:20
• @LiKao actually I changed my code, also realized I should probably use data.frame(). – Nakx Jan 31 '18 at 0:00

The design of this analysis was hard to interpret and I finally considered 4 groups instead of three with an interaction.

My final analysis was a more classic ANCOVA:

mod1 <- lmer(dependant_variable ~
treatment*evaluation_number
+pre_treatment_value
+ (1|id),
data=df)


with "treatment" being a variable with four levels t1, t2, t1t2, control.

• I don't think that this code will run... And what do you mean by "a more classic ANCOVA"? – Stefan Jan 30 '18 at 16:56
• Instead of having the weird separate treatment design (t1 and t2 in two separate vectors that can interact together for some individuals) I consider it as another treatment t1t2. This issue was mostly that I couldn't understand whether the same levels would be used twice in the analysis or not by using two different vectors. With the new analysis, the results are coherent with what I observed and with my hypothesis and this analysis is much easier to interpret. I realized I wrote aov instead of lmer, is that why you think the code won't run? I edited the answer. – Nakx Jan 30 '18 at 23:51