r语言 - 为什么我的混合模型循环不起作用？(RStudio，分频器设计)

Question

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

我有一个数据库(36行×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", "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"
))

我想对按饮食排序的每个变量进行描述性分析(正态性检验和异常值检查((这是一种治疗方法(，并运行一个混合模型，对固定效应进行方差分析和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---------nn")
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)})

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

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)}以获取所有摘要；等