如何从 me.weighted() 计算单个和组合分量的概率密度

如何解决如何从 me.weighted() 计算单个和组合分量的概率密度

我正在使用 Mclust 来估计组件成员资格的概率，但“密度”不包含在 me.weighted() 的输出中。因此，我无法绘制概率密度。下面的代码很长，因为我想清楚地说明我的目的和问题，但我清楚地指出了我的问题/问题出现的地方。我的最后一段代码是我对解决方案的尝试，但它可能只会凸显我对概率密度的无知。

在这个研究项目中，我的第一个目标是计算 1 岁鱼丰度指数，以供后续分析。为此，我想估计特定长度的 1 岁鱼的比例（即年龄长度键）。可以合理地假设较小的模式主要是 1 岁的鱼，而较大的模式是 2 岁以上的鱼。我的数据是鱼体长度（叉长，厘米）和丰度占总数的比例（即加权单变量）。请注意，省略了一些比例较小的外围大长度；因此， sum(dat.df$proportions)

我在这里的具体目的是说明叠加在鱼大小组成上的概率密度，它反映了两个年龄组。基本上，在 ggplot 代码的最后一部分中，我想用概率密度替换每个（红色）或任一（绿色）组件的成员资格估计概率，因为它会在我的手稿中形成一个很好的、信息丰富的数字。

我已阅读相关文章（墨菲；Scrucca 等；Mignan；R-Bloggers 等），但没有找到答案。

因此，我非常感谢有关如何计算每个组件的概率密度以及组件组合概率密度的任何帮助。

套餐

library(ggplot2)
library(mclust)

数据

dat.df <- data.frame(flcm = 15:33,proportion = c(0.0043,0.0114,0.0296,0.0519,0.0540,0.0403,0.0294,0.0152,0.0257,0.0793,0.1458,0.1505,0.1277,0.0909,0.0389,0.0308,0.0121,0.0101,0.0085),z1 = c(rep(1,9),rep(0,10)),z2 = c(rep(0,rep(1,10)))

绘图数据

ggplot()+
  geom_bar(aes(x=dat.df$flcm,y=dat.df$proportion),fill = "gray",position="dodge",stat="identity")+
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()

无权重（即忽略 dat.df$proportion）

拟合没有权重的混合模型

mod1 <- densityMclust(dat.df[,"flcm"],modelName = "V")

绘制概率密度

plot(mod1,what = "density",data = dat.df$flcm,breaks = 5)

带权重（即包括 dat.df$proportion）

使用权重重新拟合模型

mod1_w <- me.weighted(modelName = "V",z = cbind(dat.df$z1,dat.df$z2),weights = dat.df$proportion)

使用估计的分数隶属关系绘制数据（更新 z）

ggplot()+
  geom_bar(aes(x=dat.df$flcm,stat="identity")+
  geom_line(aes(x = dat.df$flcm,y = (mod1_w$z[,1] * dat.df$proportion)),color = "red") +
  geom_line(aes(x = dat.df$flcm,2] * dat.df$proportion)),color = "red") +
    geom_line(aes(x = dat.df$flcm,1] * dat.df$proportion) +
                  mod1_w$z[,2] * dat.df$proportion),color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()

绘制概率密度 - 这里是我的问题/疑问出现的地方

plot(mod1_w,breaks = 5)`

这是我尝试的解决方案。基本上，对于每个组件（年龄 1、年龄 2），将概率相乘并按比例缩放：

#age1 probability density
age1 <- mod1_w$z[,1]* #probability of age1 membership multiplied by
  dnorm(dat.df$flcm,mod1_w$parameters$mean[1],#probability of flcm given age1
        mod1_w$parameters$variance$sigmasq[1])* 
  sum(mod1_w$z[,1]*mod1_w$weights) #and scaled to proportional abundance of age1

#age2 probability density
age2 <- mod1_w$z[,2]* #probability of age2 membership multiplied by
  dnorm(dat.df$flcm,mod1_w$parameters$mean[2],mod1_w$parameters$variance$sigmasq[2])* #probability of flcm given age2
  sum(mod1_w$z[,2]*mod1_w$weights) #and scaled to proportional abundance of age2

#combined ages probability density
age_all <- age1 + age2

#looks bad - the probability densities don't correspond well with proportional abundance
ggplot()+
  geom_bar(aes(x=dat.df$flcm,y = age1),y = age2),y = age_all),color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()

解决方法

我正在发布我自己问题的解决方案；希望这会帮助其他人。基本上，我切换到包 mxdist，它给了我所需的输出，如以下代码所示。

library(mxdist)

#input data (dat.df is created by code in my original question)
dat.mx <- as.mixdata(dat.df[,1:2])
#preliminary plot
plot(dat.mx)
#define initial parameters
dat_parms <- data.frame(pi=c(0.3,0.7),mu=c(18,26),sigma=c(2,3))
#fit the model
fit1 <- mix(dat.mx,dat_parms,"gamma",constr=mixconstr(consigma="CCV"))
#plot default
plot(fit1)
#replot using ggplot for greater flexibility over appearance
z <- fitted(fit1)
dat.mx[dat.mx$flcm == "Inf","flcm"] <- 34
ggplot()+
  geom_bar(aes(x=dat.mx$flcm,y=dat.mx$proportion),fill = "gray",position="dodge",stat="identity")+
  geom_line(aes(x = dat.mx$flcm,y = z$joint[,1]),color = "red") +
  geom_line(aes(x = dat.mx$flcm,2]),color = "red") +
    geom_line(aes(x = dat.mx$flcm,y = z$mixed),color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()
#conditional probabilities are output
z$conditprob