Summary: In this article, we introduce a new class of nonparametric regression models, called generalized spatially varying coefficient models (GSVCMs), for data distributed over complex domains. For model estimation, we propose a nonparametric quasi-likelihood approach using the bivariate penalized spline approximation technique. We show that our estimation procedure is able to handle irregularly-shaped spatial domains with complex boundaries. Under some regularity conditions, the estimator for the coefficient function is proved to be consistent in the \(L_2\) sense and its convergence rate is established. We develop a numerically stable algorithm using penalized iteratively reweighted least squares method to estimate the coefficient functions in GSVCMs. To gain efficiency in the computation for large-scale data, we further propose a QR decomposition-based algorithm, which requires only sub-blocks of the design matrix to be computed at a time, so that it allows efficient estimation of GSVCMs for large datasets with modest computer hardware. The finite sample performance of the GSVCM and its estimation method is examined by simulations studies. The proposed method is also illustrated by an analysis of the crash data in Florida. Supplementary materials for this article are available online.