【发布时间】：2020-06-15 13:26:56
【问题描述】：
考虑这个数据框：
set.seed(123)
dat1 <- data.frame(Region = rep(c("r1","r2"), each = 100),
State = rep(c("NY","MA","FL","GA"), each = 10),
Loc = rep(c("a","b","c","d","e","f","g","h"),each = 5),
ID = rep(c(1:10), each = 2),
var1 = rnorm(200),
var2 = rnorm(200),
var3 = rnorm(200),
var4 = rnorm(200),
var5 = rnorm(200),
var6 = rnorm(200))
dat1$ID <- factor(dat1$ID)
我正在使用 mclust 包来拟合混合模型并获得自举标准错误，如下所示：
library(tidyverse)
library(mclust)
set.seed(123)
mod <- Mclust(dat1[,5:10], G=5, modelNames = "VEI", initialization = list(hcPairs = randomPairs(dat1[,5:10], seed = 123)))
plot(mod, what = "classification")
plot(mod, what = "density")
boot <- MclustBootstrap(mod, nboot = 25, type = "bs")#25 for now to speed up
接下来我想创建一个新的数据框来显示每个观察的混合概率，所以我将原始的dat1 与mod$z 结合起来。
probs <- mod$z
colnames(probs) <- paste0("Prob", 1:mod$G)
probs <- cbind(dat1, probs)
probs <- cbind(probs, cluster = mod$classification)
现在我要收集：
1) 每个集群中每个变量的平均值，以及 2) 每个平均值的标准误差（来自 summary(boot, what = "SE")$mean）。我将使用这些来创建下面的图，并返回一个显示平均值和 SE 值的表。
newdat <- cbind(dat1, cluster = mod$classification)
a <-
newdat%>%
dplyr::select_if(is.numeric)%>%
dplyr::group_by(cluster)%>%
summarise_all(mean)%>%
pivot_longer(-c("cluster"), names_to = "Vars", values_to = "mean")
b<-
as.data.frame(cbind(cluster=1:mod$G, t(summary(boot, what = "se")$mean)))%>%
pivot_longer(., -c("cluster"), names_to = "Vars", values_to = "SE")
a<-dplyr::mutate(a, "SE" = b$SE)
ggplot(a, aes(x=Vars, y=mean, group=factor(cluster))) +
geom_line(aes(colour=factor(cluster)))+
geom_point()+
geom_errorbar(aes(ymin=mean-SE, ymax=mean+SE), width = .1)+
ggtitle("Mean by cluster") +
expand_limits(y=0) +
scale_y_continuous(breaks=0:20*4) +
labs(x = "Var", y = "Cluster Average")+
theme_bw() +
theme(legend.justification=c(1,0),
legend.position=c(1,0))
library(knitr)
kable(a)
最终我想编写一个函数来执行每个步骤并立即返回输出。它将应用于结构类似于dat1 的数据帧，但它们将具有不同数量的var 列。我还需要指定在创建mod 时要使用多少个集群以及要使用哪个模型。什么是我收集存储在a 中的信息的更好方法，以便可以将相同的过程应用于任何变量组合（var. 列）和混合组件（G 在mod 中指定）从函数内部？
【问题讨论】：
-
你看过 broom 和 tidymodels 包了吗？我不完全清楚您想要做什么/代码的哪一部分需要更改，但这些可能会有所帮助
-
@camille 我试图提供我在做什么的背景以及为什么，但这可能会使我的具体问题不那么明显。我将尝试澄清。我的问题涉及从对象
mod和boot（对象a的创建）中提取均值和标准误差的方法。我本质上需要一种不同的方法来创建a而无需创建b，因为用于创建b的矩阵转置对于变量和混合分量的某些组合存在问题。 -
仅供参考，tidymodels 包可能是一条不错的路线。我使用
a<- data.frame(cbind(tidy(mod), t(summary(boot, what = "se")$mean)))获取我需要的信息，现在我只需要以与上面a相同的格式获取它
标签： r function dplyr functional-programming mclust