Расчет вероятностей переходов при изменении продолжительности цикла Маркова методом спектрального разложения исходной матрицы переходов

Задача из разряда встречающихся в практике при моделировании, опишу ее в общих чертах. Выполняем моделирование исходов на основании раннее выполненного и опубликованного исследования. В модели 3 состояния и более, есть матрица переходов между состояниями с продолжительностью цикла Маркова 1 год. Но в нашем исследовании, например по причине особенностей режима дозирования изучаемого препарата, нужно пересчитать матрицу переходов на длину цикла, равную 1 месяц. Казалось бы, задачка решается очень просто, допуская равномерную плотность инцеденса пересчитываем вероятность (год) в частоту (год), частоту делим на 12 и трансформируем обратно в вероятность (мес.). Как описано, например в Fleurence RL, Hollenbeak CS. Rates and probabilities in economic modelling: transformation, translation and appropriate application. Pharmacoeconomics. 2007;25(1):3-6. doi: 10.2165/00019053-200725010-00002. PMID: 17192114. В нашем случае:

P _мес. = 1 – e^^-z/n = 1 – e^^ln(1-P)/n = 1 – (1 – P)^^1/n = 1 – (1 – P)^^1/12

P _мес – искомая вероятность (мес.)
Р – исходная вероятность (год)
е – число Эйлера
z – частота (год)
n – число месяцев в году

Для примера возьмем часто применяемую модель Маркова с 3-мя состояниями: стабильное, прогрессия, смерть.

Соответственно, получается следующая матрица переходов:

Состояние	stab	prog	death	проверка
stab	1 – (P stab-prog) – (P stab-death)	P stab-prog	P stab-death	1
prog	0	1 – (P prog-death)	P prog-death	1
death	0	0	1	1

А именно, матрица переходов с длиной цикла 1 год:

Состояние	stab	prog	death	проверка
stab	0,3	0,2	0,5	1
prog	0	0,4	0,6	1
death	0	0	1	1

Сделаем пересчет по ранее указанной формуле: P _мес. = 1 – e^^-z/n = 1 – e^^ln(1-P)/n = 1 – (1 – P)^^1/n = 1 – (1 – P)^^1/12 и получим матрицу переходов с длиной цикла 1 мес.:

Состояние	stab	prog	death	проверка
stab	0,925450843	0,01842347	0,056125687	1
prog	0	0,926484872	0,073515128	1
death	0	0	1	1

Рассчитываем две модели Маркова: с длиной цикла 1 год (исходная матрица переходов) и длиной цикла 1 мес. (пересчитанная матрица переходов). Результат получается достаточно плачевный, пересчет вышеприведенным и очевидно напрашивающимся методом приводит к существенной ошибке (см. таблицу ниже).

		Исходная матрица			Пересчитанная матарица
год	мес	stab	prog	death	stab	prog	death
0	0	1	0	0	1	0	0
	1				0,925451	0,018423	0,056126
	2				0,856459	0,034119	0,109422
	3				0,792611	0,04739	0,159999
	4				0,733522	0,058509	0,207969
	5				0,678839	0,067721	0,25344
	6				0,628232	0,075249	0,296519
	7				0,581398	0,081292	0,33731
	8				0,538055	0,086027	0,375918
	9				0,497944	0,089615	0,412441
	10				0,460822	0,092201	0,446976
	11				0,426468	0,093913	0,479619
1	12	0,3	0,2	0,5	0,394676	0,094866	0,510459
	13				0,365253	0,095163	0,539584
	14				0,338024	0,094896	0,56708
	15				0,312824	0,094148	0,593028
	16				0,289503	0,09299	0,617507
	17				0,267921	0,091487	0,640592
	18				0,247948	0,089698	0,662355
	19				0,229464	0,087671	0,682865
	20				0,212357	0,085454	0,702189
	21				0,196526	0,083084	0,72039
	22				0,181875	0,080597	0,737528
	23				0,168317	0,078022	0,753661
