r - 计算离散时间马尔可夫链的标准误差和标准差

Question

我有一个从一种状态到另一种状态的转换计数矩阵，我想计算最大似然估计、标准误差和标准差。 “markovchain”包有一个例子，但数据是一个序列。我的数据是从 155 家公司的平衡面板数据集中获得的，因此他们提供的示例代码对我不起作用。

这是我遵循的示例:

 data(rain)
 rain$rain[1:10]

 [1] "6+"  "1-5" "1-5" "1-5" "1-5" "1-5" "1-5" "6+"  "6+"  "6+"

 #obtaining the empirical transition matrix
 createSequenceMatrix(stringchar = rain$rain)
      0 1-5  6+ 
0   362 126  60 
1-5 136  90  68 
6+   50  79 124

#fitting the DTMC by MLE
alofiMcFitMle <- markovchainFit(data = rain$rain, method = "mle", name = "Alofi")
alofiMcFitMle
$estimate
Alofi
A  3 - dimensional discrete Markov Chain defined by the following states:
 0, 1-5, 6+
 The transition matrix  (by rows)  is defined as follows:
            0       1-5        6+
0   0.6605839 0.2299270 0.1094891
1-5 0.4625850 0.3061224 0.2312925
6+  0.1976285 0.3122530 0.4901186
$standardError
             0        1-5         6+
0   0.03471952 0.02048353 0.01413498
1-5 0.03966634 0.03226814 0.02804834
6+  0.02794888 0.03513120 0.04401395
$confidenceInterval
$confidenceInterval$confidenceLevel
[1] 0.95
$confidenceInterval$lowerEndpointMatrix
            0       1-5         6+
0   0.6034754 0.1962346 0.08623909
1-5 0.3973397 0.2530461 0.18515711
6+  0.1516566 0.2544673 0.41772208
$confidenceInterval$upperEndpointMatrix
            0       1-5        6+
0   0.7176925 0.2636194 0.1327390
1-5 0.5278304 0.3591988 0.2774279
6+  0.2436003 0.3700387 0.5625151
$logLikelihood
[1] -1040.419

因为我已经有了一个计数数据矩阵，所以我不能使用上面的代码。我只想采用我的 6x6 转换计数矩阵并确定最大似然估计量、标准误差(置信区间)和标准差。有没有人有我可以效仿的例子？

Answer 1

我不确定您如何使用 markovchain 包执行此操作，但您可以编写一个函数来产生类似的结果。下面的函数通过 MLE 生成转移概率的估计值，但通过参数 Bootstrap 生成标准误差、置信下限和置信上限。您可以使用标准误差和估计值来制作正态理论置信区间 (ci = "normal") 或使用参数自举分布的相关分位数 (ci = "quantile")。 parallel 标志标识矩阵的级别是否应被视为有序的。如果订购，该模型使用订购登录来计算概率。如果无序(即 parallel = FALSE)，则使用多项登录来计算概率。

如果我理解正确，您从看起来像 createSequenceMatrix(stringchar = rain$rain) 的输出的内容开始，但没有生成它的底层值序列。这是它的工作原理。

1. 生成数据。

library(markovchain)
#> Package:  markovchain
#> Version:  0.9.0
#> Date:     2022-07-01
#> BugReport: https://github.com/spedygiorgio/markovchain/issues
data(rain)
mat <-  createSequenceMatrix(stringchar = rain$rain)

1. 编写函数来估计模型并打印结果。

