In this paper, a shape optimization problem constrained by a random elliptic partial differential equation with a pure Neumann boundary is presented. The model is motivated by applications in interface identification, where we assume coefficients and inputs are subject to uncertainty. The problem is posed as a minimization of the expectation of a random objective functional depending on the uncertain parameters. A numerical method for iteratively solving the problem is presented, which is a generalization of the classical stochastic gradient method in shape spaces. Moreover, we perform numerical experiments, which demonstrate the effectiveness of the algorithm.