i具有从一个状态到另一个状态的过渡计数的矩阵，我想计算最大似然估计，标准误差和标准偏差。" 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过渡计数矩阵，并确定最大似然估计器，标准误差（置信区间）和标准偏差。有人可以效法吗？
我不确定如何使用markovchain软件包做到这一点，但是您可以编写一个可以产生类似结果的函数。下面的函数通过MLE生成了过渡概率的估计值，但是标准误差，通过参数bootstrap的上限和上置信度边界。您可以使用标准误差和估计来制作正常理论置信区间（ci = "normal"）或使用参数引导分布的相关分位数（ci = "quantile"）。parallel标志标识是否应考虑矩阵的级别。如果有序，该模型会使用有序登录来计算概率。如果未排序（即parallel = FALSE），则使用多项式登录来计算概率。
如果我正确理解，您是从看起来像createSequenceMatrix(stringchar = rain$rain)输出的东西开始的，而没有生成它的基础值的基础序列。这就是它的工作方式。
- 生成数据。
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)
- 使功能估算模型并打印结果。
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 Probabilitiesn")
print(x$estimate, digits=digits)
cat("nStandard Errors of Transition Probabilitiesn")
print(x$se, digits=digits)
cat("nLower Confidence Boundsn")
print(x$lower, digits=digits)
cat("nUpper Confidence Boundsn")
print(x$upper, digits=digits)
cat('n')
print(x$logLik)
}
- 使用我生成的函数估算模型：
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)
- 比较脚趾
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软件包（v2.0.1）在2022-12-15创建
这两个模型的对数可能是相同的，这使人们对基础估计过程相似。标准错误是不同的 - 原始函数不使用参数bootstrap。我不确定他们做什么，因为该功能以估计模型是用C 编写的，因此我没有调查源代码。如果这种差异是不可接受的，也许其他人将能够靠近。