We present non-homogeneous Markov regression models of unknown order as a means to assess the duration of autoregressive dependence in longitudinal binary data. We describe a subject's transition probability evolving over time using logistic regression models for his or her past outcomes and covariates. When the initial values of the binary process are unknown, they are treated as latent variables. The unknown initial values, model parameters, and the order of transitions are then estimated using a Bayesian variable selection approach, via Gibbs sampling. As a comparison with our approach, we also implement the deviance information criterion (DIC) for the determination of the order of transitions. An example addresses the progression of substance use in a community sample of n = 242 American Indian children who were interviewed annually four times. An extension of the Markov model to account for subject-to-subject heterogeneity is also discussed.