:
Machine Learning for Health (ML4H) 2020

# Identifying Decision Points for Safe and Interpretable Reinforcement
Learning in Hypotension Treatment

Kristine Zhang∗<EMAIL_ADDRESS>
Yuanheng Wang∗<EMAIL_ADDRESS>
Jianzhun Du∗<EMAIL_ADDRESS>
Brian Chu<EMAIL_ADDRESS>
Harvard University  Leo Anthony Celi<EMAIL_ADDRESS>
Massachusetts Institute of Technology  Ryan Kindle<EMAIL_ADDRESS>
Massachusetts General Hospital  Finale Doshi-Velez<EMAIL_ADDRESS>
Harvard University

###### Abstract

Many batch RL health applications first discretize time into fixed intervals.
However, this discretization both loses resolution and forces a policy
computation at each (potentially fine) interval. In this work, we develop a
novel framework to compress continuous trajectories into a few, interpretable
decision points—places where the batch data support multiple alternatives. We
apply our approach to create recommendations from a cohort of hypotensive
patients dataset. Our reduced state space results in faster planning and
allows easy inspection by a clinical expert.

###### keywords:

batch RL

##  1 Introduction

Many works have investigated reinforcement learning for sequential decision-
making in health (Yu et al., 2019). A standard approach is to take already-
collected (batch) data and use that to train a policy. However, batch RL has
many subtleties, and the high stakes of healthcare settings necessitate
special attention to safety, interpretability, and robust off-policy
evaluation before polices gleaned from batch data are applied to patients
(Gottesman et al., 2019).

In this paper, we focus on the issue of time discretization in this batch RL
context. To date, most batch RL applications first discretized time into bins
and then solved for a policy decision for each bin. Not only are learnt
policies often sensitive to the binning method, but decisions must be solved
for every time step, affecting both computation and policy evaluation (Tallec
et al., 2019).

We ask and answer the question: what if we didn’t try to solve at every time
point? Indeed, especially with batch data, the _only_ times when one might be
able to give a recommendation are in states where clinicians treat similar
patients differently. We call these _decision points_ (DP). Compressing
trajectories down to high-ambiguity _decision regions_ , or areas of the state
space with many DPs, has advantages not only for efficient planning but also
for interpretability. Summarizing a continuous time and space problem in 20-40
discrete clusters enables clinicians to manually inspect the recommended
action for each cluster and thus globally validate a policy. Our contributions
are: (1) developing a general pipeline for learning this decision region MDP
and (2) demonstrating its interpretability on a hypotension management task.

##  2 Related Work

Many works have considered RL for health, including for hypotension management
in the ICU (Komorowski et al., 2018; Futoma et al., 2020; Srinivasan and
Doshi-Velez, 2020). Unlike our work, all of these approaches choose some time
discretization with a constant time step size and optimize decisions at each
of these points. Gottesman et al. (2020) propose an interpretable RL approach
to handling treatment changes or patient measurements with irregular
observations, which is similar to our methods, though we leverage the
uncertainty of the behavior policy for both state and time compression. Du et
al. (2020) take a model-based approach to continuous-time RL with
demonstrations in the healthcare environment, whereas it cannot be applied to
the batch RL setting.

More broadly, state abstraction (e.g. Li et al. (2006)) is a popular RL
problem. However, these works have been largely in online settings and focused
on aggregation based on policy-relevant information like $Q$-values or
transitions; they also only compress the state space and the trajectory
lengths remain the same. In the context of batch RL, we instead aim to
identify important states for decision-making and filter out regions where the
experts have consensus.

In this way, our methods enable efficient planning and policy validation.
Komorowski et al. (2018) also employ clustering to discretize a continuous
environment; however, their state clustering methods do not consider practice
variation and thus produce large state spaces (750 vs. 32 for us). The issue
of time discretization also remains. Other works use deep models for state
representations and policies, but these are difficult to interpret, while ours
are easy to inspect.

##  3 Background and Notation

A Markov Decision Process (MDP) is a tuple
$\langle\mathcal{X},\mathcal{A},\mathcal{T},\mathcal{R},\gamma\rangle$. We
assume a continuous state space $\mathcal{X}\in\mathbb{R}^{d}$, discrete
action space $\mathcal{A}$, and discount factor $\gamma\in(0,1]$.
$T(x^{\prime}|x,a)$ and $R(x,a)$ denote the state transition function and
reward function, respectively. A policy
$\pi:(\mathcal{X},\mathcal{A})\to[0,1]$ defines the probability of an action
given a state. The objective is to learn a policy that maximizes the expected
return, or sum of discounted rewards
$V^{\pi}=E\left[\sum_{t}\gamma^{t}r_{t}|a_{t}\sim\pi\right]$.

In batch reinforcement learning, we cannot interact with the environment
during learning, but instead begin with a pre-collected dataset of
trajectories $\mathcal{D}=\\{(x_{i},a_{i},x_{i}^{\prime},r_{i})\\}_{i=1}^{N}$.
These trajectories are sampled by following a (possibly unknown) behavior
policy $\pi_{b}$ and used to learn a new evaluation policy $\pi_{e}$. MDPs
with small discrete state and action spaces can be solved using value or
policy iteration (Sutton and Barto, 2018).