mcmat_est <- function(mat, parallel = FALSE, ci = c("normal", "quantile"), conf_level=.95, R=2500){
  ci <- match.arg(ci)
  if(!is.null(rownames(mat))){
    vals <- rownames(mat)
  }else{
    vals <- 1:nrow(mat)
  }
  state_vals <- lapply(vals, \(i)data.frame(vals = rep(vals, mat[i, ])))
  if(parallel){
    state_vals <- lapply(state_vals, \(x){x$vals <- factor(x$vals, levels=vals, ordered=TRUE); x})
    mods <- lapply(state_vals, \(x)MASS::polr(vals ~ 1, data=x, Hess = TRUE))
    draws <- lapply(mods, \(m)cbind(-Inf, MASS::mvrnorm(R, m$zeta, vcov(m)), Inf))
    qs <- lapply(draws, \(d)plogis(d))
    probs<- lapply(qs, \(x)x[,-1] - x[,-ncol(x)])
    ests <- lapply(mods, \(m)plogis(c(-Inf, m$zeta, Inf)))
    ests <- lapply(ests, \(e)e[-1]-e[-length(e)])
  }else{
    state_vals <- lapply(state_vals, \(x){x$vals <- factor(x$vals, levels=vals); x})
    mods <- lapply(state_vals, \(x)nnet::multinom(vals ~ 1, data=x, Hess = TRUE, trace=FALSE))
    draws <- lapply(mods, \(m)MASS::mvrnorm(R, c(0, coef(m)), 
                                            cbind(0, rbind(0, vcov(m)))))
    probs <- lapply(draws, \(d)prop.table(exp(d), 1))
    ests <- lapply(mods, \(m)exp(c(0, coef(m))))
    ests <- lapply(ests, \(e)e/sum(e))
  }
  logLik <- sum(sapply(mods, logLik))
  nobs <- sum(sapply(mods, \(x)attr(logLik(x), "nobs")))
  df <- sum(sapply(mods, \(x)attr(logLik(x), "df")))
  attr(logLik, "nobs") <- nobs
  attr(logLik, "df") <- df
  class(logLik) <- "logLik"
  est <- do.call(rbind, ests)
  se <- do.call(rbind, lapply(probs, \(p)apply(p, 2, sd)))
  crit_z <- 1-(1-conf_level)/2
  if(ci == "normal"){
    lwr <- est - qnorm(crit_z)*se
    upr <- est + qnorm(crit_z)*se
  }else{
    lwr <- do.call(rbind, lapply(probs, \(p)apply(p, 2, quantile, 1-crit_z)))
    upr <- do.call(rbind, lapply(probs, \(p)apply(p, 2, quantile, crit_z)))
  }
  rownames(est) <- rownames(se) <- rownames(lwr) <- rownames(upr) <- rownames(mat)
  colnames(est) <- colnames(se) <- colnames(lwr) <- colnames(upr) <- colnames(mat)
  res <- list(estimate = est, se = se, lower=lwr, upper=upr, logLik=logLik)
  class(res) <- "mcest"
  res
}
print.mcest <- function(x, ..., digits=3){
  cat("\nEstimated Transition Probabilities\n")
  print(x$estimate, digits=digits)
  cat("\nStandard Errors of Transition Probabilities\n")
  print(x$se, digits=digits)
  cat("\nLower Confidence Bounds\n")
  print(x$lower, digits=digits)
  cat("\nUpper Confidence Bounds\n")
  print(x$upper, digits=digits)
  cat('\n')
  print(x$logLik)
}

1. 使用我生成的函数估计模型:

out <- mcmat_est(mat)  
out
#> 
#> Estimated Transition Probabilities
#>         0   1-5    6+
#> 0   0.661 0.230 0.109
#> 1-5 0.463 0.306 0.231
#> 6+  0.198 0.312 0.490
#> 
#> Standard Errors of Transition Probabilities
#>          0    1-5     6+
#> 0   0.0203 0.0181 0.0135
#> 1-5 0.0281 0.0269 0.0248
#> 6+  0.0245 0.0289 0.0309
#> 
#> Lower Confidence Bounds
#>         0   1-5    6+
#> 0   0.621 0.194 0.083
#> 1-5 0.407 0.253 0.183
#> 6+  0.150 0.256 0.430
#> 
#> Upper Confidence Bounds
#>         0   1-5    6+
#> 0   0.700 0.265 0.136
#> 1-5 0.518 0.359 0.280
#> 6+  0.246 0.369 0.551
#> 
#> 'log Lik.' -1040.419 (df=6)

1. 比较 markovchainFit() 输出。

fit <- markovchainFit(data = rain$rain, method = "mle", name = "Alofi")
fit
#> $estimate
#> Alofi 
#>  A  3 - dimensional discrete Markov Chain defined by the following states: 
#>  0, 1-5, 6+ 
#>  The transition matrix  (by rows)  is defined as follows: 
#>             0       1-5        6+
#> 0   0.6605839 0.2299270 0.1094891
#> 1-5 0.4625850 0.3061224 0.2312925
#> 6+  0.1976285 0.3122530 0.4901186
#> 
#> 
#> $standardError
#>              0        1-5         6+
#> 0   0.03471952 0.02048353 0.01413498
#> 1-5 0.03966634 0.03226814 0.02804834
#> 6+  0.02794888 0.03513120 0.04401395
#> 
#> $confidenceLevel
#> [1] 0.95
#> 
#> $lowerEndpointMatrix
#>             0       1-5        6+
#> 0   0.5925349 0.1897800 0.0817850
#> 1-5 0.3848404 0.2428780 0.1763188
#> 6+  0.1428496 0.2433971 0.4038528
#> 
#> $upperEndpointMatrix
#>             0       1-5        6+
#> 0   0.7286330 0.2700740 0.1371931
#> 1-5 0.5403296 0.3693669 0.2862663
#> 6+  0.2524073 0.3811089 0.5763843
#> 
#> $logLikelihood
#> [1] -1040.419

^{由 reprex package 创建于 2022-12-15 (v2.0.1)}

这两个模型的对数似然相同，这让我们有信心基本估计过程是相似的。标准错误不同——原始函数不使用参数 Bootstrap 。我不确定他们做了什么，因为估计模型的函数是用 C++ 编写的，而且我没有调查源代码。如果这种差异是 Not Acceptable ，也许其他人将能够更接近。