We study the scaling relations between the baryonic content and total mass of groups of galaxies, as these systems provide a unique way to examine the role of non-gravitational processes in structure formation. Using Planck and ROSAT data, we conduct detailed comparisons of the stacked thermal Sunyaev-Zel'dovich (tSZ) and X-ray scaling relations of galaxy groups found in the Galaxy And Mass Assembly (GAMA) survey and the BAHAMAS hydrodynamical simulation. We use weak gravitational lensing data from the Kilo Degree Survey to determine the average halo mass of the studied systems. We analyse the simulation in the same way, using realistic weak lensing, X-ray, and tSZ synthetic observations. Furthermore, to keep selection biases under control, we employ exactly the same galaxy selection and group identification procedures to the observations and simulation. Applying this comparison, we find that the simulations reproduce the richness, size, and stellar mass functions of GAMA groups, as well as the stacked weak lensing and tSZ signals in bins of group stellar mass. However, the simulations predict X-ray luminosities that are higher than observed for this optically selected group sample. As the same simulations were previously shown to match the luminosities of X-ray-selected groups, this suggests that X-ray-selected systems may form a biased subset. Finally, we demonstrate that our observational processing of the X-ray and tSZ signals is free of significant biases. We find that our optical group selection procedure has, however, some room for improvement.