##  4 Method

The innovation of our method is that we do not consider every state to be a
potential decision point. Rather, we only consider those points with high
behavior policy variability, that is, states whose similar neighbors are
administered different treatments. While this idea is simple, its realization
requires formalizing many subtleties, which we describe below.

Figure 1: Compression process from a trajectory in the original state space to
decision clusters. White squares are non-DP states, colors denote DPs in the
same cluster. The mortality outcome is captured via two more states: A (alive)
or D (dead).

#### Decision Point Identification

We operationalize the notion of states with high clinical variation by first
training a kernel-based classifier to _predict_ clinician actions. States with
sufficient neighbors (according to the kernel) whose actions still cannot be
predicted well are then considered _decision points_.

Specifically, we consider a weighted Gaussian kernel:

$k(x,x^{\prime})=\exp(-\|(x-x^{\prime})\|_{2}^{2})$ (1)

where $k(x,x^{\prime})$ is the estimated similarity between states, and
$\bm{w}$ can be interpreted as an importance weighting over state dimensions
and $\odot$ represents an element-wise multiplication between vectors. We
learn the parameters $\bm{w}$ by backpropagating through the kernel space with
optimizing the cross-entropy loss of predicting behavior actions. Details can
be found in Appendix B.1.

The second step is to determine a set of neighbors $\\{x^{\prime}\\}$ for each
state $x$ in the dataset. We introduce two hyperparameters: the minimum kernel
similarity threshold $\delta\in(0,1]$ to be considered a neighbor, and the
minimum number of neighbors $n$ for an action to be considered valid. That is,
if a state $x$ has at least $n$ neighbors within $\delta$ performing action
$a$, then action $a$ is considered a valid action for $x$. If a state has more
than one valid action, it is considered a decision point.

#### Clustering Decision Points

The process above identifies whether each state in the dataset is a decision
point. Next, we cluster only the decision points to identify discrete
_decision regions_ that will form our reduced MDP. A top-down hierarchical
clustering approach is applied, where the Distances between decision points
are calculated using the standardized state features (dimensions of $x$). For
each intermediate cluster, we examine the difference of mean values for each
feature across actions. If the difference of these means exceed specified
thresholds, this implies the region is not homogeneous and we further split
clusters (we truly want decision regions where the optimal action is unclear).
Secondly, we examine if the current clustering is creating loops—defined as
trajectories that leave a cluster and coming back to the same cluster within
the next three time steps—as these (artificial) loops can make the agent
believe that it can “freeze time” and avoid any final consequence. If loop
percentage is high, we split that cluster.

Figure 2: Three policies obtained using three different rewards functions, in
comparison to the behavior policy.

Policy | WIS Score | ESS
---|---|---
Behavior | -1.11 | 3,884
Piecewise | -0.93 | 1,945
Mortality | -0.94 | 1,895
Terminal | -0.94 | 1,921

Table 1: OPE results from different policies. “ESS” means effective sample
size.

#### Getting a Policy

Finally, we are ready to create the MDP and solve for a policy. We convert
each patient trajectory of the form $\\{(x,a,x^{\prime})\\}$ into a new,
compressed trajectory $\\{(\bar{x},\bar{a},\bar{x}^{\prime})\\}$. We use these
transitions to learn a discrete, maximum-likelihood MDP in which only valid
actions are allowed, and then solve.

The compression process from $\\{(x,a,x^{\prime})\\}$ to
$\\{(\bar{x},\bar{a},\bar{x}^{\prime})\\}$ is summarized by Figure 1 and
formalized in Appendix B.4. At a high level, we append a new cluster
$\bar{x}_{t}$ each time a transition is observed from a different cluster or
non-DP state. Condensing contiguous DPs in the same cluster and ignoring non-
DPs is based on the idea that a “decision” is an action taken in one decision
region that causes movement to another decision region. In the final tuple of
a trajectory, we include a transition to the mortality state of the patient to
help track their outcome.

##  5 Case Study for Hypotension Management

#### Setting

Now, we apply our method to recommend treatments to manage hypotension in the
ICU. We select a cohort of 11,739 patients records from the MIMIC-III dataset
(See Appendix A for details). Our goal is to (1) recommend a set of treatment
policies that are reasonable, and (2) offer a toolkit with visualizations for
clinicians to validate and refer to when making a treatment decision for
hypotensive patients. Since this is a proof-of-concept, we only consider four
action classes: do nothing, use fluids, use vasopressors, and use both.

To check for sensitivity to reward formulation, we solved with three reward
functions: (1) Piecewise: based on the MAP and urine volume for each patient
at each timestamp. Lower map and urine values correspond to lower rewards, (2)
Mortality probability: inverse of the mortality probability of each cluster-
treatment pair in the constructed MDP, and (3) Terminal state: Based on
whether patients reach Alive ($R=1$) or Death ($R=0$) at the end of the ICU
admission. We trained the policy using value iteration and use a discount rate
of 0.98.

#### Results

