  In precision medicine, identifying optimal sequences of decision rules,
termed dynamic treatment regimes (DTRs), is an important undertaking. One
approach investigators may take to infer about optimal DTRs is via Bayesian
dynamic Marginal Structural Models (MSMs). These models represent the expected
outcome under adherence to a DTR for DTRs in a family indexed by a parameter $
\psi $; the function mapping regimes in the family to the expected outcome
under adherence to a DTR is known as the value function. Models that allow for
the straightforward identification of an optimal DTR may lead to biased
estimates. If such a model is computationally tractable, common wisdom says
that a grid-search for the optimal DTR may obviate this difficulty. In a
Bayesian context, computational difficulties may be compounded if a posterior
mean must be calculated at each grid point. We seek to alleviate these
inferential challenges by implementing Gaussian Process ($ \mathcal{GP} $)
optimization methods for estimators for the causal effect of adherence to a
specified DTR. We examine how to identify optimal DTRs in settings where the
value function is multi-modal, which are often not addressed in the DTR
literature. We conclude that a $ \mathcal{GP} $ modeling approach that
acknowledges noise in the estimated response surface leads to improved results.
Additionally, we find that a grid-search may not always yield a robust solution
and that it is often less efficient than a $ \mathcal{GP} $ approach. We
illustrate the use of the proposed methods by analyzing a clinical dataset with
the aim of quantifying the effect of different patterns of HIV therapy.
