Certain disorder-free Hamiltonians can be non-ergodic due to a strong fragmentation of the Hilbert space into disconnected sectors. Here, we characterize such systems by introducing the notion of 'statistically localized integrals of motion' (SLIOM), whose eigenvalues label the connected components of the Hilbert space. SLIOMs are not spatially localized in the operator sense, but appear localized to sub-extensive regions when their expectation value is taken in typical states with a finite density of particles. We illustrate this general concept on several Hamiltonians, both with and without dipole conservation. Furthermore, we demonstrate that there exist perturbations which destroy these integrals of motion in the bulk of the system, while keeping them on the boundary. This results in statistically localized strong zero modes, leading to infinitely long-lived edge magnetizations along with a thermalizing bulk, constituting the first example of such strong edge modes in a non-integrable model. We also show that in a particular example, these edge modes lead to the appearance of topological string order in a certain subset of highly excited eigenstates. Some of our suggested models can be realized in Rydberg quantum simulators.CONTENTS * These authors contributed equally to this work.