MSM Test Notebook 2 - Hidden Markov Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a5fdd52f",
   "metadata": {},
   "source": [
    "## The Hidden Markov model\n",
    "\n",
    "Building some first multi-state models using the supplied \"fev\" data.\n",
    "\n",
    "The \"fev\" data is a dataset of repeated measurements of forced expiratory volume in 1 second, in recipients of lung transplants. FEV measurements are used to diagnose bronchiolitis obliterans syndrome (BOS), a chronic deterioration in lung function.  FEV is measured as a % baseline value for each individual in the 1st 6 months after transplant, 100% baseline at 6 months.\n",
    "\n",
    "Stages of BOS defined as:\n",
    "\n",
    "- No BOS (>= 80% baseline FEV)\n",
    "- Mild BOS (sustained decrease below 80% baseline FEV)\n",
    "- Moderate BOS (sustained decrease below 65% baseline FEV)\n",
    "- Severe BOS (sustained decrease below 50% baseline FEV)\n",
    "- Death\n",
    " \n",
    " \n",
    "In a hidden Markov model (HMM), the state of the Markov chain is not observable. Instead, the observed data (Y) are governed by a probability distribution conditional on the unobserved states (S). P(Y|S) is the probability density matrix of the outcome conditional on the hidden state (aka the \"emission\" distribution) and P(S|S-1) is the transition probability matrix of the hidden Markov chain.  Distributions must be univariate. Restricted to situations where only 1 observation is made conditionally on an underlying Markov process. The distribution of P(S1) - the initial state - can be estimated from the data or fixed at plausible values. The distribution may depend on covariates (logistic regression). \n",
    "\n",
    "\n",
    "In medicine, these models can represent chronic staged diseases which can only be diagnosed by an error-prone marker.\n",
    "\n",
    "\n",
    "We will model BOS using the HMM because onset values fluctuate, and entry into states can be sudden.\n",
    "\n",
    "\n",
    "The HMM will have THREE states: \n",
    "\n",
    "- State 1: No BOS\n",
    "- State 2: BOS\n",
    "- State 3: Death\n",
    "\n",
    "The distribution of FEV is normal in states 1 & 2\n",
    "\n",
    "The distribution of Death is without error and identified as \"identity\" distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8b900277",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“package ‘msm’ was built under R version 3.6.3”"
     ]
    }
   ],
   "source": [
    "library(\"msm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b784d743",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>ptnum</th><th scope=col>days</th><th scope=col>fev</th><th scope=col>acute</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>1     </td><td> 190  </td><td> 95.16</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 249  </td><td>112.90</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 278  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 369  </td><td>109.68</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 466  </td><td>112.90</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 557  </td><td>114.52</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 606  </td><td> 91.94</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 613  </td><td> 83.87</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 617  </td><td>101.61</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 634  </td><td> 85.48</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 637  </td><td> 85.48</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 666  </td><td> 82.26</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 669  </td><td> 61.29</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 670  </td><td> 80.65</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 672  </td><td> 67.74</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 678  </td><td> 95.16</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td> 693  </td><td> 87.10</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 812  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 902  </td><td> 96.77</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 932  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td> 991  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1098  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1139  </td><td> 93.55</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1152  </td><td>101.61</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1188  </td><td> 93.55</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1280  </td><td> 93.55</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1342  </td><td> 82.26</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td>1345  </td><td> 93.55</td><td>1     </td></tr>\n",
       "\t<tr><td>1     </td><td>1348  </td><td> 90.32</td><td>0     </td></tr>\n",
       "\t<tr><td>1     </td><td>1351  </td><td> 91.94</td><td>0     </td></tr>\n",
       "\t<tr><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "\t<tr><td>201   </td><td>370   </td><td>109.80</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>404   </td><td>109.80</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>474   </td><td>105.88</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>565   </td><td>117.65</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>670   </td><td>117.65</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>754   </td><td>117.65</td><td>0     </td></tr>\n",
       "\t<tr><td>201   </td><td>922   </td><td>121.57</td><td>0     </td></tr>\n",
       "\t<tr><td>202   </td><td>245   </td><td>118.87</td><td>0     </td></tr>\n",
       "\t<tr><td>202   </td><td>302   </td><td>122.64</td><td>0     </td></tr>\n",
       "\t<tr><td>202   </td><td>392   </td><td>128.30</td><td>0     </td></tr>\n",
       "\t<tr><td>202   </td><td>471   </td><td>120.75</td><td>0     </td></tr>\n",
       "\t<tr><td>202   </td><td>595   </td><td>105.66</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>222   </td><td>109.43</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>283   </td><td>107.55</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>346   </td><td>111.32</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>437   </td><td>109.43</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>628   </td><td>103.77</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>724   </td><td> 75.47</td><td>1     </td></tr>\n",
       "\t<tr><td>203   </td><td>752   </td><td>103.77</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>899   </td><td> 90.57</td><td>0     </td></tr>\n",
       "\t<tr><td>203   </td><td>920   </td><td> 94.34</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>203   </td><td> 94.34</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>266   </td><td>103.77</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>321   </td><td> 94.34</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>355   </td><td> 98.11</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>441   </td><td> 94.34</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>532   </td><td> 92.45</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>621   </td><td> 86.79</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>741   </td><td> 92.45</td><td>0     </td></tr>\n",
       "\t<tr><td>204   </td><td>812   </td><td>101.89</td><td>0     </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       " ptnum & days & fev & acute\\\\\n",
       "\\hline\n",
       "\t 1      &  190   &  95.16 & 0     \\\\\n",
       "\t 1      &  249   & 112.90 & 0     \\\\\n",
       "\t 1      &  278   & 101.61 & 0     \\\\\n",
       "\t 1      &  369   & 109.68 & 0     \\\\\n",
       "\t 1      &  466   & 112.90 & 0     \\\\\n",
       "\t 1      &  557   & 114.52 & 0     \\\\\n",
       "\t 1      &  606   &  91.94 & 1     \\\\\n",
       "\t 1      &  613   &  83.87 & 1     \\\\\n",
       "\t 1      &  617   & 101.61 & 1     \\\\\n",
       "\t 1      &  634   &  85.48 & 0     \\\\\n",
       "\t 1      &  637   &  85.48 & 0     \\\\\n",
       "\t 1      &  666   &  82.26 & 1     \\\\\n",
       "\t 1      &  669   &  61.29 & 1     \\\\\n",
       "\t 1      &  670   &  80.65 & 1     \\\\\n",
       "\t 1      &  672   &  67.74 & 1     \\\\\n",
       "\t 1      &  678   &  95.16 & 1     \\\\\n",
       "\t 1      &  693   &  87.10 & 0     \\\\\n",
       "\t 1      &  812   & 101.61 & 0     \\\\\n",
       "\t 1      &  902   &  96.77 & 0     \\\\\n",
       "\t 1      &  932   & 101.61 & 0     \\\\\n",
       "\t 1      &  991   & 101.61 & 0     \\\\\n",
       "\t 1      & 1098   & 101.61 & 0     \\\\\n",
       "\t 1      & 1139   &  93.55 & 0     \\\\\n",
       "\t 1      & 1152   & 101.61 & 0     \\\\\n",
       "\t 1      & 1188   &  93.55 & 0     \\\\\n",
       "\t 1      & 1280   &  93.55 & 0     \\\\\n",
       "\t 1      & 1342   &  82.26 & 1     \\\\\n",
       "\t 1      & 1345   &  93.55 & 1     \\\\\n",
       "\t 1      & 1348   &  90.32 & 0     \\\\\n",
       "\t 1      & 1351   &  91.94 & 0     \\\\\n",
       "\t ... & ... & ... & ...\\\\\n",
       "\t 201    & 370    & 109.80 & 0     \\\\\n",
       "\t 201    & 404    & 109.80 & 0     \\\\\n",
       "\t 201    & 474    & 105.88 & 0     \\\\\n",
       "\t 201    & 565    & 117.65 & 0     \\\\\n",
       "\t 201    & 670    & 117.65 & 0     \\\\\n",
       "\t 201    & 754    & 117.65 & 0     \\\\\n",
       "\t 201    & 922    & 121.57 & 0     \\\\\n",
       "\t 202    & 245    & 118.87 & 0     \\\\\n",
       "\t 202    & 302    & 122.64 & 0     \\\\\n",
       "\t 202    & 392    & 128.30 & 0     \\\\\n",
       "\t 202    & 471    & 120.75 & 0     \\\\\n",
       "\t 202    & 595    & 105.66 & 0     \\\\\n",
       "\t 203    & 222    & 109.43 & 0     \\\\\n",
       "\t 203    & 283    & 107.55 & 0     \\\\\n",
       "\t 203    & 346    & 111.32 & 0     \\\\\n",
       "\t 203    & 437    & 109.43 & 0     \\\\\n",
       "\t 203    & 628    & 103.77 & 0     \\\\\n",
       "\t 203    & 724    &  75.47 & 1     \\\\\n",
       "\t 203    & 752    & 103.77 & 0     \\\\\n",
       "\t 203    & 899    &  90.57 & 0     \\\\\n",
       "\t 203    & 920    &  94.34 & 0     \\\\\n",
       "\t 204    & 203    &  94.34 & 0     \\\\\n",
       "\t 204    & 266    & 103.77 & 0     \\\\\n",
       "\t 204    & 321    &  94.34 & 0     \\\\\n",
       "\t 204    & 355    &  98.11 & 0     \\\\\n",
       "\t 204    & 441    &  94.34 & 0     \\\\\n",
       "\t 204    & 532    &  92.45 & 0     \\\\\n",
       "\t 204    & 621    &  86.79 & 0     \\\\\n",
       "\t 204    & 741    &  92.45 & 0     \\\\\n",
       "\t 204    & 812    & 101.89 & 0     \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| ptnum | days | fev | acute |\n",
       "|---|---|---|---|\n",
       "| 1      |  190   |  95.16 | 0      |\n",
       "| 1      |  249   | 112.90 | 0      |\n",
       "| 1      |  278   | 101.61 | 0      |\n",
       "| 1      |  369   | 109.68 | 0      |\n",
       "| 1      |  466   | 112.90 | 0      |\n",
       "| 1      |  557   | 114.52 | 0      |\n",
       "| 1      |  606   |  91.94 | 1      |\n",
       "| 1      |  613   |  83.87 | 1      |\n",
       "| 1      |  617   | 101.61 | 1      |\n",
       "| 1      |  634   |  85.48 | 0      |\n",
       "| 1      |  637   |  85.48 | 0      |\n",
       "| 1      |  666   |  82.26 | 1      |\n",
       "| 1      |  669   |  61.29 | 1      |\n",
       "| 1      |  670   |  80.65 | 1      |\n",
       "| 1      |  672   |  67.74 | 1      |\n",
       "| 1      |  678   |  95.16 | 1      |\n",
       "| 1      |  693   |  87.10 | 0      |\n",
       "| 1      |  812   | 101.61 | 0      |\n",
       "| 1      |  902   |  96.77 | 0      |\n",
       "| 1      |  932   | 101.61 | 0      |\n",
       "| 1      |  991   | 101.61 | 0      |\n",
       "| 1      | 1098   | 101.61 | 0      |\n",
       "| 1      | 1139   |  93.55 | 0      |\n",
       "| 1      | 1152   | 101.61 | 0      |\n",
       "| 1      | 1188   |  93.55 | 0      |\n",
       "| 1      | 1280   |  93.55 | 0      |\n",
       "| 1      | 1342   |  82.26 | 1      |\n",
       "| 1      | 1345   |  93.55 | 1      |\n",
       "| 1      | 1348   |  90.32 | 0      |\n",
       "| 1      | 1351   |  91.94 | 0      |\n",
       "| ... | ... | ... | ... |\n",
       "| 201    | 370    | 109.80 | 0      |\n",
       "| 201    | 404    | 109.80 | 0      |\n",
       "| 201    | 474    | 105.88 | 0      |\n",
       "| 201    | 565    | 117.65 | 0      |\n",
       "| 201    | 670    | 117.65 | 0      |\n",
       "| 201    | 754    | 117.65 | 0      |\n",
       "| 201    | 922    | 121.57 | 0      |\n",
       "| 202    | 245    | 118.87 | 0      |\n",
       "| 202    | 302    | 122.64 | 0      |\n",
       "| 202    | 392    | 128.30 | 0      |\n",
       "| 202    | 471    | 120.75 | 0      |\n",
       "| 202    | 595    | 105.66 | 0      |\n",
       "| 203    | 222    | 109.43 | 0      |\n",
       "| 203    | 283    | 107.55 | 0      |\n",
       "| 203    | 346    | 111.32 | 0      |\n",
       "| 203    | 437    | 109.43 | 0      |\n",
       "| 203    | 628    | 103.77 | 0      |\n",
       "| 203    | 724    |  75.47 | 1      |\n",
       "| 203    | 752    | 103.77 | 0      |\n",
       "| 203    | 899    |  90.57 | 0      |\n",
       "| 203    | 920    |  94.34 | 0      |\n",
       "| 204    | 203    |  94.34 | 0      |\n",
       "| 204    | 266    | 103.77 | 0      |\n",
       "| 204    | 321    |  94.34 | 0      |\n",
       "| 204    | 355    |  98.11 | 0      |\n",
       "| 204    | 441    |  94.34 | 0      |\n",
       "| 204    | 532    |  92.45 | 0      |\n",
       "| 204    | 621    |  86.79 | 0      |\n",
       "| 204    | 741    |  92.45 | 0      |\n",
       "| 204    | 812    | 101.89 | 0      |\n",
       "\n"
      ],
      "text/plain": [
       "     ptnum days fev    acute\n",
       "1    1      190  95.16 0    \n",
       "2    1      249 112.90 0    \n",
       "3    1      278 101.61 0    \n",
       "4    1      369 109.68 0    \n",
       "5    1      466 112.90 0    \n",
       "6    1      557 114.52 0    \n",
       "7    1      606  91.94 1    \n",
       "8    1      613  83.87 1    \n",
       "9    1      617 101.61 1    \n",
       "10   1      634  85.48 0    \n",
       "11   1      637  85.48 0    \n",
       "12   1      666  82.26 1    \n",
       "13   1      669  61.29 1    \n",
       "14   1      670  80.65 1    \n",
       "15   1      672  67.74 1    \n",
       "16   1      678  95.16 1    \n",
       "17   1      693  87.10 0    \n",
       "18   1      812 101.61 0    \n",
       "19   1      902  96.77 0    \n",
       "20   1      932 101.61 0    \n",
       "21   1      991 101.61 0    \n",
       "22   1     1098 101.61 0    \n",
       "23   1     1139  93.55 0    \n",
       "24   1     1152 101.61 0    \n",
       "25   1     1188  93.55 0    \n",
       "26   1     1280  93.55 0    \n",
       "27   1     1342  82.26 1    \n",
       "28   1     1345  93.55 1    \n",
       "29   1     1348  90.32 0    \n",
       "30   1     1351  91.94 0    \n",
       "...  ...   ...  ...    ...  \n",
       "5867 201   370  109.80 0    \n",
       "5868 201   404  109.80 0    \n",
       "5869 201   474  105.88 0    \n",
       "5870 201   565  117.65 0    \n",
       "5871 201   670  117.65 0    \n",
       "5872 201   754  117.65 0    \n",
       "5873 201   922  121.57 0    \n",
       "5874 202   245  118.87 0    \n",
       "5875 202   302  122.64 0    \n",
       "5876 202   392  128.30 0    \n",
       "5877 202   471  120.75 0    \n",
       "5878 202   595  105.66 0    \n",
       "5879 203   222  109.43 0    \n",
       "5880 203   283  107.55 0    \n",
       "5881 203   346  111.32 0    \n",
       "5882 203   437  109.43 0    \n",
       "5883 203   628  103.77 0    \n",
       "5884 203   724   75.47 1    \n",
       "5885 203   752  103.77 0    \n",
       "5886 203   899   90.57 0    \n",
       "5887 203   920   94.34 0    \n",
       "5888 204   203   94.34 0    \n",
       "5889 204   266  103.77 0    \n",
       "5890 204   321   94.34 0    \n",
       "5891 204   355   98.11 0    \n",
       "5892 204   441   94.34 0    \n",
       "5893 204   532   92.45 0    \n",
       "5894 204   621   86.79 0    \n",
       "5895 204   741   92.45 0    \n",
       "5896 204   812  101.89 0    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "07ed41be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "msm(formula = fev ~ days, subject = ptnum, data = fev, qmatrix = three.q,     hmodel = hmodel1, hcovariates = list(~acute, ~acute, NULL),     hconstraint = list(acute = c(1, 1)), death = 3, method = \"BFGS\")\n",
       "\n",
       "Maximum likelihood estimates\n",
       "Baselines are with covariates set to their means\n",
       "\n",
       "Transition intensities\n",
       "                  Baseline                          \n",
       "State 1 - State 1 -7.038e-04 (-8.333e-04,-0.0005945)\n",
       "State 1 - State 2  6.275e-04 ( 5.201e-04, 0.0007572)\n",
       "State 1 - State 3  7.631e-05 ( 3.967e-05, 0.0001468)\n",
       "State 2 - State 2 -8.010e-04 (-1.013e-03,-0.0006336)\n",
       "State 2 - State 3  8.010e-04 ( 6.336e-04, 0.0010126)\n",
       "\n",
       "Hidden Markov model, 3 states\n",
       "State 1 - normal distribution\n",
       "Parameters: \n",
       "       Estimate       LCL       UCL\n",
       "mean  97.933879 97.272340 98.595419\n",
       "sd    16.185028 15.777826 16.602739\n",
       "acute -8.791774 -9.951421 -7.632128\n",
       "\n",
       "State 2 - normal distribution\n",
       "Parameters: \n",
       "       Estimate       LCL       UCL\n",
       "mean  51.752521 50.693089 52.811954\n",
       "sd    17.676279 17.082764 18.290415\n",
       "acute -8.791774 -9.951421 -7.632128\n",
       "\n",
       "State 3 - identity distribution\n",
       "Parameters: \n",
       "      Estimate LCL UCL\n",
       "which      999  NA  NA\n",
       "\n",
       "\n",
       "-2 * log-likelihood:  51597.89 \n",
       "[Note, to obtain old print format, use \"printold.msm\"]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "three.q <- rbind(c(0, exp(-6), exp(-9)), c(0, 0, exp(-6)), c(0, 0, 0))\n",
    "hmodel1 <- list(hmmNorm(mean = 100, sd = 16), hmmNorm(mean = 54, sd = 18), hmmIdent(999))\n",
    "fev1.msm <- msm(fev ~ days, subject = ptnum, data = fev, qmatrix = three.q, hmodel = hmodel1, \n",
    "               hcovariates = list(~ acute, ~ acute, NULL), hconstraint = list(acute = c(1, 1)), death = 3, \n",
    "               method = \"BFGS\")\n",
    "fev1.msm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04dbaa76",
   "metadata": {},
   "source": [
    "- In State 1 (No BOS), FEV mean is 98% of baseline (plus or minus 16%).\n",
    "- In State 2 (BOS), FEV mean is 52% of baseline (plus or minus 18%).\n",
    "- FEV is estimated to be 9% lower within 14 days of acute illness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "50c6ee89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>estimates</th><th scope=col>SE</th><th scope=col>L</th><th scope=col>U</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>State 1</th><td>1420.777 </td><td>122.3945 </td><td>1200.0469</td><td>1682.107 </td></tr>\n",
       "\t<tr><th scope=row>State 2</th><td>1248.390 </td><td>149.3067 </td><td> 987.5232</td><td>1578.169 </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       "  & estimates & SE & L & U\\\\\n",
       "\\hline\n",
       "\tState 1 & 1420.777  & 122.3945  & 1200.0469 & 1682.107 \\\\\n",
       "\tState 2 & 1248.390  & 149.3067  &  987.5232 & 1578.169 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | estimates | SE | L | U |\n",
       "|---|---|---|---|---|\n",
       "| State 1 | 1420.777  | 122.3945  | 1200.0469 | 1682.107  |\n",
       "| State 2 | 1248.390  | 149.3067  |  987.5232 | 1578.169  |\n",
       "\n"
      ],
      "text/plain": [
       "        estimates SE       L         U       \n",
       "State 1 1420.777  122.3945 1200.0469 1682.107\n",
       "State 2 1248.390  149.3067  987.5232 1578.169"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sojourn.msm(fev1.msm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28197488",
   "metadata": {},
   "source": [
    "The average estimate and confidence intervals for -1/qrr - the average onset and progression in days:\n",
    "\n",
    "BOS State 1: 1,421 days or about 3 years after transplantation\n",
    "\n",
    "BOS State 2: 1,248 days or about an additional 3 years later\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e779798b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating the most likely true series of states underlying the data using the Viterbi algorithm\n",
    "\n",
    "#\n",
    "# viterbi.msm(fev1.msm)\n",
    "# Doesn't work yet:\n",
    "# https://www.rdocumentation.org/packages/msm/versions/1.6.8/topics/viterbi.msm"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}