Starburst galaxies are efficient $\gamma$-ray producers, because their high supernova rates generate copious cosmic ray (CR) protons, and their high gas densities act as thick targets off which these protons can produce neutral pions and thence $\gamma$-rays. In this paper we present a first-principles calculation of the mechanisms by which CRs propagate through such environments, combining astrochemical models with analysis of turbulence in weakly ionised plasma. We show that CRs cannot scatter off the strong large-scale turbulence found in starbursts, because efficient ion-neutral damping prevents such turbulence from cascading down to the scales of CR gyroradii. Instead, CRs stream along field lines at a rate determined by the competition between streaming instability and ion-neutral damping, leading to transport via a process of field line random walk. This results in an effective diffusion coefficient that is nearly energy-independent up to CR energies of $\sim 1$ TeV. We apply our computed diffusion coefficient to a simple model of CR escape and loss, and show that the resulting $\gamma$-ray spectra are in good agreement with the observed spectra of the starbursts NGC 253, M82, and Arp 220. In particular, our model reproduces these galaxies' relatively hard GeV $\gamma$-ray spectra and softer TeV spectra without the need for any fine-tuning of advective escape times or the shape of the CR injection spectrum.