0

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

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

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

我已阅读相关文章（Murphy；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,9), 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", 
            data = dat.df$flcm,
            z = cbind(dat.df$z1, dat.df$z2),
            weights = dat.df$proportion)

用估计的分数成员绘制数据（更新的 z）

ggplot()+
  geom_bar(aes(x=dat.df$flcm, y=dat.df$proportion),
           fill = "gray", position="dodge", 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, 
                y = (mod1_w$z[,2] * dat.df$proportion)), 
            color = "red") +
    geom_line(aes(x = dat.df$flcm, 
                y = (mod1_w$z[,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, what = "density", data = dat.df$flcm, breaks = 5)`

这是我尝试的解决方案。基本上，对于每个成分（age1、age2），乘以概率并缩放到比例丰度：

#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=dat.df$proportion),
           fill = "gray", position="dodge", stat="identity")+
  geom_line(aes(x = dat.df$flcm, 
                y = age1), 
            color = "red") +
  geom_line(aes(x = dat.df$flcm, 
                y = age2), 
            color = "red") +
    geom_line(aes(x = dat.df$flcm, 
                y = age_all), 
              color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()
4

1 回答 1

0

我正在发布我自己的问题的解决方案；希望这会对其他人有所帮助。基本上，我切换到包 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, 
                y = z$joint[,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
于 2021-02-25T21:12:35.273 回答