We propose a numerical method for time-dependent, semilinear partial differential equations (PDEs) with random parameters and random initial data. The method is based on an operator splitting approach. The linear part of the right-hand side is discretized by a stochastic Galerkin method in the stochastic variables and a pseudospectral method in the physical space, whereas the nonlinear part is approximated by a stochastic collocation method in the stochastic variables. In this setting both parts of the random PDE can be propagated very efficiently. The Galerkin method and the collocation method are combined with sparse grids in order to reduce the computational costs. This approach is discussed in detail for the Lugiato-Lefever equation, which serves as a motivating example throughout, but also applies to a much larger class of random PDEs. For such problems our method is computationally much cheaper than standard stochastic Galerkin methods, and numerical tests show that it outperforms standard stochastic collocation methods, too.