  We present a robust approach to incorporating models for the time-variable
broadening of the pulse profile due to scattering in the ionized interstellar
medium into profile-domain pulsar timing analysis. We use this approach to
simultaneously estimate temporal variations in both the dispersion measure (DM)
and scattering, together with a model for the pulse profile that includes
smooth evolution as a function of frequency, and the pulsar's timing model. We
show that fixing the scattering timescales when forming time-of-arrival
estimates, as has been suggested in the context of traditional pulsar timing
analysis, can significantly underestimate the uncertainties in both DM, and the
arrival time of the pulse, leading to bias in the timing parameters. We apply
our method using a new, publicly available, GPU accelerated code, both to
simulations, and observations of the millisecond pulsar PSR J1643$-$1224. This
pulsar is known to exhibit significant scattering variability compared to
typical millisecond pulsars, and we find including low-frequency ($< 1$ GHz)
data without a model for these scattering variations leads to significant
periodic structure in the DM, and also biases the astrometric parameters at the
$4\sigma$ level, for example, changing proper motion in right ascension by
$0.50 \pm 0.12$. If low frequency observations are to be included when
significant scattering variations are present, we conclude it is necessary to
not just model those variations, but also to sample the parameters that
describe the variations simultaneously with all other parameters in the model,
a task for which profile domain pulsar timing is ideally suited.
