为什么我的混合模型循环不工作？ （RStudio，跨界设计）
Why isn't my mixed model loop working? (RStudio, Crossover design)
我不明白为什么我的循环不工作。
我有一个数据库（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", "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---------\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()
此代码 运行s 但没有完成工作，因为它为我提供了每个变量的错误结果，因为它们与我逐一分析每个变量时得到的结果不同。
我还想将这个按饮食（治疗）代码排序的常态性测试添加到循环中。不知是否可行
aggregate(formula = VARIABLENAME ~ Diet,
data = Seleccio,
FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})
非常感谢所有愿意帮助我的人，任何帮助将不胜感激
我不认为我可以 运行 只有 6 个观察值的模型，所以我找不到为什么你的循环不 return 与一个一个地做一样.也许问题出在 cat(colnames(Seleccio)[i]) 上：您只需要 Var 名称，对于 i=1、2 和 3，该代码将 return “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)} 获取所有摘要；等等
我不明白为什么我的循环不工作。
我有一个数据库（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", "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---------\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()
此代码 运行s 但没有完成工作，因为它为我提供了每个变量的错误结果，因为它们与我逐一分析每个变量时得到的结果不同。 我还想将这个按饮食（治疗）代码排序的常态性测试添加到循环中。不知是否可行
aggregate(formula = VARIABLENAME ~ Diet,
data = Seleccio,
FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})
非常感谢所有愿意帮助我的人，任何帮助将不胜感激
我不认为我可以 运行 只有 6 个观察值的模型，所以我找不到为什么你的循环不 return 与一个一个地做一样.也许问题出在 cat(colnames(Seleccio)[i]) 上：您只需要 Var 名称，对于 i=1、2 和 3，该代码将 return “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)} 获取所有摘要；等等