将方差分析处理*时间交互提取到数据框中

Question

我正在分析一项随机对照试验，并正在研究治疗*时间相互作用测试。我希望将方差分析的 p 值提取到数据框中，以便导出到 excel。目前，我的代码使用 p 值作为维度 [1:4] 的数值向量对对象进行编程。但是，当我将其复制到 excel 时，数据被转录到每行一个单元格中，值由空格分隔，而不是每个 p 值占据自己的单元格。

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)
}

Answer 1

几点。

最好不要使用 plist 作为索引，而是使用 length(plist).

seq_len()

我们可以将p-values存储在一个矩阵中，我们一般称之为out。然后我们将列名分配给 out，这样更容易理解哪个 p-value 属于哪个固定效应。
我们观察到其中一个模型很难估计随机效应的方差，因为它 returns boundary (singular) fit: see ?isSingular。如果您使用自己的 (non-simulated) 数据遇到此问题，则需要您的注意。有关详细信息，请参阅 this page。

## 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

我们可以将 out 保存为 .csv 文件以备后用 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

将方差分析处理*时间交互提取到数据框中

Extract ANOVA treat*time interactions into a dataframe

r

export

anova