Running our algorithm on these data results in 33 clusters, where the largest
cluster includes 117,730 decision points and the smallest cluster includes 13
decision points. Figure 4 in Appendix B.3 shows the mean value of each feature
for clusters containing high mortality rate ($\geq 10\%$) and enough data
($\geq 10$ patients treatment points for at least two types of treatments);
this figure was easily inspected by intensivists.

Figure 2 shows the optimal policy for each reward function. We also include
the behavior policy for comparison, which is the empirical proportion of each
action taken in the historical treatment data. The fact that each policy from
different rewards functions is slightly different but largely similar to
behavior policy demonstrates the robustness of our methods, and the fact it
can be so succinctly summarized allowed us to discuss each case with
clinicians; in each case, they could posit whether the policy recommendation
made sense or could have been due to some confound. Such discussion would not
be possible with a black-box policy.

To quantitatively evaluate the evaluation policy learned using the
aforementioned rewards, we use _weighted importance sampling_ (WIS) for off
policy evaluation (OPE), whose estimate is biased but consistent
(asymptotically correct) with much lower variance. We apply standard weight
clipping by capping weights at 95 percentile. Table 1 reports the values of
WIS estimates. It shows all three evaluation policies achieve similar rewards,
which is higher than the reward of the behavior policy.

##  6 Conclusion

In this work, we present a pipeline for identifying the decision regions in an
RL environment and employing them to construct a summary MDP. This summary MDP
can then be used for efficient planning and policy validation by human
experts. We demonstrate the value of this new approach to interpretability in
healthcare settings using the example of hypotension data. Our proposed
framework is highly generalizable and can be fine-tuned for other applications
in healthcare and beyond.

## References

  * Du et al. (2020) Jianzhun Du, Joseph Futoma, and Finale Doshi-Velez.  Model-based reinforcement learning for semi-markov decision processes with neural odes.  _arXiv preprint arXiv:2006.16210_ , 2020.
  * Futoma et al. (2020) Joseph Futoma, Michael C Hughes, and Finale Doshi-Velez.  Popcorn: Partially observed prediction constrained reinforcement learning.  _arXiv preprint arXiv:2001.04032_ , 2020.
  * Gottesman et al. (2019) Omer Gottesman, Fredrik Johansson, Matthieu Komorowski, Aldo Faisal, David Sontag, Finale Doshi-Velez, and Leo Anthony Celi.  Guidelines for reinforcement learning in healthcare.  _Nat Med_ , 25(1):16–18, 2019.
  * Gottesman et al. (2020) Omer Gottesman, Joseph Futoma, Yao Liu, Sonali Parbhoo, Leo Anthony Celi, Emma Brunskill, and Finale Doshi-Velez.  Interpretable off-policy evaluation in reinforcement learning by highlighting influential transitions, 2020.
  * Johnson et al. (2016) Alistair EW Johnson, Tom J Pollard, Lu Shen, H Lehman Li-wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark.  Mimic-iii, a freely accessible critical care database.  _Scientific data_ , 3:160035, 2016.
  * Komorowski et al. (2018) Matthieu Komorowski, Leo A. Celi, Omar Badawi, Anthony C. Gordon, and A. Aldo Faisal.  The artificial intelligence clinician learns optimal treatment strategies for sepsis in intensive care.  _Nature Medicine_ , 24(11):1716–1720, 2018.  10.1038/s41591-018-0213-5.  URL https://doi.org/10.1038/s41591-018-0213-5.
  * Li et al. (2006) Lihong Li, Thomas J. Walsh, and Michael L. Littman.  Towards a unified theory of state abstraction for mdps.  In _In Proceedings of the Ninth International Symposium on Artificial Intelligence and Mathematics_ , pages 531–539, 2006.
  * Nguyen et al. (2017) Tu Dinh Nguyen, Trung Le, Hung Bui, and Dinh Phung.  Large-scale online kernel learning with random feature reparameterization.  In _Proceedings of the 26th International Joint Conference on Artificial Intelligence_ , IJCAI’17, page 2543–2549. AAAI Press, 2017.  ISBN 9780999241103.
  * Rahimi and Recht (2008) Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  In _Advances in neural information processing systems_ , pages 1177–1184, 2008.
  * Srinivasan and Doshi-Velez (2020) Srivatsan Srinivasan and Finale Doshi-Velez.  Interpretable batch irl to extract clinician goals in icu hypotension management.  _AMIA Jt Summits Transl Sci Proc_ , pages 636––645, 2020.
  * Sutton and Barto (2018) Richard S. Sutton and Andrew G. Barto.  _Reinforcement Learning: An Introduction_.  A Bradford Book, Cambridge, MA, USA, 2018.  ISBN 0262039249.
  * Tallec et al. (2019) Corentin Tallec, Léonard Blier, and Yann Ollivier.  Making deep q-learning methods robust to time discretization, 2019.
  * Yu et al. (2019) Chao Yu, Jiming Liu, and Shamim Nemati.  Reinforcement learning in healthcare: A survey.  _arXiv preprint arXiv:1908.08796_ , 2019.

##  Appendix A Data Processing

