I don't think the fit can be improved unless you assume a different model (like higher order polynomial or whatever). That looks like a linear fit like what you'd get from
coefficients = polyfit(x, y, 1)
a = -coefficients(1);
b = coefficients(2);
and the fit is what it is. You can't make the line get any closer to the points since the line it gives you is already the closest overall to the points that is possible.
Now if you want to go to a higher order, you could do
coefficients = polyfit(x, y, 3)
fittedY = polyval(coefficients, x)
You can see that it's no more complicated than the linear fit you wanted to "keep it simple". So just pick the best one. But you can't pick the best model on only 15 sample points. I'd get data from a few hundred patients/subjects and then look at the scatterplot and see what you think the best model would be.