Most of the non-linear transceivers, which are based on Tomlinson Harashima (TH) precoding and have been proposed in the literature for two-way relay networks, assume perfect channel state information (CSI). In this paper, we propose a novel and robust TH precoding scheme for two-way relay networks with multiple antennas at the transceiver and the relay nodes. We assume imperfect CSI and that the channel uncertainty is bounded by a spherical region. Furthermore, we consider the sum of the mean square error as the objective function, under a limited power constraint for transceiver and relay nodes. Simulations are provided to evaluate the performance and to validate the efficiency of the proposed scheme.