为什么我的混合模型循环不能工作？(RStudio，交叉设计)

Question

我搞不懂为什么我的循环不起作用。

我有一个数据库( 36行x51列，其名称为"Seleccio")，由3个因素组成(前3栏:动物(12只动物)、饮食(3种饮食)和周期(3期))和48个变量(许多临床参数)，每栏36条。这是一个3x3交叉设计，所以我想实现一个混合模型，包括动物随机效应，也包括周期和饮食固定效应，以及它们之间的相互作用。

数据示例(但行和列较少)：

  Animal Diet  Period  Var1  Var2  Var3  Var4  Var5  Var6
  <chr>  <chr> <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A      A     A         11    55   1.2 0.023    22     3
2 B      A     A         13    34   1.6 0.04     23     4
3 C      B     A         15    13   1.1 0.052    22     2
4 A      B     B         10    22   1.5 0.067    27     4
5 B      C     B          9    45   1.4 0.012    24     2
6 C      C     B         13    32   1.5 0.014    23     3

> dput(sample[1:9,])
structure(list(Animal = c("A", "B", "C", "A", "B", "C", NA, NA, 
NA), Diet = c("A", "A", "B", "B", "C", "C", NA, NA, NA), Period = c("A", 
"A", "A", "B", "B", "B", NA, NA, NA), Var1 = c(11, 13, 15, 10, 
9, 13, NA, NA, NA), Var2 = c(55, 34, 13, 22, 45, 32, NA, NA, 
NA), Var3 = c(1.2, 1.6, 1.1, 1.5, 1.4, 1.5, NA, NA, NA), Var4 = c(0.023, 
0.04, 0.052, 0.067, 0.012, 0.014, NA, NA, NA), Var5 = c(22, 23, 
22, 27, 24, 23, NA, NA, NA), Var6 = c(3, 4, 2, 4, 2, 3, NA, NA, 
NA)), row.names = c(NA, -9L), class = c("tbl_df", "tbl", "data.frame"
))

我想对每个变量进行描述性分析(正态性检验和异常值检验)，并对固定效果进行方差分析( ANOVA )和图基检验。

我可以一个一个地进行分析，但这需要很长时间，我已经多次尝试创建一个for循环来自动化对所有变量的分析，但是它没有工作(我对R非常陌生)。

到目前为止我得到的是：

sink("output.txt") # to store the output into a file, as it would be to large to be shown in the console
vars <-as.data.frame(Seleccio[,c(4:51)])
fact <-Seleccio$Diet
dim(vars)
for (i in 1:length(vars)) { 
  variable <- vars[,i]
  lme_cer <- lme(variable ~ Period*Diet, random = ~1 | Animal, data = Seleccio) # the model
  cat("\n---------\n\n")
  cat(colnames(Seleccio)[i]) # the name of each variable, so I don't get lost in the text file
  cat("\n")
  print(boxplot(vars[,i]~fact)$out) #checking for outliers
  print(summary(lme_cer))
  print(anova(lme_cer)) 
  print(lsmeans(lme_cer, pairwise~Diet, adjust="tukey"))
}
sink()

这段代码运行但没有完成这项工作，因为它给出了每个变量的错误结果，因为它们与我逐个分析每个变量时得到的结果不同。我还想在循环中添加这个按Diet (治疗)代码排序的正常测试。我想知道这是否可能。

aggregate(formula = VARIABLENAME ~ Diet,
          data = Seleccio,
          FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})

非常感谢所有愿意帮助我的人，任何帮助都将不胜感激。

Answer 1

我不认为我能运行的模型只有6个观察，所以我找不出为什么你的循环将不返回相同的做法，一个一个。也许cat(colnames(Seleccio)[i])的问题在于:您只需要Var名称，对于i=1，2和3，该代码将返回“动物”、“饮食”和“周期”，从而混淆了比较结果的方式。使用cat(colnames(vars)[i])可能会纠正这种情况。如果您找到了一种包含更多Seleccio观察的方法，我可能会提供更多帮助。

我建议您创建一个列表来存储输出：

vars <- as.data.frame(Seleccio[,c(4:51)])
fact <- Seleccio$Diet
dim(vars)
output = list() #Create empty list
for (i in 1:length(vars)) {
  var = colnames(vars)[i] 
  output[[var]] = list() #Create one entry for each variable
  variable <- vars[,i]
  lme_cer <- lme(variable ~ Period*Diet, random = ~1 | Animal, data = Seleccio) # the model
  #Fill that entry with each statistics:
  output[[var]]$boxplot = boxplot(vars[,i]~fact)$out #checking for outliers
  output[[var]]$summary = summary(lme_cer)
  output[[var]]$anova = anova(lme_cer)
  output[[var]]$lsmeans = lsmeans(lme_cer, pairwise~Diet, adjust="tukey")
  output[[var]]$shapiro = aggregate(formula = variable ~ Diet, data = Seleccio,
            FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})
}

这样，您就可以在R环境中获得结果，并有更好的可视化选项:执行输出$ Var1并获取Var1的所有结果，这些结果应该适合控制台；执行for(i in output){print(i$summary)}以获取所有摘要，等等。