2	24	0,09	0,14	0,77	0,155769	0,075388	0,768844
	25				0,144156	0,072715	0,783128
	26				0,13341	0,070025	0,796565
	27				0,123464	0,067335	0,809201
	28				0,11426	0,06466	0,82108
	29				0,105742	0,062011	0,832247
	30				0,097859	0,059401	0,84274
	31				0,090564	0,056837	0,852599
	32				0,083812	0,054327	0,861861
	33				0,077564	0,051877	0,870559
	34				0,071782	0,049492	0,878726
	35				0,06643	0,047177	0,886393
3	36	0,027	0,074	0,899	0,061478	0,044932	0,89359
	37				0,056895	0,042762	0,900343
	38				0,052654	0,040666	0,90668
	39				0,048728	0,038647	0,912625
	40				0,045096	0,036703	0,918201
	41				0,041734	0,034836	0,92343
	42				0,038623	0,033044	0,928334
	43				0,035743	0,031326	0,932931
	44				0,033079	0,029682	0,93724
	45				0,030613	0,028109	0,941278
	46				0,028331	0,026607	0,945063
	47				0,026218	0,025173	0,948609
4	48	0,0081	0,035	0,9569	0,024264	0,023805	0,951931
	49				0,022455	0,022502	0,955043
	50				0,020781	0,021262	0,957957
	51				0,019232	0,020081	0,960687
	52				0,017798	0,018959	0,963243
	53				0,016471	0,017893	0,965635
	54				0,015243	0,016881	0,967875
	55				0,014107	0,015921	0,969972
	56				0,013055	0,015011	0,971934
	57				0,012082	0,014148	0,97377
	58				0,011181	0,01333	0,975488
	59				0,010348	0,012556	0,977096
5	60	0,00243	0,01562	0,98195	0,009576	0,011824	0,9786

Решение было найдено в Chhatwal J, Jayasuriya S, Elbasha EH. Changing Cycle Lengths in State-Transition Models: Challenges and Solutions. Medical Decision Making. 2016;36(8):952-964. doi:10.1177/0272989X16656165. В качестве такового возможно использование спектрального разложение матрицы переходов (или разложение матрицы на основе собственных векторов или eigen decomposition). Спектральная декомпозиция не всегда приводит к ожидаемым результатам – могут быть получены отрицательные вероятности перехода или комплексные числа (что делать в этом случае – см. вышеуказанные источник). Однако, в данном примере ограничимся ситуацией, когда применение спектральной декомпозиции приводит к необходимому результату и разберем выполнение данной процедуры в программной среде R. В сухом остатке нам необходимо получить корень 12 степени (12 месяцев в одном году) исходной матрицы переходов или, другими словами, возвести исходную матрицу в 1/12 степень. Формула решения выглядит следующим образом:

P^1/12 = V * D^1/n * V^-1 = V * D^1/12 * V^-1

P^1/12 – искомая матрица переходов с шагом один месяц
V — матрица, столбцы которой являются собственными векторами исходной матрицы
D — диагональная матрица с соответствующими собственными значениями на главной диагонали
V^-1 – обратная матрица векторов

В программной среде R решением выглядит следующим образом:

1. В переменную mtrx записываем нашу исходную матрицу

input <- c(0.3, 0.2, 0.5,
           0, 0.4, 0.6,
           0, 0, 1)
mtrx <- matrix(input, nrow = 3, ncol = 3, byrow = T)

2. В переменную Х записываем результаты спектральной декомпозиции исходной матрицы:

X <- eigen(mtrx)
#> X
#eigen() decomposition
#$values  – это собственные значения по диагонали матрицы
#[1] 1.0 0.4 0.3

#$vectors – это собственные векторы матрицы
#          [,1]      [,2] [,3]
#[1,] 0.5773503 0.8944272    1
#[2,] 0.5773503 0.4472136    0
#[3,] 0.5773503 0.0000000    0

3. Применяем формулу P^1/12 = V * D^1/n * V^-1 = V * D^1/12 * V^-1

X$vectors %*% diag(X$values)^(1/12) %*% solve(X$vectors)

Где:

X$vectors — собственные векторы матрицы

diag(X$values) — диагональная матрица с соответствующими собственными значениями на главной диагонали

solve(X$vectors) — обратная матрица векторов

%*% — знак умножения матриц

Итог:

