This work is concerned with various weighted four dimensional matrix summability methods in modular function spaces associating with generalized difference operator involving (p, q)-gamma function. Following very recent results of Kadak (2017) , we first introduce a new type of difference operator on double function sequences and, also define modularly weighted A-statistical convergence and modularly statistically <(N)(A)over tilde> -summability by means of weighted four dimensional regular matrix A. We also present some important inclusion relations between newly proposed summability methods. Moreover, based on the concept of modularly statistically (N) over tilde (A)-summability, we prove a Korovkin type approximation theorem for functions of two variables in modular spaces. Finally, in order to show that our proposed method and its applications to approximation results are stronger than the existing methods, we display an illustrative example using bivariate (p, q)-Bernstein Kantorovich operators. (C) 2018 Elsevier Inc. All rights reserved.