I am analyzing an RCT and am looking at treat*time interaction tests. I wish to extract the p-values for an ANOVA into a data frame for exporting into excel. Presently, my code programs an object with the p-values as a numerical vector of dimensions [1:4]. However, when I copy this into excel, the data is transcribed into one cell per line with the values separated by spaces rather than each p-value occupying its own cell.

library(tidyverse)
library(rstatix)
library(lme4)

set.seed(42)

n <- 1000
dat1 <- data.frame(id=1:n,
                  treat = factor(sample(c('Trt','Ctrl'), n, rep=TRUE, prob=c(.5, .5))),
                  time = factor("T1"),
                  outcome1=rbinom(n = 1000, size = 1, prob = 0.3),
                  st=runif(n, min=24, max=60),
                  qt=runif(n, min=.24, max=.60),
                  zt=runif(n, min=124, max=360)
                  )
dat2 <- data.frame(id=1:n,
                   treat = dat1$treat,
                   time = factor("T2"),
                   outcome1=dat1$outcome1,
                   st=runif(n, min=24, max=60),
                   qt=runif(n, min=.24, max=.60),
                   zt=runif(n, min=124, max=360)
)
dat <- rbind(dat1,dat2)

id <- dat$id
st <- dat$st
qt <- dat$qt
zt <- dat$zt
treat <- dat$treat
time <- dat$time
plist<- list("st","qt", "zt")
for (i in plist){ 
  model <- lmer(paste(i, "~ (treat*time)", "  (1|id)"), data=dat) 
  anovamodel <- (Anova(model, type=3))
  grpxtime <- anovamodel$`Pr(>Chisq)`
  print(grpxtime)
} 

CodePudding user response：

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