Thanks! That worked.
You can use LinearModelFit:
data = Flatten[Table[{x, y, x + Sin[x + y]}, {x, 5}, {y, 5}], 1]; lmf[x_, y_] = Normal@LinearModelFit[data, {x, y}, {x, y}]; Plot3D[x + Sin[x + y] - lmf[x, y], {x, 0, 5}, {y, 0, 5}]