#> X$vectors %*% diag(X$values)^(1/12) %*% solve(X$vectors)
#          [,1]       [,2]       [,3]
#[1,] 0.9045379 0.04389393 0.05156816
#[2,] 0.0000000 0.92648487 0.07351513
#[3,] 0.0000000 0.00000000 1.00000000

Запишем итог в красивую таблицу и применим эту матрицу в Эксель.

Состояние	stab	prog	death	проверка
stab	0,9045379	0,04389393	0,05156816	1
prog	0	0,92648487	0,07351513	1
death	0	0	1	1

Результаты с использованием матрицы переходов шагом 1 месяцев в результате применения спектрального разложения матрицы.

		Исходная матрица			Пересчитанная матарица
год	мес	stab	prog	death	stab	prog	death
0	0	1	0	0	1	0	0
	1				0,904538	0,043894	0,051568
	2				0,818189	0,080371	0,10144
	3				0,740083	0,110376	0,149541
	4				0,669433	0,134747	0,19582
	5				0,605527	0,154225	0,240248
	6				0,547723	0,169466	0,282811
	7				0,495436	0,181049	0,323515
	8				0,44814	0,189486	0,362373
	9				0,40536	0,195227	0,399413
	10				0,366664	0,198667	0,434669
	11				0,331661	0,200157	0,468182
1	12	0,3	0,2	0,5	0,3	0,2	0,5
	13				0,271361	0,198465	0,530173
	14				0,245457	0,195786	0,558757
	15				0,222025	0,192167	0,585808
	16				0,20083	0,187785	0,611385
	17				0,181658	0,182795	0,635546
	18				0,164317	0,177331	0,658352
	19				0,148631	0,171507	0,679862
	20				0,134442	0,165423	0,700135
	21				0,121608	0,159163	0,719229
	22				0,109999	0,1528	0,737201
	23				0,099498	0,146395	0,754107
2	24	0,09	0,14	0,77	0,09	0,14	0,77
	25				0,081408	0,133658	0,784933
	26				0,073637	0,127406	0,798957
	27				0,066607	0,121272	0,812121
	28				0,060249	0,11528	0,824471
	29				0,054497	0,10945	0,836053
	30				0,049295	0,103796	0,846909
	31				0,044589	0,098329	0,857082
	32				0,040333	0,093057	0,86661
	33				0,036482	0,087987	0,875531
	34				0,033	0,08312	0,883881
	35				0,029849	0,078458	0,891693
3	36	0,027	0,074	0,899	0,027	0,074	0,899
	37				0,024423	0,069745	0,905832
	38				0,022091	0,06569	0,912219
	39				0,019982	0,06183	0,918187
	40				0,018075	0,058162	0,923763
	41				0,016349	0,054679	0,928971
	42				0,014789	0,051377	0,933834
	43				0,013377	0,048249	0,938374
	44				0,0121	0,045289	0,942611
	45				0,010945	0,042491	0,946564
	46				0,0099	0,039848	0,950252
	47				0,008955	0,037353	0,953692
4	48	0,0081	0,035	0,9569	0,0081	0,035	0,9569
	49				0,007327	0,032783	0,959891
	50				0,006627	0,030694	0,962678
	51				0,005995	0,028729	0,965277
	52				0,005422	0,02688	0,967698
	53				0,004905	0,025142	0,969954
	54				0,004437	0,023509	0,972055
	55				0,004013	0,021975	0,974012
	56				0,00363	0,020536	0,975834
	57				0,003283	0,019185	0,977531
	58				0,00297	0,017919	0,979111
	59				0,002686	0,016732	0,980581
5	60	0,00243	0,01562	0,98195	0,00243	0,01562	0,98195

Так намного лучше или даже идеально!

На рисунке ниже приведена визуализация двух подходов: неправильного («стандартного») и правильного (eigen decomposition или спектральное разложение исходной матрицы переходов).

Дорогие читатели, в этом материале мы познакомились с методом, который может позволить корректно пересчитать вероятности перехода при изменении длины цикла Маркова для моделей с 3 и более состояниями.

Дорогой читатель, буду рад твоим комментариям, предложениям и вопросам!