Regression problems typically assume the training data are independent samples, but in many modern applications samples come from individuals connected by a network or some other form of baseline information. In the case that the population is divided into discrete and disjoint classes with no notion of inter-class connection, hierarchical linear modelling can address such issues. We propose a series of optimization schemes that create generalized linear models, which incorporate both the local neighborhoods and global structure of the population network Laplacian to define localized regressions that smoothly vary with respect to the network. Depending on the context of the data and function, this method can be used to learn local function gradients on a manifold, model and denoise time varying processes with drift terms determined by some initial condition on a network, and perform co-clustering of a matrix or tensor. We also discuss incorporating the local regression coefficient estimates to redefine the network geometry and build a function adapted diffusion process.