Fifty years of pulsar data has led to the discovery of emission and rotation variability on timescales of months and years; we have developed techniques to identify this long timescale variability. Individual observations may be too noisy to identify subtle changes in a pulse profile; we use Gaussian process regression to model noisy observations and produce a continuous map of pulse profile variability. Generally, multiple observing epochs are required to obtain the pulsar spin frequency derivative. Gaussian process regression is, therefore, also used to monitor this rate of spindown. We have applied variability detection techniques to both millisecond and long period pulsar datasets. I will discuss the techniques used and present the most interesting results from the pulsars analysed.