###  A.1 MIMIC Data Cleaning

Our hypotension data is obtained from the publicly available MIMIC-III dataset
for ICU patients (Johnson et al., 2016), with records from 15,653 unique ICU
stays. Patients are considered hypotensive if their systolic blood pressure
(SBP) drops below 90 mmHg.

Preprocessing was performed and is described in greater detail by Futoma et
al. (2020). Trajectories were initially discretized into hourly bins starting
at $t=1$ hour into ICU admission and were truncated either at discharge or at
$t=72$ hours (as patients staying beyond 72 hours tend to be abnormal cases).

The initial state space consists of 111 clinical features including patient
measurements, vital signs, and records of past treatment actions. These
dimensions contain some redundancy with indicators and quantitative measure
for the same variables.

For actions, we track the amount of either fluids or vasopressors administered
to each patient at each time step. For this pipeline, we primarily consider
the four general actions corresponding to “no treatment”, “fluids only”,
“vasopressors only”, or “fluids and vasopressors both given”. The learned
policy is intended to help physicians decide whether to give each type of
treatment at all, as there are existing guidelines for the quantity to give
once a choice has been made. Actions are unbalanced in this dataset, with
around 80% of the states receiving no treatment. To analyze treatment impact
on outcome, we define a patient’s mortality to be true if they died in the
hospital within 30 days of ICU admission.

###  A.2 Dataset Split

Because we lack access to a simulator for MIMIC, the efficacy of our methods
must be tested on the batch data. We split the trajectories in the original
dataset 75-25 into train and test sets, yielding 555,196 transition tuples
from 11,739 trajectories in the train set and 185,671 from 3,914 trajectories
in the test set.

The train set is used for initial data analysis, MDP compression, and
planning. The test set is reserved for OPE.

###  A.3 State Feature Selection

Due to the high dimensionality of the starting MIMIC state space, we
preemptively select for features that are most relevant to action prediction.
The objective was to increase the efficiency of kernel-learning and avoid
collinearity between redundant features. We construct a random forest
predictor from states to actions, implemented using a random forest classifier
with 100 estimators, depth 3, and balanced class weights. We employ only the
top 20 features ranked by random forest feature importance in the kernel
learning algorithm. Figure 4 lists all continuous features within this top 20
and the features clinicians consider to be important.

##  Appendix B Implementation Details

###  B.1 Kernel Learning

We are interested in evaluating the Gaussian kernel distance between states,
defined as

$k(x,x^{\prime})=\exp(-\|(x-x^{\prime})\|_{2}^{2})$ (2)

Naively computing the pairwise distance between $N$ states of dimension $d$ is
computationally expensive, given $\mathcal{O}(N^{2}d)$ complexity and the
large MIMIC-III dataset. Rahimi and Recht (2008) showed that this process can
be accelerated through the use of _Random Fourier Features_ (RFF). After
applying a randomized feature map $z:\mathbb{R}^{d}\to\mathbb{R}^{D}$ from
state vectors to a low-dimensional Euclidean inner product space, kernel
evaluation can be estimated by the inner product of the randomized features:

$k(x,x^{\prime})\approx z(x)^{\top}z(x^{\prime})$ (3)

We can therefore approximately train kernel machines by projecting to the
Fourier space and applying fast linear methods such as regression. We evaluate
the efficacy of a set of transformations $\bm{w}$ by performing multi-class
logistic regression to predict actions. We define $z(x_{i})=RFF(\bm{w}\odot
x_{i})$ to be the covariates such that $k(x,x^{\prime})\approx
z(x)^{\top}z(x^{\prime})$. The targets are $a_{i}$ with regression
coefficients $\bm{V}$. The objective function is the corresponding cross-
entropy loss in Equation 4.

$\begin{gathered}\min_{\bm{w},\bm{V}}-\sum_{i=1}^{N}\sum_{a\in\mathcal{A}}\mathbb{I}(a_{i}=a)\hat{p}(a|x_{i})\\\
\text{with}\quad\hat{p}(a|x)=\text{softmax}(z(x)^{\top}\bm{V})_{a}\end{gathered}$
(4)

We propose Algorithm 1 that uses gradient descent to simultaneously optimize
$\bm{w}$ and $\bm{V}$. Nguyen et al. (2017) showed that the gradient of the
cross-entropy loss can be calculated over $\bm{w}$ despite the random feature
mapping. We also apply a reparameterization trick, optimizing over
$\bm{u}=\log\bm{w}$ to ensure that $\bm{w}$ is positive.

Input : Set of state-action tuples $\\{x_{i},a_{i}\\}_{i=1}^{N}$; learning
rate $\eta$, number of epochs $M$, mini-batch size $B$.

Output : Learned kernel weights $\bm{w}$.

Initialize
$\bm{u}\leftarrow\bm{0}^{d},\bm{V}\leftarrow\bm{1}^{D\times|\mathcal{A}|}$;

Generate randomized parameters $\omega\in\mathbb{R}^{d},b\in\mathbb{R}$
according to the RFF algorithm;

for  _$l\leftarrow 1$ to $M$_ do

