Therapeutic outcomes in patients with metastatic colorectal cancer (mCRC) receiving bevacizumab treatment are highly variable, and a reliable predictive factor is not available. Progression‐free survival (PFS) and overall survival (OS) were recorded from an observational, prospective study after 5 years of follow‐up, including 46 patients with mCRC receiving bevacizumab treatment. Three vascular endothelial growth factor (VEGF)‐A and two intercellular adhesion molecule‐1 genes polymorphisms, age, gender, weight, dosing scheme, and co‐treatments were collected. Given the relatively small number of events (37 [80%] for the PFS and 26 [57%] for the OS), to study the effect of these covariates on PFS and OS, a covariate analysis was performed using statistical and supervised machine learning techniques, including Cox regression, penalized Cox regression techniques (least absolute shrinkage and selection operator [LASSO], ridge regression, and elastic net), survival trees, and survival forest. The predictive performance of each method was evaluated in bootstrapped samples, using prediction error curves and the area under the curve of the receiver operating characteristic. The LASSO penalized Cox‐regression model showed the best overall performance. Nonlinear mixed effects (NLME) models were developed, and a conventional stepwise covariate search was performed. Then, covariates identified as important by the LASSO model were included in the base NLME models developed for PFS and OS, resulting in improved models as compared to those obtained with the stepwise covariate search. It was shown that having gene polymorphisms in VEGFA (rs699947 and rs1570360) and ICAM1 (rs1799969) are associated with a favorable clinical outcome in patients with mCRC receiving bevacizumab treatment.