The y_pred vector should hold all predictions for your observations present in the X_test dataset, which should be 756 observations. If you want to use your model on the whole dataset you can simply use the .predict() method on your X_train dataset:
# predict on your training dataset
# predict on your test dataset
Seems like Mixture density networks provide a great solution.
Gaussian mixtures can be numerically compared to both the linear regression MSE loss approach and the softmax cross-entropy loss approach via negative log likelihood.
Using lm is not the right approach to model a binary outcome. You would use a Logit in this case (see some example here and see why not lm here).
However, there are (at least) two more issues:
You have a highly unbalanced target
You may have "noisy" features
You should check if some oversampling of the minority class or using SMOTE ...