Superconducting transmon qubits exhibit considerable coherence variability. Both the decoherence and its variability are generally attributed in large part to coupling of the qubit to two-level systems (TLSs) in the device materials. Much progress has already been made in understanding TLSs and their role in qubit decoherence. Building on this understanding, we implement a device-scale numerical model of a qubit interacting with the surrounding TLSs. TLSs are distributed according to the Standard Tunneling Model, and their contributions to qubit decoherence are summed incoherently. With judicious approximations based on existing literature, we are able to eliminate all but a few free parameters in the model, and even those are constrained by available data. In this way we generate ensembles of nominally identical qubits that differ only in the statistical selection of the millions of TLSs close enough to interact. The results explain the observed statistical variability of qubits, and quantify the relative importance of the few TLSs that interact most strongly, vs the background of many weak interactions. We also compare with qubit spectroscopy results. Finally, we address time-domain variability of the decoherence.