Symbolic regression (SR) is a function identification process, the task of which is to identify the relationship between two sets of variables and express the relationship in mathematical models. SR is named to emphasize the fact that it has the ability to find the shape of the model and the coefficients in the model simultaneously. As one of the most successful application of Genetic Programming (GP), it has been widely researched in the GP community. GP based SR does not require any predefined model and has a flexible representation. These advantages make GP base SR an attractive and powerful alternative to classical regression techniques. Despite the many successful stories, GP based SR generally has poor generalisation ability which degrades its reliability. Meanwhile, this drawback also hampers the application of SR to science and real-world modeling. Therefore, this research aims to develop a new GP approach to SR that can evolve models exhibiting good generalisation ability. The generalisation will be measured on the same domain as well as over different but similar domains. Another issue is that the available data for regression problems might be high-dimensional. When searching in a large feature space, convergence of GP will be slow and the evolved models often accompany with high risk of overfitting. However, feature selection which is desired for high-dimensional data is seldom studied for regression problems, let alone research for SR. Therefore, this research will develop a novel feature selection method for GP based SR.