The next step is missing for me:
ops = lm(reshape(x, length(x), 1), y) y_linreg_fitted = coef(ops) .* x