Superhydrophobic surfaces (SHSs) have the potential to reduce drag at solid boundaries. However, multiple independent studies have recently shown that small amounts of surfactant, naturally present in the environment, can induce Marangoni forces that increase drag, at least in the laminar regime. To obtain accurate drag predictions, one must solve the mass, momentum, bulk surfactant and interfacial surfactant conservation equations. This requires expensive simulations, thus preventing surfactant from being widely considered in SHS studies. To address this issue, we propose a theory for steady, pressure-driven, laminar, two-dimensional flow in a periodic SHS channel with soluble surfactant. We linearise the coupling between flow and surfactant, under the assumption of small concentration, finding a scaling prediction for the local slip length. To obtain the drag reduction and interfacial shear, we find a series solution for the velocity field by assuming Stokes flow in the bulk and uniform interfacial shear. We find how the slip and drag depend on the six dimensionless groups characterizing surfactant transport near SHSs, the gas fraction and the ratio of the interface length to the channel half-height. Our model agrees with numerical simulations spanning orders of magnitude in each dimensionless group. The simulations also provide the three constants in the scaling theory. Our model significantly improves predictions relative to a surfactant-free one, which can otherwise overestimate slip and underestimate drag by several orders of magnitude. Our slip length model can provide the boundary condition in other simulations, thereby accounting for surfactant effects without having to solve the full problem.