Extract ANOVA treat*time interactions into a dataframe

A couple of points.

  • It is better not to use plist as the index, but rather the output of a call to seq_len() with length(plist).
  • We can store the p-values in a matrix, which we generically call out. We then assign column names to out, so that it is easier to appreciate which p-value belongs to which fixed effect.
  • We observe that one of the models has a hard time to estimate the variance of the random effects, as it returns boundary (singular) fit: see ?isSingular. This requires your attention if you encouter this with your own (non-simulated) data. Refer to this page for more information.
## snip ##
plist<- list("st","qt", "zt")
X <- model.matrix(~ treat * time, data = dat)
out <- matrix(rep(0L, length(plist) * dim(X)[2L]), ncol = 4L)
for (i in seq_len(length(plist))){
  model <- lmer(paste(plist[i], "~ (treat*time)", "+ (1|id)"), data=dat) 
  anovamodel <- (Anova(model, type=3))
  out[i, ] <- anovamodel$`Pr(>Chisq)`
}
colnames(out) <- colnames(model.matrix(model))
# --------------------------------------------------
> out
     (Intercept)  treatTrt    timeT2 treatTrt:timeT2
[1,]           0 0.3149593 0.7015615       0.3278066
[2,]           0 0.7774849 0.3511975       0.9013959
[3,]           0 0.5941231 0.1599605       0.9484378

We can save out as .csv file for later use in Excel.

# specify the folder where the file is to be stored
projdir <- 'my_directory'
write.csv(out, file = file.path(projdir, 'my_pvals.csv'))

# -----------------------------------------------------------------------
## write.csv2(out, file = file.path(projdir, 'my_pvals.csv'))
## to use a comma for the decimal point and a semicolon for the separator
## the Excel convention for CSV files in some Western European locales