In this article, I present a new multivariate regression model for analyzing outcomes with network dependence. The model is capable to account for two types of outcome dependence including the mean dependence that allows the outcome to depend on selected features of a known dependence network and the error dependence that allows the outcome to be additionally correlated based on patterned connections in the dependence network (e.g., according to whether the ties are asymmetric, mutual, or triadic). For example, when predicting a group of students' smoking status, the outcome can depend on the students' positions in their friendship network and also be correlated among friends. I show that analyses ignoring the mean dependence can lead to severe bias in the estimated coefficients while analyses ignoring the error dependence can lead to inefficient inferences and failures in recognizing unmeasured social processes. I compare the new model with related models such as multilevel models, spatial regression models, and exponential random graph models and show their connections and differences. I propose a two-step, feasible generalized least squares estimator to estimate the model that is computationally fast and robust. Simulations show the validity of the new model (and the estimator) while four empirical examples demonstrate its versatility. Associated R package "fglsnet" is available for public use.