为什么我的混合模型循环不起作用？ RStudio，跨界设计

编程问答 2022-04-29

问题描述

我不明白为什么我的循环不起作用。

我有一个数据库（36 行 x 51 列，名称为“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",NA,NA),Diet = c("A",Period = c("A",Var1 = c(11,13,15,10,9,Var2 = c(55,34,22,45,32,Var3 = c(1.2,1.6,1.1,1.5,1.4,Var4 = c(0.023,0.04,0.052,0.067,0.012,0.014,Var5 = c(22,23,27,24,Var6 = c(3,4,2,3,NA)),row.names = c(NA,-9L),class = c("tbl_df","tbl","data.frame"
))

我想对按饮食（治疗）排序的每个变量进行描述性分析（正态性检验和异常值检查），并运行混合模型并对固定效应进行方差分析和 Tukey 检验。>

我可以一个一个的做分析，但是需要很多时间，我尝试了几次创建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()

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

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

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

解决方法

我不认为我可以只用 6 个观察来运行模型，所以我找不到为什么你的循环不会返回与一个一个一样的结果。也许问题出在 cat(colnames(Seleccio)[i]) 上：您只需要 Var 名称，而对于 i=1、2 和 3，该代码将返回“Animal”、“Diet”和“Period”，从而弄乱了您的方式重新比较结果。使用 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 环境中获得结果，并有更好的可视化选项：执行 output$Var1 并获取 Var1 的所有结果，它应该适合控制台；执行 for(i in output){print(i$summary)} 以获取所有摘要；等

crossoverfor-looploopsnlmerr