for  _$j\leftarrow 1$ to $\lfloor N/B\rfloor$_ do

Select mini-batch of states $X_{j}$ and actions $A_{j}$;

$\bm{w}\leftarrow\exp(\bm{u})$

Project $z(x)\leftarrow\text{RFF}((\bm{w}\odot x;\omega,b)$ for all $x\in
X_{j}$;

$L\leftarrow$ cross-entropy loss when predicting $A_{j}$ with features $Z_{j}$
and weights $\bm{V}$;

$\bm{V}\leftarrow\bm{V}-\eta\frac{\partial L}{\partial\bm{V}}$;

$\bm{u}\leftarrow\bm{u}-\eta\frac{\partial L}{\partial\bm{u}}$;

end for

end for

Return $\bm{w}=\exp(\bm{u})$;

Algorithm 1 Learning Weighted Gaussian Kernel Figure 3: $\epsilon$
conceptually represents distance between points (which inversely approximate
similarity $\delta$). Higher similarity threshold$\delta$ translates to a
smaller $\epsilon$-ball, which is equivalent to a smaller neighborhood.

###  B.2 DP Hyperparameter Choice

To identify decision points, we have two main hyperparameters: $\delta$, the
minimum similarity from a point $x$ to neighbors, and $n$ the number of
neighbors with the same action (Figure 3). Corresponding actions for points
present in this $\epsilon$-ball are considered to be allowable actions, and
decision points are defined as points that exist in this $\epsilon$-ball where
neighbors do not dominantly take a sole action. Since we are more interested
to measure the similarity between every point in the kernel space, we instead
consider points within the similarity threshold $\delta$. Also, we require a
decision point must have a minimum number of neighbors $n$ in this
$\delta$-specified region. We perform the grid search over a set of candidate
values on a sampled validation dataset. For each pair of $\delta,n$, we
evaluate the overall AUC score for a sampled subset of points, whose labels
are determined by the empirical distribution of labels present in the region.
Conducting such experiment multiple rounds, we finally choose the similarity
threshold $\delta=0.95$ and $n=20$ that achieves the highest AUC score.

Input : Trajectory $\\{(x_{i},a_{i})\\}_{i=1}^{L}$ with corresponding cluster
labels $\\{c_{i}\\}_{i=1}^{L}$ (with $c_{i}=0$ for non-DPs); mortality state
$I_{m}$, action summary function $h$.

Output : Lists of abstracted states and actions
$\\{\bar{x}_{j}\\},\\{\bar{a}_{j}\\}$.

Initialize $\bar{X}\leftarrow[],\bar{A}\leftarrow[],c_{0}\leftarrow
0,A\leftarrow[]$ ;

for  _$l\leftarrow 1$ to $L$_ do

if  _$c_{l} >0$_ then

$A$.append($a_{l}$) ;

_// entering new decision region_

if  _$c_{l}\neq c_{0}$_ then

$\bar{X}$.append($c_{l}$);

end if

 _// leaving previous decision region_

if  _$l=L$ or $c_{l+1}\neq c_{l}$_ then

$\bar{A}$.append($h(A)$) ;

$A\leftarrow[]$;

end if

end if

end for

$\bar{X}$.append($I_{m}$);

Return $\bar{X},\bar{A}$;

Algorithm 2 MDP Trajectory Compression Figure 4: The feature means for points
in each decision region. Sorted by descending order of Mean Arterial Pressure
(MAP) value.

###  B.3 Decision Region Clustering

In order to form decision regions efficiently, we utilize the hierarchical
clustering implementations in scipy. We take the output from scipy linkage
function call, which is a matrix specifying the order of points being merged
at each step. We iterate through the matrix in a top-down manner, and
determine if the children of each split meets the specified condition, i.e.,
the differences of the means of feature values for points in each child are
smaller than the predetermined threshold. If not, we continue iterate through
the matrix until conditions are satisfied or maximum number of splits are
reached.

###  B.4 MDP Compression Details

Algorithm 2 formalizes the MDP compression method. Recall that we convert each
patient trajectory $\\{(x,a,x^{\prime})\\}$ into a compressed trajectory
$\\{(\bar{x},\bar{a},\bar{x}^{\prime})\\}$, where all
$\bar{x},\bar{x}^{\prime}\in\mathcal{C}$ are decision clusters and
$\bar{a}\in\bar{\mathcal{A}}$ are summarized actions. The summary function $h$
should compress multiple actions in $\mathcal{A}$ that take place while the
patient state remains in the same cluster.

Given a set of compressed trajectories, we estimate a behavior policy
$\bar{\pi}_{b}$ and transition function
$\bar{T}(\bar{x}^{\prime}|\bar{x},\bar{a})$ from the new dataset, as
formalized in Equation 5 and 6. Following standard state abstraction protocol
(Li et al., 2006), the behavior policy is estimated as the proportion of times
an action was taken from a given cluster, and the transition function is
estimated as the proportion of times a given cluster-action pair led to
another cluster.

$\bar{\pi}_{b}(\bar{a}^{*}|\bar{x}^{*})=\frac{\sum_{i}\mathbb{I}(\bar{x}_{i}=\bar{x}^{*},\
\bar{a}_{i}=\bar{a}^{*})}{\sum_{i}\mathbb{I}(\bar{x}_{i}=\bar{x}^{*})}$ (5)
$\bar{T}(\bar{x}^{\prime*}|\bar{x}^{*},\bar{a}^{*})=\frac{\sum_{i}\mathbb{I}(\bar{x}_{i}=\bar{x}^{*},\
\bar{a}_{i}=\bar{a}^{*},\
\bar{x}^{\prime}=\bar{x}^{\prime*})}{\sum_{i}\mathbb{I}(\bar{x}_{i}=\bar{x}^{*},\
\bar{a}_{i}=\bar{a}^{*})}$ (6)

This process can be thought of as simultaneous state and temporal abstraction.
We perform state abstraction to map decision points into high-uncertainty
clusters, then temporal abstraction to condense all actions that do not cause
transitions out of decision regions. In this way, we reduce both the size of
the state space and the length of trajectories.


# Statistical Modeling of Atmospheric Mean Temperature in sub$-$Sahel West
Africa

Falaiye O. A  Department of Physics, Faculty of Physical Sciences, University
of Ilorin, Nigeria  Sukam Y. M  Department of Physics, Faculty of Physical
Sciences, University of Ilorin, Nigeria  Abimbola O. J  Department of Physics,
Faculty of Science, Federal University Lafia, Nigeria

## I. Introduction

Preciptable water vapor (PWV) is an important atmospheric parameter essential
in both weather and climatic prediction. Knowledge of the variability of PWV
is also very important in astronomy as water vapor interact with the incoming
electromagnetic waves in the atmosphere.

GPS meteorology offers a real-time continuous measurement of PWV. From the GPS
data, it is possible to estimate the zenith total delay (ZTD). The zenith
hydrostatic delay (ZHD), which is the delay error caused by the dry component
of the atmosphere could be estimated using models such as the Saastamoinen [7]
model. The delay error introduced by the water vapor in the atmosphere, called
the zenith wet delay (ZWD) is calculated thus:

$ZWD=ZTD-ZHD$ (1)

The ZWD is usually transformed into the PWV using the transformation equation
as given below

$PWV=\Pi\times ZWD$ (2)

The transformation constant $\Pi$, is given by

$\Pi=\left[10^{-6}\left(\frac{k_{3}}{T_{m}}+k^{{}^{\prime}}_{2}\right)\rho_{v}R_{v}\right]^{-1}$
(3)

here, $\rho_{v}$ is the liquid water density = $1000\ kgm^{-3}$, $R_{v}$ is
water vapor gas constant = $461.524\ JK^{-1}kg^{-1}$, $k_{3}=377600\
K^{2}hPa$, $k^{{}^{\prime}}_{2}=22.1\ KhPa^{-1}$ and $T_{m}$ is called
weighted atmospheric mean temperature determination of which is very crucial
to accurate transformation of ZWD into the PWV.

The weighted mean atmospheric temperature $T_{m}$, could be defined as

$T_{m}=\frac{\int^{p_{1}}_{p_{2}}Td(lnp)}{\int^{p_{1}}_{p_{2}}d(lnp)}=\frac{\int^{\infty}_{z}\frac{e}{T}dz}{\int^{\infty}_{z}\frac{e}{T^{2}}dz}$
(4)

where $p_{1}$ and $p_{2}$ are pressures at two different layers of the
atmosphere. $T$ is the ambient temperature. $e$ is the vapor pressure while
$z$ is the height above the ground. From Eq. 4. $T_{m}$ could be estimated
from the radiosonde data taken at different layers of the atmosphere. However,
such radiosonde data are usually not readily available across the globe,
especially Africa, hence several attempts have been made to statistically
relate $T_{m}$ with the surface temperature $T_{s}$ which is readily available
as a basic meteorological parameter at any weather observing station. Bevis et
al.[1] used $8700$ radiosonde profiles from $13$ stations in the United States
of America (USA) to obtain a statistical model of the form $T_{m}=aT_{s}+b$.
The Bevis et al.[1] model was adopted in many studies (e.g., Raju et al.[2];
Fernandez et al.[3]; Musa et al.[4]; Abimbola et al.[5]) across the globe
despite the fact that it was obtained for the USA region. Mendes et al.[6] and
Solbrig[8], using radiosonde data, obtained similar linear relation for
Germany, Raju et al.[2] for Indian sub-continent, Shoji[9] for Japan and
Isioye et al.[10] for West Africa. Meanwhile, Schuler et al.[11], using a
numerical weather prediction data, obtained a model which is more global in
outlook relating $T_{m}$ to $T_{s}$ as shown in Eq. 5

$T_{m}=0.65T_{s}+86.9$ (5)

From satellite data: which is of higher resolution spatio-temporally than the
radiosonde data, it is the aim of this work to find a suitable statistical
model for the evaluation of atmospheric mean temperature $T_{m}$ in West
African region.

## II. Methodology

### i. Study Area

This study was conducted for locations in West African region as shown in
Figure 1. The selected locations in West Africa for this research work are:
Dakar ($14.76$ ${}^{o}N$ Latitude, $17.36$ ${}^{o}W$ Longitude), Conakry
($9.64\;^{o}N$ Latitude, $13.57\;^{o}W$ Longitude), Bamako ($12.63\;^{o}N$
Latitude, $8.00\;^{o}W$ Longitude), Abidjan ($5.36\;^{o}N$ Latitude,
$4.00\;^{o}W$ Longitude), Niamey ($13.51\;^{o}N$ Latitude, $12.12\;^{o}E$
Longitude), and Abuja ($90.07\;^{o}N$ Latitude, $7.39\;^{o}E$ Longitude).

Figure 1. Map of West Africa showing Locations of study area (Encyclopeadia
Britannica, $2018$).

For ease of analysis of the data sets and climatic consideration, the study
locations were further divided into ($1$) Hinterland, that is those locations
in the interior of the study region and ($2$) Coastal, that is those regions
close to the Atlantic Ocean. These divisions are shown in Table 1.

Table 1. Division of the Study Area into Hinterland and Coastal Regions. Station |  | Country | Longitude | Latitude | Period of data
---|---|---|---|---|---
Dakar | Coastal Region | Senegal | $17.366\;^{o}W$ | $14.765\;^{o}N$ | $2004-2016$
Conakry | Guinea | $13.578\;^{o}W$ | $09.641\;^{o}N$ | $2004-2016$
Abidjan | Cote d’ivoire | $04.008\;^{o}W$ | $05.360\;^{o}N$ | $2004-2016$
Bamako | Hinter- land Region | Mali | $08.003\;^{o}W$ | $12.639\;^{o}N$ | $2004-2016$
Niamey | Niger | $12.125\;^{o}E$ | $13.512\;^{o}N$ | $2004-2016$
Abuja | Nigeria | $07.399\;^{o}E$ | $09.077\;^{o}N$ | $2004-2016$

### ii. Data Analysis

The data used in this work was obtained from the archive of the Climate
Monitoring Satellite Application Facilities (CM-SAF) of the European
Meteorological Satellites (EUMETSAT) at https://www.cmsaf.eu. The data was
obtained as a monthly average for a time period spanning from $2004$ to $2016$
in Network Common Data Format (netCDF). The data obtained at six different
pressure levels (level$1$ = $200$ $hPa$; level$2$ = $300$ $hPa$; level$3$ =
$500$ $hPa$; level$4$ = $700$ $hPa$; level$5$ = $850$ $hPa$ and level$6$
(surface) = $1000$ $hPa$). are the specific humidity $q_{m}$ ($g/kg$) and the
ambient temperature $T$ ($K$). The vapor pressure $e$ (in $hPa$) was estimated
from the specific humidity using:

$e=\frac{0.5q_{m}}{10^{-3}q_{m}+0.622}$ (6)

Panoply (v$4.3.1$), open source software from NOAA was used for the processing
of the data while Jupyter notebook, an environment for running Python code,
was used for data plotting, curve fittings and other statistical analysis.

Goodness of fit and correlational analysis of the derived empirical equations
were done using coefficient of determination ($R^{2}$), mean-biased error
($MBE$) and root-mean-square error ($RMSE$).

## III. Results and Discussion

The derived empirical equations relating the mean atmospheric temperature
$T_{m}$ to the surface temperature $T_{s}$ for each of the locations
considered in this work within the West African region are shown in Table 2.
All the stations considered show good linear correlations between $T_{m}$ and
$T_{s}$Except for Conakry with a below average coefficient of determination.
It will be observed from Table 2 that the coastal region generally does not
show good linear correlation between $T_{m}$ and $T_{s}$, whereas the
hinterland region shows very good linear correlation: also it was observed
that the ambient surface temperature at the coast is generally much lower than
that at the hinterland of West Africa.

Table 2. Statistics for each of the stations in West Africa. Location |  |  Empirical Model $T_{m}=aT_{s}+b$ |  Coefficient of Determination
---|---|---|---
Name |  | $a$ | $b$ | $R^{2}$
Dakar | Coastal Region | $0.57$ | $123.15$ | $0.575$
Conakry | $1.19$ | $63.28$ | $0.442$
Abidjan | $0.87$ | $32.68$ | $0.741$
Bamako | Hinter- land Region | $0.68$ | $89.96$ | $0.646$
Niamey | $0.81$ | $52.23$ | $0.763$
Abuja | $0.85$ | $40.88$ | $0.902$

Figures 2 and 3 show the combined plot and linear fit to the dataset of
coastal and hinterland respectively. For these figures empirical equations Eq.
7 and Eq. 8 were derived for the coastal and hinterland respectively. It will
be observed again from Figures 2 and 3 as well as the derived Eqs. 7 and 8
that the model fit to the hinterland data has better correlation than the fit
to the coastal data.

Figure 2. Coastal Region data plot and linear fit.

$T_{m}=0.60T_{s}+112.73\;\;\;\;\;R^{2}=0.459$ (7)

Figure 3. Hinterland Region data plot and linear fit.

$T_{m}=0.72T_{s}+77.89\;\;\;\;\;R^{2}=0.752$ (8)

Figure 4 shows a plot of all the combined data for West Africa. The
statistical model derived from Figure 4 is given in Eq. 9 with the
corresponding coefficient of determination:

$T_{m}=0.617T_{s}+108.49\;\;\;\;\;R^{2}=0.709$ (9)

Isioye et. al.[10] using radiosonde data covering $2009$ to $2013$ for West
Africa obtained a corresponding statistical model as given in Eq. 10:

$T_{m}=0.5743T_{s}+116.60\;\;\;\;\;R^{2}=0.436$ (10)

It will be noted that Eq. 9 shows a better statistical performance than Eq.
10, though the two statistical models could be observed to be quite similar. A
statistical comparison of Eq. 9 with Eq. 10 yields root-mean-square error
($RMSE$) of $1.189$ and mean-biased error ($MBE$) of $0.0953$, further showing
that the two statistical models are comparable. Satellite data has wider
spatio-temporal coverage than radiosonde data hence, the better coefficient of
determination $R^{2}$, observed for Eq. 9 as compared to Eq. 10.

Figure 4. A combined West African data plot and linear fit.

## IV. Conclusion

Using satellite data from the Satellite Application Facility (CM-SAF) of the
EUMETSAT a suitable statistical model has been derived to estimate weighted
atmospheric mean temperature $T_{m}$. The statistical model derived is a
simple linear model of $T_{m}$ as a function of surface temperature $T_{s}$.

The derived model was compared with a similar statistical model which had
earlier been derived from the radiosonde data in West Africa. The
correspondence between the two models were found to be significant.

Acknowledgment
This work was done using data from EUMETSAT’s Satellite Application Facility
on Climate Monitoring (CM SAF). DOI: $10.5676$/EUM_SAF_CM/WVT_ATOVS/V$001$.

## References

  * [1] Bevis, M., Businger, S., Herring, T. A., Rocken, C., Anthes, R. A. and Ware R. H. (1992).  GPS Meteorology: Remote Sensing of Atmospheric Water Vapor using the Global Positioning System.  Journ. Geophy. Res., 97. doi:10.1029/92JD01517.
  * [2] Raju, C. S., Saha, K., Thampi, B. V. and Parameswaran, K. (2007).  Empirical Model for Mean Temperature for Indian Zone and Estimation of Precipitable Water Vapor from ground based GPS Measurements.  Annal of Geophysics, 25, 1935-1948.
  * [3] Fernandez, W., KhaChatryan, V., Tumasyan, A. (2010).  Observation of Long range, Near-side Angular Correlations in Proton-proton Collisions at the LHC  Journal of High Energy Physics, 2010:91.
  * [4] Musa, T. A., Amir, S., Othman, R., Ses, S., Omar, K., Abdullah, K., Lim, S. and Rizos, C. (2011).  GPS Meteorology in a low Latitude Region: Remote Sensing of Atmospheric Vapor over the Malaysian Peninsula.  Journal of Atmospheric Physics, 73:2410-2422.
  * [5] Abimbola, O. J., O. A. Falaiye and J. Omojola (2017).  Estimation of precipitable water vapor in Nigeria using NIGNET GNSS/GPS, NCEP-DOE Reanalysis II and surface meteorological data.  J. Phys. Sci., 28(2), 19 – 29. doi: 10.21315/jps2017.28.2.2.
  * [6] Mendes, V. B., Prates, G., Santos, L. and Langley, R. B. (2000).  An Evaluation of the Accuracy of Models of the Determination of the Weighted Mean Temperature of the Atmosphere.  Proceedings of ION, National Technical Meeting, Jan. 26-28, 2000,.2.2. Pacific Hotel Disneyland, Anaheim, CA.
  * [7] Saastamoinen, J. (1972).  Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. In Henriksen, W. et al. The use of artificial satellites for geodesy, geophysics monograph series, vol. 15. Washinghton DC:  American Geophysical Union, 247, https://doi.org/10.1029/GM015p0247.
  * [8] Solbrig, P. (2000).  Untersuchungenuber die Nutzungnumerischer Wettermodellezur wasser dampfbestimmungmitHilfe des Global Positioning System.  Diploma Thesis, Institude of Geodesy and Navigation, University FAF Munich,Germany, Feb. 2000.
  * [9] Shoji, Y. (2010).  Accurate Estimation of Precipitable Water Vapor using Ground based GPS Observation Network and its data assimilation into a mesoscale numerical weather prediction model.  PhD Thesis, Kyoto University.
  * [10] Isioye, O. A., L. Combrinck and J. Botai (2016).  Modelling weighted mean temperature in the West African region: implications for GNSS meteorology.  Meteorol. Appl. 23, 614 – 632. doi: 10.1002/met.1584.
  * [11] Schuler, T., Posfay, A., Hein, G. W. and Biberger, R. ($2001$).  A global analysis of the mean atmospheric temperature for GPS water vapor estimation.  Inproc. 14th Int. Tech. Meeting Satell. – C5: 558 Atmospheric Effects, 10NGPS, P.559.