多项式hmm使用R

问题描述

我有以下类型的数据：

df <- data.frame(id = rep(1,each = 40),year = seq(1961,2000),x1 = rbinom(40,size = 1,prob = 1 - 0.6) * rpois(40,lambda = 4),X2 = rbinom(40,prob = 1 - 0.7) * rpois(40,X3 = rbinom(40,lambda = 5),X4 = rbinom(40,lambda = 6))

正如您在我的数据中看到的，有四个计数变量。我想估计一个多项式 HMM，因为我期望有一个潜在变量 C，它有 3 个可能的状态影响 Pr(X_t=xt)（假设每个向量 X_t 在马尔可夫链 C_t 上是相互独立的条件） .例如，我希望如果 C_t=1 我们会观察到一个更像这样的 X_t 向量 (4,1,0)，如果 C_T=2 一个更像这样的 X_t 向量 (0,0)如果 C=3，则更有可能观察到这样的 X_t 向量 (0,5)。

我还没有找到能够估计这种类型模型的包，所以目前，我使用的是 depmixS4 包。

library(depmixS4)
mod<-depmix(list(X1 ~ 1,X2 ~ 1,X3 ~ 1,X4 ~ 1),data=df,nstates=3,family=list(poisson(),poisson(),poison()))

但是，根据我的理论预期，我不确定这是否是正确的模型。我可以以不同的方式使用 depmix 以更适合我的模型吗？

解决方法

您可以简单地对响应使用多项分布。我假设您的意思是让 X1,...,X4 指代单个分类变量的四个级别？然后这些变量中的每一个都包含特定年份的级别计数？一种选择是拟合以下模型：

> library(depmixS4)
Loading required package: nnet
Loading required package: MASS
Loading required package: Rsolnp
Loading required package: nlme
> set.seed(1)
> df <- data.frame(id = rep(1,each = 40),+                  year = seq(1961,2000),+                  X1 = rbinom(40,size = 1,prob = 1 - 0.6) * rpois(40,lambda = 4),+                  X2 = rbinom(40,prob = 1 - 0.7) * rpois(40,+                  X3 = rbinom(40,lambda = 5),+                  X4 = rbinom(40,lambda = 6))
> # matrix for single multinomial response variable
> X <- as.matrix(df[,c("X1","X2","X3","X4")])
> 
> # formulate model
> mod<-depmix(X ~ 1,data=df,nstates=3,+             family=multinomial("identity"))
> 
> # fit model
> fmod <- fit(mod)
converged at iteration 22 with logLik: -161.5714 
> 
> # show results
> summary(fmod)
Initial state probabilities model 
pr1 pr2 pr3 
  0   0   1 

Transition matrix 
        toS1  toS2  toS3
fromS1 0.290 0.355 0.356
fromS2 0.132 0.328 0.540
fromS3 0.542 0.191 0.267

Response parameters 
Resp 1 : multinomial 
    Re1.pr1 Re1.pr2 Re1.pr3 Re1.pr4
St1   0.092    0.56   0.288   0.061
St2   0.033    0.00   0.423   0.544
St3   0.608    0.00   0.392   0.000

hidden-markov-models multinomial r r